# Roots of unity

Roberto Bagnara bagnara at cs.unipr.it
Fri Jan 18 10:19:16 CET 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
>
> 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?

> 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
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;
}

```