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
>> <[email protected]>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 <[email protected]>
>>> 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 <[email protected]> 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 <[email protected]>
>>>>>> 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 <[email protected]> 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
