[Jprogramming] Cool roots
John Randall
randall at andromeda.rutgers.edu
Wed Dec 12 17:30:56 UTC 2012
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.canada at 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 <HenryHRich at 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 <km at 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.canada at 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 <km at 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
More information about the Programming
mailing list