Hey, the expression parser must be working!

Paul F. Dubois pauldubois at home.com
Fri Apr 14 00:02:19 CEST 2000


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)







More information about the Python-list mailing list