Precision errors and quadratic equation

Solving quadratic equation is a quite common question in mathematics. But because of the floating-point precision issue, for same special coefficients, it will produce numerical errors.

Precision errors and quadratic equation

Let's consider a quadratic equation, $x^2 + bx + c = 0$ and a Python function which receives $b$ and $c$, the coefficients of a monic quadratic polynomial, $x^2 + b x + c$ and returns the pair of its roots.

For example, given a quadratic $x^2 - 2x + 1$, the function should return a pair of $(1, 1)$. Of course, in floating point, the answers may differ from an exact unity.

import cmath
def solve_quad(b, c):
    root = cmath.sqrt(b**2 - 4*c)
    return ((-b-root)/2, (-b+root)/2)

This function should be very straightforward based on the quadratic formula. And it should be able to solve most of the quadratic equations.

Let's test this function on several examples. Based on Vieta's formulas, consider $x_1$ and $x_2$ are two roots of th quadratic formula. We have $x_1x_2 = c$ and $x_1 + x_2 = -b$.

from numpy import allclose
variants = [{'b': 4.0, 'c': 3.0},
            {'b': 3.0, 'c': 1.0},
            {'b': 2, 'c': 4.0},
            {'b': 1e15, 'c': 1e10},
            {'b': 1e15, 'c': 1e30},
            {'b': 1e15, 'c': 1e14},
            {'b': 1e15, 'c': 3.0},
            {'b': -1e10, 'c': 4.0}]
for var in variants:
    x1, x2 = solve_quad(**var)
    print(allclose(x1*x2, var['c']))
True
True
True
False
True
False
False
False

Notice that when $b$ is very small, this function always returns a correct result. But, when $b$ becomes bigger, the result becomes unstable.

Let's review the quadratic formula again to understand what has happened.

We need to evaluate the square root of $b^2 - 4c$. But if $b^2$ is much bigger than $4c$, we will have a precision issue in the square root term. Thus, we notice in the results, when $b = 1e+15, c = 1e+30$, the function outputs a correct result.

Because of the precision error, $-b + \sqrt{b^2 - 4c}$ becomes 0, which outputs the wrong result of the equation. Let's take a closer took of the square root term.

Since $b^2$ is much greater than $4c$, we have

$$\sqrt{b^2 - 4c} = \lim_{\epsilon \to 0^+}{|b| - \epsilon}$$

Since we always have, $b^2 > 4c$, thus, $\epsilon \to 0$. Thus

$$x_1 = -b + \frac{\epsilon}{2}, x_2 = -\frac{\epsilon}{2},\text{for}\ b > 0$$

$$x_1 = -b - \frac{\epsilon}{2}, x_2 = \frac{\epsilon}{2},\ \text{for}\ b< 0$$

Moreover,

$$\frac{-b - \sqrt{b^2 - 4c}}{2} = \frac{(-b - \sqrt{b^2 - 4c})(-b + \sqrt{b^2 - 4c})}{2(-b + \sqrt{b^2 - 4c})} = \frac{2c}{-b + \sqrt{b^2 - 4c}}$$

$$\frac{-b + \sqrt{b^2 - 4c}}{2} = \frac{(-b + \sqrt{b^2 - 4c})(-b - \sqrt{b^2 - 4c})}{2(-b - \sqrt{b^2 - 4c})} = \frac{2c}{-b - \sqrt{b^2 - 4c}}$$

Aagain, we substitute $\sqrt{b^2 - 4c}$ by $\lim_{\epsilon \to 0^+}{|b| - \epsilon}$, we will have

$$x_1 = \frac{2c}{-2b+\epsilon}, x_2 = -\frac{2c}{\epsilon}, \text{for}\ b > 0$$

$$x_1 = \frac{2c}{-2b-\epsilon}, x_2 = \frac{2c}{\epsilon}, \text{for}\ b < 0$$

Let's combine these two equations together, then,

$$x_1 = -b+\frac{\epsilon}{2} = -\frac{2c}{\epsilon}, x_2 = \frac{2c}{-2b+\epsilon} = -\frac{\epsilon}{2}, \text{for}\ b > 0$$

$$x_1 = -b-\frac{\epsilon}{2} = \frac{2c}{\epsilon}, x_2 = \frac{2c}{-2b-\epsilon} = \frac{\epsilon}{2}, \text{for}\ b < 0$$

Since $\epsilon \to 0^+$, we have $lim_{\epsilon \to 0^+}{\frac{2c}{\epsilon}} = \infty$ and $lim_{\epsilon \to 0^+}{\frac{\epsilon}{2} = 0}$. These two limits are unbounded, which introduces precision erros into the result. Fortunately, we can use the other formula, which gives more precise result.

def solve_quad_precise(b, c):
    root = cmath.sqrt(b**2 - 4*c)
    if b > 0:
        x1 = 2*c/(-b-root)
        x2 = (-b-root)/2
    else:
        x1 = 2*c/(-b+root)
        x2= (-b+root)/2
    return x1, x2
for var in variants:
    x1, x2 = solve_quad_precise(**var)
    print(allclose(x1*x2, var['c']))
True
True
True
True
True
True
True
True