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