[PURRS-devel] Re: Roots of unity

18 Jan 2002

      Bob McElrath wrote:
...
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 think it would be unreasonable to require that (-1)^(1/3) = -1,
since the other solutions are as good as -1.
Thus I believe GiNaC does the right thing by leaving the symbolic
expression (-1)^(1/3) untouched.
On the other hand, a little experiment with ginsh shows that
...
1^(1/2);
1
so that here GiNaC does choose one root.  I believed it does
so on the grounds that even roots of non-negative real numbers
are usually defined that way.  However, still with
ginsh, we get also
...
(-1)^(1/2);
I
so it seems this tradition is pushed a little too far.
I think these issues should be clarified with the GiNaC team.
Perhaps you are criticizing the behavior of evalf() when applied to
(-1)^(1/3) and I agree with you.  Since odd roots of real
numbers always have a real solution, why not defining evalf() so that
it does "the Right Thing"?
The point is perhaps: how much would it cost?
GiNaC people?
Please forgive me if I failed to understand your point.
...
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();
}
If that is your problem, I believe the piece of code below will solve it.
Beware: it is extracted from a library we are writing so that it only
received very little testing.  I hope it helps.
All the best
Roberto
-- 
Prof. Roberto Bagnara
Computer Science Group
Department of Mathematics, University of Parma, Italy
http://www.cs.unipr.it/~bagnara/
mailto:bagnara@cs.unipr.it

/* Definitions for the algebraic equation solver.
   Copyright (C) 2002 Roberto Bagnara bagnara@cs.unipr.it

This file is part of the Parma University's Recurrence Relation
Solver (PURRS).

The PURRS is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.

The PURRS is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
USA.

For the most up-to-date information see the PURRS site:
http://www.cs.unipr.it/purrs/ . */

#include <ginac/ginac.h>
#include <cln/complex.h>

typedef GiNaC::ex GExpr;
typedef GiNaC::symbol GSymbol;
typedef GiNaC::lst GList;
typedef GiNaC::numeric GNumber;

using namespace GiNaC;

static GExpr
cubic_root(const GExpr& e) {
  static GExpr one_third = GExpr(1)/3;
  return pow(e, one_third);
}

/*!
  Solve the equation \f$x^3 + a_1 x^2 + a_2 x + a_3 = 0\f$
  and return <CODE>true</CODE> if and only if all the solutions are real.
  \f$x_1\f$ is guaranteed to be a real solution.
*/
bool
solve_equation_3(const GNumber& a1, const GNumber& a2, const GNumber& a3,
    	 GExpr& x1, GExpr& x2, GExpr& x3) {
  GExpr Q = (3*a2 - pow(a1, 2)) / 9;
  GExpr R = (9*a1*a2 - 27*a3 -2*pow(a1, 3)) / 54;
  GExpr d = pow(Q, 3) + pow(R, 2);
  GExpr a1_div_3 = a1/3;
  if (d < 0) {
    GExpr sqrt_minus_Q = sqrt(-Q);
    GExpr theta = acos(-R/Q*sqrt_minus_Q);
    x1 = -a1_div_3 + 2*sqrt_minus_Q*cos(theta/3);
    x2 = -a1_div_3 + 2*sqrt_minus_Q*cos((theta+2*Pi)/3);
    x3 = -a1_div_3 + 2*sqrt_minus_Q*cos((theta+4*Pi)/3);
  }
  else {
    GExpr sqrt_d = sqrt(d);
    GExpr S = cubic_root(R + sqrt_d);
    GExpr T = cubic_root(R - sqrt_d);
    GExpr S_plus_T = S + T;
    GExpr t1 = -S_plus_T/2 - a1_div_3;
    GExpr t2 = (S - T)*I*sqrt(GExpr(3))/2;
    x1 = S_plus_T - a1_div_3;
    x2 = t1 + t2;
    x3 = t1 - t2;
  }
  // The roots are all real if and only if d <= 0.
  return d <= 0;
}