Roots of unity

[Roberto Bagnara]

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 at cs.unipr.it


/* Definitions for the algebraic equation solver.
   Copyright (C) 2002 Roberto Bagnara <bagnara at 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;
}