Numerically Stable Method for Solving Quadratic Equations
The commonly used formula for the solutions of a quadratic does not
provide for the most accurate computation of both roots when faced with
the limitations of finite precision arithmetic. One of the two roots is found
with lower precision than the other due to round-off when two quantities
of the same sign and similar magnitude are subtracted from one another.
By multiplying
ax2 + bx + c = 0 (1)
by 4a and completing the square one gets
(2ax + b)2 + (4ac − b2) = 0 (2)
and so finds
x = −b ± √
b2 − 4ac
2a
(3)
By instead completing the square in
a + b 1
x
+ c
1
x2 = 0 (4)
one finds
1
x = −b ∓ √
b2 − 4ac
2c
(5)
or
x = 2c
−b ∓ √
b2 − 4ac
(6)
Given limitations of computer arithmetic, one or the other of these (eq. (3)
or eq. (6)) may provide more accuracy for a particular root. When quantities
of the same sign are subtracted, some loss in precision may be expected.
This is a particular concern here if ac is relatively small compared
to b2, in which case b has about the same magnitude as √
b2 − 4ac.
This suggests that one use one of the above equations for one root,
and use the other equation for the other root:
x1 = −b − √
b2 − 4ac
2a
& x2 = 2c
−b − √
b2 − 4ac
when b ≥ 0 (7)
and
x1 = 2c
−b + √
b2 − 4ac
& x2 = −b + √
b2 − 4ac
2a
when b < 0 (8)
We can check these results by noting that x1x2 = c/a and x1+x2 = −b/a.
Note that no more work is involved using eq. (7) or eq. (8) than blindly
using either eq. (3) or eq. (6) for both solutions of the quadratic.