I reinvented the Durand-Kerner method (
http://en.wikipedia.org/wiki/Durand-Kerner_method ) in 1992. I think it cool. I
don't know if it is better or worse than p. .
- Bo
>________________________________
> Fra: John Randall <rand...@andromeda.rutgers.edu>
>Til: programm...@jsoftware.com
>Sendt: 18:30 onsdag den 12. december 2012
>Emne: Re: [Jprogramming] Cool roots
>
>As Henry points out, the quadratic formula is not the best way to get
>accurate numerical results in all cases.
>
>The two easiest ways of losing numerical precision are dividing by a small
>quantity (which happens with the quadratic formula when a is small) and
>subtracting nearly equal quantities (which happens if b^2 is large
>compared to 4ac). Here are some examples of the latter.
>
>We define 3 verbs which estimate the solutions of the quadratic equation
>with a=1, c=1, b=-y+%y, which has solutions y,%y. The verbs return the
>relative error (|(actual-estimate)%actual) in the estimates.
>
>Warning: these verbs are not cool.
>
>pdot: uses p.
>traditional: uses the quadratic formula as is
>modified: uses the "rationalizing the numerator" technique described by Henry
>
>pdot=:3 : 0
>b=.-+/s=.y,%y
>t=.1{:: p. 1,b,1
>|(s-t)%s
>)
>
>traditional=:3 : 0
>b=.-+/s=.y,%y
>d=.%: _4+*: b
>t=.-:-b+(,~ -) d
>|(s-t)%s
>)
>
>modified=:3 : 0
>b=.-+/s=.y,%y
>d=.%: _4+*: b
>t=.(,%)-:d-b
>|(s-t)%s
>)
>
> pdot"0 ] 10^6+i.6
> 0 0
>1.00583e_14 1.00585e_14
>1.49012e_16 1.65436e_16
> 0 0
> 0 0
> 0 1.61559e_16
> traditional"0 ] 10^6+i.6
>0 7.61449e_6
>0 0.00348485
>0 0.490116
>0 1
>0 1
>0 1
> modified"0 ] 10^6+i.6
>0 0
>0 0
>0 0
>0 0
>0 0
>0 0
>
>As we would expect, there is special magic in pdot. The traditional
>formula goes wrong in expected ways, and the modified one stays good.
>
>Incidentally, the quadratic formula assumes a is nonzero. Describing
>solutions of equations where a can be zero is more complicated. Even when
>solving the linear equation ax=b you have the possibilities:
>
>a~:0 : x=b%a
>a=0, b=0 : infinitely many
>a=0, b~:0 : no solutions
>
>Best wishes,
>
>John
>
>
>Henry Rich wrote:
>> I have to disagree. That's way cool.
>>
>> I assumed you used Newton's method to shine up the final values of the
>> roots, but I never knew you tried rational results too.
>>
>> It is possible to produce accurate results for quadratics without such
>> heroic measures. You just can't use the quadratic formula like you
>> learned it in highschool.
>>
>> There's the usual quadratic formula:
>>
>> (2*a) %~ (-b) (+,-) %: (b^2) - 4*a*c
>>
>> and there's a less-well-known alternative form:
>>
>> (2*c) % (-b) (-,+) %: (b^2) - 4*a*c
>>
>> Each version produces two roots. The roots created are more accurate
>> when +-(-b) is positive, so you can get one accurate root from each
>> formula.
>>
>> To calculate, use the combined form
>>
>> ((a %~ ]) , (c % ])) _0.5 * (%: (*:b) - 4*a*c) ((* *) + ]) b
>>
>> Accuracy has a coolness all its own.
>>
>> Henry Rich
>>
>> On 12/12/2012 3:13 AM, Roger Hui wrote:
>>>> ... produce more accurate results in some difficult cases.
>>>
>>> c=: p. <20$1.5
>>> c
>>> 3325.26 _44336.8 280799 _1.1232e6 3.18239e6 _6.78911e6 1.13152e7
>>> _1.50869e7
>>> 1.63441e7 _1.45281e7 1.0654e7 _6.45695e6 3.22847e6 _1.3245e6 441501
>>> _117734
>>> 24527.8 _3847.5 427.5 _30 1
>>> p. c
>>> ┌─┬───────────────────────────────────────────────────────────────────────────────┐
>>> │1│3r2 3r2 3r2 3r2 3r2 3r2 3r2 3r2 3r2 3r2 3r2 3r2 3r2 3r2 3r2 3r2
>>> 3r2 3r2
>>> 3r2 3r2│
>>> └─┴───────────────────────────────────────────────────────────────────────────────┘
>>>
>>>
>>>
>>> On Tue, Dec 11, 2012 at 11:29 PM, Roger Hui
>>> <rogerhui.can...@gmail.com>wrote:
>>>
>>>> You may or may not know that p. employs some extraordinary measures
>>>> which
>>>> produce more accurate results in some difficult cases. But those
>>>> extraordinary measures are not "cool". For example:
>>>>
>>>> w=: p. <1+i.20 NB. Wilkinson's
>>>> polynomial<http://en.wikipedia.org/wiki/Wilkinson_polynomial>
>>>> w
>>>> 2432902008176640000 _8752948036761600000 13803759753640704000
>>>> _12870931245150988800 8037811822645051776 _3599979517947607200
>>>> 1206647803780373360 _311333643161390640 63030812099294896
>>>> _10142299865511450 1307535010540395 _135585182899530 11310276995381
>>>> _7561...
>>>>
>>>> p. w
>>>> ┌─┬──────────────────────────────────────────────────┐
>>>> │1│20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1│
>>>> └─┴──────────────────────────────────────────────────┘
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> On Tue, Dec 11, 2012 at 3:40 PM, Henry Rich <henryhr...@nc.rr.com>
>>>> wrote:
>>>>
>>>>> None of these cute ways is very accurate in tough cases:
>>>>>
>>>>> 0j15 ": (2*a) %~ (-b) (+,-) %: (b^2) - 4*a*c [ 'a b c' =. 1e_6 1e6
>>>>> 1e_6
>>>>> 0.000000000000000 _1000000000000.000000000000000
>>>>>
>>>>> But p. does better:
>>>>>
>>>>> 0j15 ": 1 {:: p. c,b,a
>>>>> _1000000000000.000000000000000 _0.000000000001000
>>>>>
>>>>> Henry Rich
>>>>>
>>>>>
>>>>>
>>>>> On 12/11/2012 6:29 PM, Roger Hui wrote:
>>>>>
>>>>>> There are some cheeky (or is it cheesy?) versions:
>>>>>>
>>>>>> (2*a) %~ (-b) (+,-) %: (b^2) - 4*a*c NB. Kip Murray
>>>>>> (2*a) %~ - b (+,-) %: (b^2) - 4*a*c
>>>>>> (+:a) %~ - b (+,-) %: (*:b) - 4*a*c
>>>>>> -: a %~ - b (+,-) %: (*:b) - 4*a*c
>>>>>>
>>>>>>
>>>>>>
>>>>>> On Tue, Dec 11, 2012 at 11:37 AM, km <k...@math.uh.edu> wrote:
>>>>>>
>>>>>> It appears this could be translated into J as the rather cool
>>>>>>>
>>>>>>> (2*a) %~ (-b) (+,-) %: (b^2) - 4*a*c
>>>>>>>
>>>>>>> Sent from my iPad
>>>>>>>
>>>>>>>
>>>>>>> On Dec 11, 2012, at 12:59 PM, Roger Hui <rogerhui.can...@gmail.com>
>>>>>>> wrote:
>>>>>>>
>>>>>>> Example from the Iverson and McDonnell *Phrasal
>>>>>>>> Forms*<http://www.jsoftware.**com/papers/fork.htm<http://www.jsoftware.com/papers/fork.htm>>paper
>>>>>>>> (which
>>>>>>>> introduced fork):
>>>>>>>>
>>>>>>>> (-b)(+,-)√((b*2)-4×a×c)÷2×a
>>>>>>>>
>>>>>>>> √ is a postulated APL primitive, spelled %: in J.
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> On Tue, Dec 11, 2012 at 10:49 AM, km <k...@math.uh.edu> wrote:
>>>>>>>>
>>>>>>>> What is the coolest way of programming the quadratic formula in
>>>>>>>> J? We
>>>>>>>>>
>>>>>>>> are
>>>>>>>
>>>>>>>> finding the roots of polynomial c + x*(b + x*a) without using p. .
>>>>>>>> I
>>>>>>>>>
>>>>>>>> offer
>>>>>>>
>>>>>>>>
>>>>>>>>> roots
>>>>>>>>> 3 : 0
>>>>>>>>> 'a b c' =. y
>>>>>>>>> q =. %: (b^2) - 4*a*c
>>>>>>>>> (2*a) %~ (-b) + q,-q
>>>>>>>>> )
>>>>>>>>> roots 1 3 2
>>>>>>>>> _1 _2
>>>>>>>>> roots 1 0 1
>>>>>>>>> 0j1 0j_1
>>>>>>>>> roots 1 _2 1
>>>>>>>>> 1 1
>>>>>>>>>
>>>>>>>>> partly as problem definition. I am looking for cool roots verbs!
>>>>>>>>>
>>>>>>>>> Kip Murray
>>>>>>>>>
>>>>>>>>> Sent from my iPad
>>>>>>>>>
>>>>>>>>> ------------------------------**------------------------------**
>>>>>>>>> ----------
>>>>>>>>> For information about J forums see http://www.jsoftware.com/**
>>>>>>>>> forums.htm <http://www.jsoftware.com/forums.htm>
>>>>>>>>>
>>>>>>>> ------------------------------**------------------------------**
>>>>>>>> ----------
>>>>>>>> For information about J forums see http://www.jsoftware.com/**
>>>>>>>> forums.htm <http://www.jsoftware.com/forums.htm>
>>>>>>>>
>>>>>>> ------------------------------**------------------------------**
>>>>>>> ----------
>>>>>>> For information about J forums see http://www.jsoftware.com/**
>>>>>>> forums.htm <http://www.jsoftware.com/forums.htm>
>>>>>>>
>>>>>>> ------------------------------**------------------------------**
>>>>>> ----------
>>>>>> For information about J forums see
>>>>>> http://www.jsoftware.com/**forums.htm<http://www.jsoftware.com/forums.htm>
>>>>>>
>>>>>> ------------------------------**------------------------------**
>>>>> ----------
>>>>> For information about J forums see
>>>>> http://www.jsoftware.com/**forums.htm<http://www.jsoftware.com/forums.htm>
>>>>>
>>>>
>>>>
>>> ----------------------------------------------------------------------
>>> For information about J forums see http://www.jsoftware.com/forums.htm
>>>
>> ----------------------------------------------------------------------
>> For information about J forums see http://www.jsoftware.com/forums.htm
>
>
>----------------------------------------------------------------------
>For information about J forums see http://www.jsoftware.com/forums.htm
>
>
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
