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
> 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;
}
More information about the GiNaC-list
mailing list