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:

1. $b^2 \gg \left\| 4ac \right\|$ which we’ll fix with an algebraic rewrite
2. $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:

$\text{sgn}\left(x\right) = \begin{cases} 1 & x \geq 0 \\[2ex] -1 & x < 0 \end{cases}$

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$)