Aleksandar Donev wrote:

> Program received signal SIGFPE, Arithmetic exception.
> 661 collision_time = current_time + C /
> (-B+Sqrt(B**2-A*C))

and then in another letter wrote:

> I've read somewhere that what I wrote was the better way, as opposed
> to the usual formula (-B+Sqrt(B**2-A*C))/A (my A and C have a factor
> of 2 inside), but I don't remember where I read this and what the
> explanation was...I would appreciate if anyone does know of such
> better approaches.

The two aren't equivalent. They select opposite roots. To select the
same root, and to see the correct alternative mathematical formulation,
multiply the last term by 1 == (-B-Sqrt(B**2-A*C))/(-B-Sqrt(B**2-A*C)) to
get (-B-Sqrt(B**2-A*C))/A. (Notice the "-" in front of the Sqrt.)

This form and your original one are for one root of A*x**2 + 2*B*x + C ==
0. The other root has the opposite sign on the Sqrt.

The form you're using is the wrong one for the root you're using if B >
0. It guarantees cancellation, and therefore loss of precision: As
abs(A*C) becomes ever smaller compared to B**2, A*C is subtracted from
later and later digits of B**2, and eventually from digits that aren't
represented. So B**2-A*C gets closer and closer to B**2. This isn't a
problem when B < 0, but when B > 0, more and more of the leading digits
of -B+Sqrt(B**2-A*C) cancel. When the leading digits cancel, the
low-order digits are filled with zeros, which aren't necessarily correct,
so you have poor relative accuracy. Eventually you get a computational
if not a mathematical zero. If B < 0, -B+Sqrt(B**2-A*C) -> -2*B as A*C
shrinks, and therefore has no cancellation, but -B-Sqrt(B**2-A*C) has the
problems discussed here. So you shouldn't use the same formulation in
both the B > 0 and B < 0 cases.

When B > 0 you should use (-B-Sqrt(B**2-A*C))/A. When B < 0, you should
use the form you have. Neither has cancellation, no matter the size of
A*C. When B < 0, as A*C gets small compared to B**2, yours approaches
-C/(2*B); when B > 0 the one here approaches -2*B/A. If you want the
other root, the reasoning about the sign of B is reversed.

When B > 0, in the limit as A approaches zero the root you're selecting
goes off to the infinity with the sign of -A, but the other root is
-C/(2*B). When B < 0 the situation is reversed. If A == 0.0 can occur
in your physical problem, you should have a special branch for it. You
should also wonder whether you're selecting the correct root in the
general case. Or maybe you need to check both roots in the A /= 0 case.

If A==0 .and. B==0, the problem and the roots aren't defined. Otherwise
if C == 0 .and. (A==0 .or. B==0) the root(s) are zero.

Otherwise if B == 0.0 can occur frequently, provide a special branch for
it, where the solutions are +/- SQRT(-C/A). If C == 0.0 can occur
frequently, provide a special branch for it, where one solution is
-2.0*B/A (the other one is 0.0).

This illustrates that even for the "simple" problem of calculating the
roots of a quadratic equation, a carefully written library-quality
procedure isn't the obvious student exercise one might suspect it to be.

--
Van Snyder | What fraction of Americans believe
[log in to unmask] | Wrestling is real and NASA is fake?
Any alleged opinions are my own and have not been approved or disapproved
by JPL, CalTech, NASA, Sean O'Keefe, George Bush, the Pope, or anybody else.