From: "James Giles" <jamesgiles@worldnet.att.net>
Newsgroups: sci.math.num-analysis,comp.lang.fortran
Subject: Re: Most stable formula for evaluating roots of quadratic equation
Date: 16 Nov 1998 19:38:53 GMT

Jos Bergervoet wrote in message ...
>Michael Buck <buck@ike.uni-stuttgart.de> writes:
>
>>I am looking for the numerically most stable, i.e. most accurate
>>formulation for evaluating the roots
>>of a quadratic equation (2nd order polynom)
>
>Use the code listed below.
...
>      Subroutine quadrc(a,b,c,w1,w2)
>!    -- Solves a x^2 + b x + c = 0; w1 becomes root with smaller abs. val.
>        implicit complex*16(a-z)
>
>        b2 = -b/2
>        d = sqrt(b2*b2-a*c)
>        s = b2+d
>        s2= b2-d
>        if(abs(s2) .gt. abs(s)) s=s2
>        w2 = s/a
>        w1 = 0
>        if(c.ne.0) w1=c/s
>      end

I wasn't going to answer this thread because I figured there would be
many essentially correct solutions posted, but I believe it is always
desirable to give the reasons for the changes you suggest to someone.

First, given that the original statement of the problem was that A and
B were greater than zero and C was less than zero, there is no need
to introduce the COMPLEX data type into his problem.  The expression
within the SQRT call will always be positive.

Second, the introduction of B2 saves a little time (and slightly decreases
the chance that B**2 will overflow) but does nothing to increase the
accuracy of the solution.  I'm not saying that it's a bad idea: the saved
time is no doubt worth something.  But it should be explained that
the modification if for speed.

Third, you use B2*B2 instead of B2**2.  This should be
what the compiler does anyway for small integer exponents.  But
if it isn't, the single multiply is almost certainly more accurate than
the exp(2*log(B2)) that might usually be the alternative.  It's
always a good idea to multiply rather than exponentiate if you
can.

Fourth, the introduction of S and S2 is again mostly a time saver.
The normal solution divides by A to get the first root, and then
inverts that root and multiplies it by C/A to get the second root.
This second division by A essentially produces the S value you
had before.  There is possibly a little loss of precision from this,
but not significant in this case (see below).

Fifth, the IF statement to choose that largest value solution as
the principal root and computing the other from it resolves the
issue of "front cancellation" when subtracting nearly equal
values.  But again, this should be a moot point in this case
(next paragraph).  If he _did_want the second root, this would
be a very important step.  Subtracting two (nearly equal) values
which are both the result of calculation is a good way to _promote_
error.

Given that A and B are positive, C is negative, and he only wanted
the (most) negative root, the naive solution is:

   R1 = (-B - SQRT(B*B - 4*A*C)) / (2*A)

This is also the largest magnitude root of the two, given the prior
knowledge of the values.  In spite of your changing variables,
this is the result your program would produce as well (assuming that
multiplications by 4 and 2 are exact as well as your division by 2).
Since the other root is never needed, the whole discussion above
about S and S2 is not really relevant in this case.  Nor is the
issue of selecting the largest root as the principal one.  (If he decides
he needs the other root, the statement R2 = C/(A*R1) can be used.
This still avoids the "front cancellation" from the subtract.)

Horst Kraemer also gave a solution.  In it, he also did a lot
of substitution of variables.  He defined P as B/(2*A) and
Q as C/A.  The solution was then given for:

   X^2 + 2PX + Q = 0

The solution also involved some time saving and a test for the
largest magnitute root.  But, in the end it should again produce the
same solution as R1 above for the case in question.

In the case, as described, if the solution given by the naive
solution is not adequate, you might try using extended precision
(say double the precision that the variables A, B, and C are
declared) for all the intermediate calculations in the solution.
This would probably be about the best you could do (assuming
the extended precision SQRT function was "good to the last ULP")

If the system you're on doesn't have extended precision (or if the
cost of an extended precision SQRT is prohibitive), you might
try:

   R1_PRIME = R1 - (A*R1*R1 + B*R1 + C) / (2*A*R1 + B)

This is a Newton's method iteration applied to the quadratic formula.
Assuming that the problem really is a poor approximation of R1 from
the naive formula, this should correct the approximation quite precisely
(assuming that R1 is close at all).  It may even be cheaper to do the
naive calculation and the Newton's iteration both in normal precision
than to do a single expression of the naive solution in extended precision.
Applying the various time saving variable substitutions may also help.

Maybe I've confused the issue now, with too much information.  But
I really do think that answers should give explanations and not just
code to blindly follow.  Anyway, good luck.

--
J. Giles