Solving a cubic equation
Question: As part of a program I'm writing, I need to solve a cubic equation exactly
(rather than using a numerical root finder):

    
    
    a*x**3 + b*x**2 + c*x + d = 0.
    

I'm trying to use the equations from
[here](http://en.wikipedia.org/wiki/Cubic_function#Root-finding_formula).
However, consider the following code (this is Python but it's pretty generic
code):

    
    
    a =  1.0
    b =  0.0
    c =  0.2 - 1.0
    d = -0.7 * 0.2
    
    q = (3*a*c - b**2) / (9 * a**2)
    r = (9*a*b*c - 27*a**2*d - 2*b**3) / (54*a**3)
    
    print "q = ",q
    print "r = ",r
    
    delta = q**3 + r**2
    
    print "delta = ",delta
    
    # here delta is less than zero so we use the second set of equations from the article:
    
    rho = (-q**3)**0.5
    
    # For x1 the imaginary part is unimportant since it cancels out
    s_real = rho**(1./3.)
    t_real = rho**(1./3.)
    
    print "s [real] = ",s_real
    print "t [real] = ",t_real
    
    x1 = s_real + t_real - b / (3. * a)
    
    print "x1 = ", x1
    
    print "should be zero: ",a*x1**3+b*x1**2+c*x1+d
    

But the output is:

    
    
    q =  -0.266666666667
    r =  0.07
    delta =  -0.014062962963
    s [real] =  0.516397779494
    t [real] =  0.516397779494
    x1 =  1.03279555899
    should be zero:  0.135412149064
    

so the output is not zero, and so x1 isn't actually a solution. Is there a
mistake in the Wikipedia article?

ps: I know that numpy.roots will solve this kind of equation but I need to do
this for millions of equations and so I need to implement this to work on
arrays of coefficients.


Answer: Wikipedia's notation `(rho^(1/3), theta/3)` does not mean that `rho^(1/3)` is
the real part and `theta/3` is the imaginary part. Rather, this is in polar
coordinates. Thus, if you want the real part, you would take `rho^(1/3) *
cos(theta/3)`.

I made these changes to your code and it worked for me:

    
    
    theta = arccos(r/rho)
    s_real = rho**(1./3.) * cos( theta/3)
    t_real = rho**(1./3.) * cos(-theta/3)
    

(Of course, `s_real = t_real` here because `cos` is even.)

