January
10th,
2020

Micro-post: Solve the quadratic equation without catastrophic cancellation using FMA and no branching.

Given $a$, $b$ and $c$ we want to find one or both roots $r_0$ and $r_1$:

\[ax^2 + bx + c = a\left(\left(x-r_0\right)\left(x-r_1\right)\right)\]

Recall we have two sources of catastrophic cancellation:

- $b^2 \gg \left\| 4ac \right\| $ which we’ll fix with an algebraic rewrite
- $b^2 \approx 4ac $ which we’ll fix by computing in extended precision.

Calling $r_1$ the larger and $r_0$ the smaller magnitude root we can rewrite the standard equation as follows:

\[\begin{align*} r_1 & = \frac{-\left(b + \text{sgn}\left(b\right)\sqrt{b^2-4ac}\right) }{2a} \\ & = \frac{b + \text{sgn}\left(b\right)\sqrt{b^2-4ac}}{-2a} \\ \\ r_0 & = \frac{2c}{-\left(b + \text{sgn}\left(b\right)\sqrt{b^2-4ac}\right) } \\ & = \frac{-2c}{b + \text{sgn}\left(b\right)\sqrt{b^2-4ac}} \end{align*}\]

where $\text{sgn}$ is the sign function:

The remaining part is computing $b^2-4ac$ in extended precision which is a special case of extended precision $ab+cd$. This is covered in detail in *“Further analysis of Kahan’s algorithm for the accurate computation of 2x2 determinants”*

Toy code which assumes two real roots (including $r_0=r_1$)