Roots of unity

18 Jan 2002

      ginsh:
    > evalf((-1)^(1/3));
    0.5+0.86602540378443864673*I

While this answer is technically correct:
    exp(I*pi/3) = -1

It's not what I was expecting, and the cause of some headache today.  I was
trying to write a little routine like RootOf to solve a cubic polynomial, and
was getting three complex answers.  (at least one of them must be real if you
start with a real polynomial)

How should I go about enforcing that (-1)^(1/3) = -1?  I'm interested in
polynomials over the reals, so the complex solutions don't interest me.  This
behaviour of evalf is preventing me from extracting the real solution.

Here's my little function if anyone is interested:

// Find the roots of a polynomial.
// The polynomial can be at most 4th order.
lst RootOf(const ex& poly, const ex& x) {
    ex p = poly;
    ex D,Q,R,S,T;
    p.expand();
    ex a0 = poly.coeff(x,0);
    ex a1 = poly.coeff(x,1);
    ex a2 = poly.coeff(x,2);
    ex a3 = poly.coeff(x,3);
    ex a4 = poly.coeff(x,4);
    switch(p.degree(x)) {
        case 0:
            return lst();
        case 1:
            return lst(-a0/a1);
        case 2:
            return lst((-a1 + sqrt(a1*a1-4*a2*a0))/(2*a2), (-a1 - sqrt(a1*a1-4*a2*a0))/(2*a2));
        case 3:
            a0 /= a3;
            a1 /= a3;
            a2 /= a3;
            Q = (3*a1-a2*a2)/9;
            R = (9*a2*a1-27*a0-2*a2*a2*a2)/54;
            D = pow(Q,3) + pow(R,2);
            S = pow(R+sqrt(D),numeric(1,3));
            T = pow(R-sqrt(D),numeric(1,3));
            // These formulas come from: http://mathworld.wolfram.com/CubicEquation.html
            return lst(
                    -numeric(1,3)*a2+S+T,
                    -numeric(1,3)*a2-numeric(1,2)*(S+T)+numeric(1,2)*I*sqrt(ex(3))*(S-T),
                    -numeric(1,3)*a2-numeric(1,2)*(S+T)-numeric(1,2)*I*sqrt(ex(3))*(S-T)
            );
        //case 4:
            // http://mathworld.wolfram.com/QuarticEquation.html
            // It seems this is defined ambiguously. eq. 30 wants "a real root" but there
            // may be as many as 3!  Can there 12 roots!??!
        default:
            throw(std::domain_error("RootOf: Cannot take the root of a polynomial with degree > 4"));
    }
    return lst();
}

Cheers,
-- Bob

Bob McElrath (rsmcelrath@students.wisc.edu) 
Univ. of Wisconsin at Madison, Department of Physics