Just in case anybody was nervous about whether Python arithmetic expressions
have been adequately tested, the following works. Cool, huh? (Composed using
Maple 6 via Fortran codegeneration and a text editor)
from math import sqrt
def cmplx (x, y):
return x + 1.0j * y
def cubic_roots (a, b, c, d):
"Roots of a*x**3+b*x**2+c*x+d"
t0 = 1/a*(36*c*b*a-108*d*a**2-8*b**3+12*sqrt(3.0)*sqrt(4*c**3*a-c\
**2*b**2-18*c*b*a*d+27*d**2*a**2+4*d*b**3)*a)**(1.0/3.0)/6-2.0/\
3.0*(3*c*a-b**2)/a/(36*c*b*a-108*d*a**2-8*b**3+12*sqrt(3.0)*sqrt\
(4*c**3*a-c**2*b**2-18*c*b*a*d+27*d**2*a**2+4*d*b**3)*a)**(1.0/3.0)-b/a/3
t1 = -1/a*(36*c*b*a-108*d*a**2-8*b**3+12*sqrt(3.0)*sqrt(4*c**3*a-\
c**2*b**2-18*c*b*a*d+27*d**2*a**2+4*d*b**3)*a)**(1.0/3.0)/12+(3*\
c*a-b**2)/a/(36*c*b*a-108*d*a**2-8*b**3+12*sqrt(3.0)*sqrt(4*c**3*\
a-c**2*b**2-18*c*b*a*d+27*d**2*a**2+4*d*b**3)*a)**(1.0/3.0)/3-b/\
a/3+cmplx(0.0,1.0/2.0)*sqrt(3.0)*(1/a*(36*c*b*a-108*d*a**2-8*b\
**3+12*sqrt(3.0)*sqrt(4*c**3*a-c**2*b**2-18*c*b*a*d+27*d**2*a**2+\
4*d*b**3)*a)**(1.0/3.0)/6+2.0/3.0*(3*c*a-b**2)/a/(36*c*b*a-108\
*d*a**2-8*b**3+12*sqrt(3.0)*sqrt(4*c**3*a-c**2*b**2-18*c*b*a*d+27\
*d**2*a**2+4*d*b**3)*a)**(1.0/3.0))
t2 = -1/a*(36*c*b*a-108*d*a**2-8*b**3+12*sqrt(3.0)*sqrt(4*c**3*a-\
c**2*b**2-18*c*b*a*d+27*d**2*a**2+4*d*b**3)*a)**(1.0/3.0)/12+(3*\
c*a-b**2)/a/(36*c*b*a-108*d*a**2-8*b**3+12*sqrt(3.0)*sqrt(4*c**3*\
a-c**2*b**2-18*c*b*a*d+27*d**2*a**2+4*d*b**3)*a)**(1.0/3.0)/3-b/\
a/3+cmplx(0.0,-1.0/2.0)*sqrt(3.0)*(1/a*(36*c*b*a-108*d*a**2-8*\
b**3+12*sqrt(3.0)*sqrt(4*c**3*a-c**2*b**2-18*c*b*a*d+27*d**2*a**2\
+4*d*b**3)*a)**(1.0/3.0)/6+2.0/3.0*(3*c*a-b**2)/a/(36*c*b*a-108\
*d*a**2-8*b**3+12*sqrt(3.0)*sqrt(4*c**3*a-c**2*b**2-18*c*b*a*d+27\
*d**2*a**2+4*d*b**3)*a)**(1.0/3.0))
return [t0, t1, t2]
if __name__ == "__main__":
print cubic_roots(1.0, -1.0, 1.0, -1.0)
