mathtalk-ga, I tried to boil down your answer at
<http://answers.google.com/answers/threadview?id=433886> to something
I thought I could use to write a Python script:

x^3 + a x^2 + b x + c = 0 (I assume a,b,c to be real)

M  =  b - a^2/3

N  =  -(2(a^3)/27 - ab/3 + c)

z is a root of z^3 + Mz = N  if and only if
  
z = w - M/(3w) and w^3 is a root of (w^3)^2 - N(w^3) - (M^3)/27 = 0

I will need to compute the 3 cube roots of w^3 (which I can now do),
and then solve for z. There will be 3 identical pairs of z's.
Discarding one member of each pair, what remain will be the 3 roots of
the cubic equation.

Is my understanding correct? If not, please correct it.

Thanks,

donkan

---

Request for Question Clarification by mathtalk-ga on 12 Dec 2004 07:43 PST

Hi, donkan-ga:

You did leave out the part where I said:

  x is a root of x^3 + a x^2 + b x + c = 0
  
  if and only if
  
  x = z - a/3 and z is a root of z^3 + Mz = N

Maybe for the sake of your Python program, it makes sense to describe
the outer and inner problems as subroutines.

The outer routine takes three "real" numbers a,b,c and (we'll assume)
prints three "complex" roots of x^3 + a x^2 + b x + c = 0.

To do so requires a call to an inner routine, which takes two "real"
numbers M,N and returns three "complex" roots of z^3 + Mz = N.

All the outer routine needs to do is compute M,N (from a,b,c) as
inputs to the inner routine, and then subtract a/3 from the roots
returned by the inner routine to produce the roots for the "original"
problem.

If this is the sort of clarification you need, I'll post as an Answer
detailed comments on how to write the inner routine, which solves z^3
+ Mz = N.

regards, mathtalk-ga

---

Clarification of Question by donkan-ga on 12 Dec 2004 07:56 PST

mathtalk-ga,

Yes, please.

donkan

---

Clarification of Question by donkan-ga on 12 Dec 2004 08:37 PST

mathtalk,

You didn't mention w^3. Don't I have to compute w at some point?

donkan

Hi, donkan-ga:

Below is an implementation of the analytical formulas for solving the
cubic equation, using three functions:

  solveCubic (module name)
    solveCubicMonic(a,b,c)
    solveCubicReduced(m,n)
    cubeRoots(x)

Each function returns a triple of complex numbers.  Take a look and
see if you have any questions.  Since Python treats "whitespace" as
significant in some ways, look carefully to see that the posting here
in this textbox hasn't "broken" the code in some way (I'll check it
also).

regards, mathtalk-ga

>>>>>>>>>>> begin solveCubic.py >>>>>>>>>>>>>>>>>>

# -*- coding: UTF-8 -*-

import math, cmath

def solveCubicMonic(a,b,c):
    """Solve a monic real cubic polynomial.
    
    Input three real coefficients a,b,c and return
    three complex roots of Uł + aU˛ + bU + c = 0.
    """
    m = b - (a*a)/3
    n = (b + m + m)*a/9 - c
    x,y,z = solveCubicReduced(m,n)
    x,y,z = x - a/3, y - a/3, z - a/3
    return (x,y,z)

def solveCubicReduced(m,n):
    """Solve a reduced monic real cubic polynomial.

    Input two real coefficients m,n and return
    three complex roots of uł + mu = n.
    """
    if m == 0:
        return cubeRoots(n)
    if n == 0:
        return (0,sqrt(-m + 0j),-sqrt(-m + 0j))
    m = m/3 + 0j
    n = n/2
    w = n + cmath.sqrt(n*n + m*m*m)
    x,y,z = cubeRoots(w)
    return (x - m/x, y - m/y, z - m/z)

def cubeRoots(x):
    """Find three complex cube roots of real or complex x."""
    z = x + 0j
    s = 1
    if z.real < 0:
        z,s = -z,-s
    t = math.atan2(z.imag, z.real)/3
    a = s * (abs(z)**(1./3))
    w = complex(-0.5,math.sqrt(0.75))
    w_ = w.conjugate()
    u = a * complex(math.cos(t),math.sin(t))
    return (u, u*w, u*w_)

---

Clarification of Answer by mathtalk-ga on 17 Dec 2004 06:13 PST

Okay, I found one mistake.  In the middle of the program:

        return (0,sqrt(-m + 0j),-sqrt(-m + 0j))

should qualify function sqrt with module name cmath:

        return (0,cmath.sqrt(-m + 0j),-cmath.sqrt(-m + 0j))

as it is a couple of lines below in the same function.

regards, mathtalk-ga

Wonderful, mathtalk-ga. You appear to have learned a lot of Python to do this.
I'd like to give you 10 stars to go along with my tip. 

Thank you!!

BTW When I ran what you gave me, IDLE suggested a different coding,   
         "# -*- coding: cp1252 -*-".
You can see your work at <http://www.rcblue.com/Python/solveCubicEquation.py>

donkan

The interested reader may want to check out these links:

[python -- the home page]
http://www.python.org/

The Python programming language was created at CWI in Amsterdam in
1990, and has become a very popular interpretive platform for rapid
prototyping because of its object-oriented and extensible design. 
"Python" is a humorous homage to Monty Python.  The most recent major
release is Python 2.4 (Nov. 30, 2004), which is provided for Windows
platforms starting with this release as an .msi (Microsoft installer)
file.

[Dive Into Python]
http://diveintopython.org/

This site hosts an online and up-to-date copy of the book Dive Into
Python by Mark Pilgrim.  The text includes recommendations on choosing
a Python implementation.

[ActivePython -- freeware Python]
http://www.activestate.com/Products/ActivePython/

Windows users may be interested in Mark Hammond's PythonWin IDE, which
is bundled with the Windows version available here (and at
SourceForge).

[A Brief Introduction to Python for Scientific Computing]
http://volsung.gotdns.org/gem-website/python_intro.html

A concise introductory overview with links to related sites, from the
point of view of using Python for scientific computing/numerical
methods.

[Python Complex Number Arithmetic]
http://faculty.ed.umuc.edu/~swalsh/Python/PCNA.html

As an interpreted language it is often unnecessary to declare
variables, but this raises the question of how the interpreter should
know whether real or complex arithmetic is meant.  An example here
illustrates writing complex numbers using "j" as the imaginary square
root of -1 to implicitly signal that a complex value is intended, i.e.
4 + 0j.  Here "j" functions as a modifier rather than as an
independent name, sort of like using ".0" to force an integer to be
treated as a floating point value.

[5.8 cmath -- Mathematical functions for complex numbers]
http://docs.python.org/lib/module-cmath.html

The "base" module cmath in Python does not have "convenience" function
cbrt for complex arguments, though it does implement sqrt for complex
numbers (and base module math implements cbrt for real ones).  I'll
discuss a workaround in the Answer posted above.

regards, mathtalk-ga

And here's a Python function, written by Tim Peters, for computing the
3 cube roots of a complex number:
<http://www.rcblue.com/Python/croots.txt>

donkan