Quadratic Formula

We all remember the quadratic formula from high school,  x = ( -b ± √ b² - 4ac ) / 2a, but if it's so familiar to us all of us, why do I keep coming across such bad implementations in code?

Polynomials

Lets start with a slightly more general problem, we're interested in finding all the real roots of these polynomials in x :

k₀  +  k₁ x  =  0

k₀  +  k₁ x  +  k₂ x²  =  0

k₀  +  k₁ x  +  k₂ x²  +  k₃ x³  =  0

k₀  +  k₁ x  +  k₂ x²  +  k₃ x³  +  k₄ x⁴   =  0

k₀  +  k₁ x  +  k₂ x²  +  k₃ x³  +  k₄ x⁴  +  k₅ x⁵  =  0

We know from the Fundamental Theorem Of Algebra that, for example, a quintic will have at most 5 real roots, so lets go ahead and write out the function prototypes :

(We also know from Galois Theory that quintics might require a numerical approximation, but that's a post for another day!)

The Quadratic Formula, Revisited

descriminant = kLinear² - 4*kQuadratic*kConstant

x = ( -kLinear ± √ descriminant ) / 2 / kQuadratic

So the first case that I hardly ever see handled, is when kQuadratic vanishes, leading to a divide-by-zero exception.  Actually, this case is very easy to handle if we have a family of functions:

Another important case is when kConstant is zero.  We can see graphically that the polynomial goes through the origin :

kConstant is zero means that 0.0 is a root.

Here's the code:

Note the sneaky 1+ and +1

Catastrophic Cancellation

The last problem we need to fix is when the square root of the descriminant is close in magnitude to kLinear.  Graphically :

This happens surprisingly often in physics calculations.

To fix this, consider the product of the roots :

rootProduct = root[0] * root[1]

rootProduct =  (-kLinear + √descriminant)/2/kQuadratic *

(-kLinear - √descriminant)/2/kQuadratic

rootProduct = ( -kLinear + √ descriminant ) * ( -kLinear - √ descriminant ) / 4 / kQuadratic²

Now something cool happens! Remember that (x+y) * (x-y) = x² - y²

rootProduct = ( (-kLinear)² - √ descriminant² ) / 4 / kQuadratic²

rootProduct = ( kLinear² - descriminant ) / 4 / kQuadratic²

Expanding:

rootProduct = ( kLinear² - kLinear² + 4*kQuadratic*kConstant) / 4 / kQuadratic²

rootProduct = ( kQuadratic*kConstant) / kQuadratic²

rootProduct = kConstant / kQuadratic

What this means is that we can use the larger root to calculate the smaller root by dividing by their product, and keep all the accuracy.

Appendix - Quadratic formula in C++

If you're just here for the code, the following two functions are hereby licensed under CC0. 

---

#include <math.h>
#define EPSILON (1e-10)

// From http://missingbytes.blogspot.com/2012/04/quadratic-formula.html
int SolveLinear(double rootDest[1],double kConstant,double kLinear){
    if(-EPSILON<kLinear&&kLinear<EPSILON){
        return 0;
    }
    rootDest[0]=-kConstant/kLinear;
    return 1;
}

// From http://missingbytes.blogspot.com/2012/04/quadratic-formula.html
int SolveQuadratic(double rootDest[2],double kConstant,double kLinear,double kQuadratic){
    if(-EPSILON<kQuadratic&&kQuadratic<EPSILON){
        return SolveLinear(rootDest,kConstant,kLinear);
    }
    if(-EPSILON<kConstant&&kConstant<EPSILON){
        rootDest[0]=0.0f;
        return 1+SolveLinear(rootDest+1,kLinear,kQuadratic);
    }
    double descriminant=kLinear*kLinear-4*kQuadratic*kConstant;
    if(descriminant<0){
        return 0;
    }
    if(descriminant==0){// no EPSILON check!!
        rootDest[0]=-kLinear/(2*kQuadratic);
        return 1;
    }
    double adder=sqrt(descriminant);
    if(adder*kLinear>0){
        adder=-adder;
    }
    double rootLarge=(-kLinear+adder)/(2*kQuadratic);
    rootDest[0]=rootLarge;
    rootDest[1]=kConstant/(kQuadratic*rootLarge);
    //std::sort(rootDest,rootDest+2);
    return 2;
}

---

Even More Math

They tell me that each mathematical formula in your blog post exactly halves your potential readership.

  How does that make you feel?

Why not talk about your feelings in the comments below?