Numerically stable algorithm for solving the quadratic equation when $a$ is very small or $0$

Solving $a x^2 + bx +c=0$ for $x$ gives

$$x = \frac{-b \pm \sqrt{b^2-4ac}}{2a} \text{, for } a \ne 0$$

But for $a = 0$ we get

$$x=-\frac{c}{b}$$

How to implement a numerically stable algorithm for computing $x$ from $a,b,c$ that allows $a$ to be close to zero or zero?

• Look at cs.berkeley.edu/~wkahan/Qdrtcs.pdf. Numerics aside, the exact solution can be extremely sensitive to small variations. Take $a=1, b=0$ and $c<0$ close to zero. Then the solutions are $\pm \sqrt{c}$ which is very sensitive for small $c$. Jul 13 '14 at 21:56

The previous answers correctly identify that there are two quadratic formulas, and that each has a different range of numeric stability; however, they miss the other subtle instabilities.

To recap: the standard quadratic formula, \begin{align*} x=\frac{-b\pm\sqrt{b^2-4\,a\,c}}{2\,a}, \end{align*} is not unique. A second formula can be generated by multiplying the numerator and denominator by the "conjugate" of the numerator, i.e., \begin{align*} x&=\frac{-b\pm\sqrt{b^2-4\,a\,c}}{2\,a}\times\frac{-b\mp\sqrt{b^2-4\,a\,c}}{-b\mp\sqrt{b^2-4\,a\,c}}= \frac{2\,c}{-b\mp\sqrt{b^2-4\,a\,c}}, \end{align*} where $\mp=-(\pm)$. These two expressions are analytically equivalent; however, each is numerically unstable for certain values of the coefficients.

The instability discussed by the other posts is if the term, $4\,a\,c$, is small compared to $b^2$ and the sign of the radical term and $b$ are the same, then catastropic cancellation occurs.

Fortunately, the two quadratic formulas have opposite signs on the radical term for the same roots, thus it is possible to avoid catastrophic cancellation by selecting the stable form. I.e.,

\begin{align*} x_1 &= \frac{-b-{\rm sign}(b)\sqrt{b^2-4\,a\,c}}{2a}& x_2 = \frac{c}{a\,x_1} \end{align*}

For completeness, there are at least three other numerical instabilities in this formula:

1. When $a=0$: thus the equation is linear and has at most one root given by $-c/b$. It should be noted that if $b=0$ as well, the equation is actually $c=0$ and leaves $x$ undefined/unconstrained!
2. When $a\neq0$ and $c=0$: thus one root is at $x=0$. In this case, the second form of the quadratic equation will yield a NaN ($0/0$) for the second root. It is best therefore to determine the second root by factoring out the root at zero to give $-b/a$.
3. Overflow: the floating point representation might over/under-flow during calculation and this is most critical for the $b$ coefficient. Defining $\mathcal{R}_{max}$ as the largest representable floating point number ($\approx1.8\times10^{308}$ for double precision) then $b^2$ will overflow to infinity at $|b|>\sqrt{\mathcal{R}_{max}}$ according to the IEEE 754 specification (which is "only" $\approx1.304\times10^{154}$). You might not define this as an instability, more an out of range error, but I think the reduced upper range of $b$ is surprising and certainly caught me out!
4. ??: I'm pretty sure I've missed some more, floating point arithmetic can be very subtle.

I have written a compile-time algebra library and you can take a look at my implementation of the quadratic formula. Below it is a "more-stable" implementation of the cubic formula, shamelessly stolen from D. Herbison-Evans. For higher-order polynomials I rely on bisection techniques based around Budan's theorem which isn't nearly as difficult to make stable!

• This is an excellent answer that answered many of my questions. Thanks! May 21 '18 at 19:52
• This is a very good answer, and the code was useful. I used it for making my Rust implementation more stable. Let me add if I may a few tips. First, the division by $a$ might overflow. Thus, instead of checking for $0$, I do the division and check finiteness (this will be a theme). Similarly, while $b^2$ is the most likely reason for arg to overflow, it is not the only one. Thus, I also calculate that and check for finiteness. Finally, there's a chance $c/ax_1$ could fail, so I check that too. Sep 3 '19 at 22:39

Sometimes the Citardauq Formula $$x=\frac{2c}{-b\mp \sqrt{b^2-4ac}}$$ can be useful. One chooses the sign that gives no cancellation in the denominator.

• That's a nice way to generalize a quadratic to the $a=0$ case. Do you know of a similar version for cubics? Jul 13 '14 at 23:43
• No, in any case cubics are hopelessly unpleasant by "formula." Though I once made some consulting money by using the Vieta trigonometric solution in an honest to goodness application. Jul 13 '14 at 23:54

Note that only for one from $x_1=\frac{-b-\sqrt{b^2-4ac}}{2a}$ and $x_2\frac{-b+\sqrt{b^2-4ac}}{2a}$ you do "close subtraction", so you can calculate one of them and the second using Viete'a furmula: $x_1 x_2=\frac{c}{a}$.

• Could you elaborate, why only one of those roots causes a problem? Mar 30 '16 at 12:53
• subtraction of values that are close to each other destroys numeric accuracy. if $b$ is positive and $b^2 >> 4ac$ then $x_2$ subtracts two numbers very close to each other; if $b$ is negative and $b^2 >> 4ac$ then $x_1$ subtracts two numbers very close to each other. May 18 '16 at 17:43
• This is not an answer to my question, where I am looking for a solution that avoids to divide by a. Apr 6 '18 at 16:19
• True, but it’s a good answer to half the question, namely when a ≈ 0. I, too, need a formula (or code) which goes smoothly through a=0, but I can see this being useful on another day. Upvoted (to be clear, I also upvoted the right answer). Sep 26 '20 at 9:46