The standard quadratic formula, $x = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a}$, is a cornerstone of algebra. In pure mathematics, it is infallible. In computational practice, however, it harbours a subtle but critical flaw, especially when implemented using finite-precision arithmetic like IEEE 754 floating-point numbers.

The problem emerges from subtraction, specifically when two nearly equal numbers are subtracted. This phenomenon, known as catastrophic cancellation , can lead to a drastic loss of significant digits, polluting the result with noise from rounding errors.

Consider the quadratic equation where the term $4ac$ is much smaller in magnitude than $b^2$. In this scenario, the discriminant $\Delta = b^2 - 4ac$ is very close to $b^2$, making $\sqrt{\Delta}$ very close to $|b|$.

This approximation creates a problem for one of the roots. If $b > 0$, the calculation of the root $x_1 = \frac{-b + \sqrt{b^2 - 4ac}}{2a}$ involves subtracting two nearly identical numbers in the numerator, triggering catastrophic cancellation. The same issue occurs for the other root, $x_2$, if $b < 0$.

Let's analyse the condition number for the roots of $ax^2+bx+c=0$. The condition number of a problem measures how much the output value can change for a small change in the input data. A problem with a high condition number is called ill-conditioned. For the roots $x$, the condition number with respect to changes in $c$ is given by:

When $4ac$ is small, the discriminant $\sqrt{b^2-4ac} \approx |b|$. One root, let's call it the large root, will be approximately $-b/a$. The other root, the small root, will be approximately $-c/b$. For the small root, the condition number becomes large, indicating that the problem of finding this root is inherently sensitive to small errors in the coefficients. But the failure of the standard formula is an issue of numerical instability, not just ill-conditioning.

To avoid this loss of precision, we can use a different formulation for the problematic root. The key is to rationalise the numerator, transforming the problematic subtraction into an addition.

We know that for a quadratic equation, the product of the roots $x_1 x_2$ is equal to $c/a$. This relationship, known as Vieta's formulas, provides the basis for a stable algorithm. The procedure is as follows:

1. First, compute the discriminant $\Delta = b^2 - 4ac$.
2. Compute the first root using the standard formula, but choose the sign inside the numerator to avoid subtraction. This means if $b$ is positive, we use $-b - \sqrt{\Delta}$. If $b$ is negative, we use $-b + \sqrt{\Delta}$. This can be written compactly using the sign function, $\text{sgn}(b)$.

This calculation for $x_1$ is always stable as it involves adding quantities of the same sign (or subtracting quantities of opposite signs).

3. Finally, compute the second root using the product-of-roots relationship to avoid cancellation altogether.

This two-step method ensures that both roots are computed with high relative accuracy, sidestepping the pitfalls of catastrophic cancellation that plague the naive implementation of the quadratic formula.

Let's test your understanding of these numerical pitfalls.

Understanding these stability issues is crucial for writing robust numerical code, ensuring that theoretical algorithms perform reliably in the real world of finite-precision machines.