title: Solving a cubic equation
body: <p>As part of a program I'm writing, I need to solve a cubic equation exactly (rather than using a numerical root finder):</p>
  
  <pre><code>a*x**3 + b*x**2 + c*x + d = 0.
  </code></pre>
  
  <p>I'm trying to use the equations from <a href="https://en.wikipedia.org/wiki/Cubic_equation" rel="nofollow noreferrer">here</a>. However, consider the following code (this is Python but it's pretty generic code):</p>
  
  <pre><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
  </code></pre>
  
  <p>But the output is:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>so the output is not zero, and so x1 isn't actually a solution. Is there a mistake in the Wikipedia article?</p>
  
  <p>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.</p>