Solving a cubic equation
I've looked at the Wikipedia article and your program.
I also solved the equation using Wolfram Alpha and the results there don't match what you get.
I'd just go through your program at each step, use a lot of print statements, and get each intermediate result. Then go through with a calculator and do it yourself.
I can't find what's happening, but where your hand calculations and the program diverge is a good place to look.
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.)