January 2002 - PURRS-devel - BUGSENG Mailman

Re: Roots of unity
by Roberto Bagnara 05 Feb '02

05 Feb '02

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(a)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; }

3 4

purrs/doc purrs.tex
by Alessandro Zaccagnini 31 Jan '02

31 Jan '02

CVSROOT: /cvs/purrs Module name: purrs Changes by: zaccagnini(a)cs.unipr.it 2002-01-31 14:52:46 Modified files: doc : purrs.tex Log message: 1. Added comments on the papers by Lueker and by Cohen & Katcoff; 2. Corrected more misprints; 3. Added paragraph on the solution of the equation of degree 4 (notation consistent with code); 4. Added an example of reduction of a 5-variable system of linear recurrences from Cohen & Katcoff into a single equation of order 5; 5. Renamed section 2.2 to 2.1.5 (which is more appropriate) and expanded to include the general solution of a linear recurrence of arbitrary (but finite) order; 6. Changed notation in sections 2.2.3 and 2.2.4 (for the fundamental solution) for consistency with the notation in section 2.2.5; 7. Added reference to continued fractions in general; 8. Removed a comment from section 4.1.3, because it was wrong. Patches: http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/doc/purrs.tex.diff?cvsroot=…

1 0

purrs/src alg_eq_solver.cc
by Roberto Bagnara 31 Jan '02

31 Jan '02

CVSROOT: /cvs/purrs Module name: purrs Changes by: roberto(a)cs.unipr.it 2002-01-31 12:27:59 Modified files: src : alg_eq_solver.cc Log message: Doxygen comments are not allowed within function bodies. Ordinary C++ comments should be used instead. Patches: http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/src/alg_eq_solver.cc.diff?c…

1 0

purrs/src alg_eq_solver.cc
by Alessandro Zaccagnini 31 Jan '02

31 Jan '02

CVSROOT: /cvs/purrs Module name: purrs Changes by: zaccagnini(a)cs.unipr.it 2002-01-31 12:15:28 Modified files: src : alg_eq_solver.cc Log message: Added appropriate comments in the body of solve_equation_3 for the handling of cubic roots so as to avoid taking roots of negative numbers. Patches: http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/src/alg_eq_solver.cc.diff?c…

1 0

purrs/infrastructure README
by Roberto Bagnara 31 Jan '02

31 Jan '02

CVSROOT: /cvs/purrs Module name: purrs Changes by: roberto(a)cs.unipr.it 2002-01-31 07:39:41 Modified files: infrastructure : README Log message: Explain how to build NTL. Patches: http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/infrastructure/README.diff?…

1 0

purrs/src alg_eq_solver.cc
by Roberto Bagnara 30 Jan '02

30 Jan '02

CVSROOT: /cvs/purrs Module name: purrs Changes by: roberto(a)cs.unipr.it 2002-01-30 23:00:08 Modified files: src : alg_eq_solver.cc Log message: Make sure cubic_root() is always invoked on non-negative real numbers. Patches: http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/src/alg_eq_solver.cc.diff?c…

1 0

purrs/src Makefile.am
by Roberto Bagnara 30 Jan '02

30 Jan '02

CVSROOT: /cvs/purrs Module name: purrs Changes by: roberto(a)cs.unipr.it 2002-01-30 22:59:14 Modified files: src : Makefile.am Log message: Disable interval computations. Patches: http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/src/Makefile.am.diff?cvsroo…

1 0

purrs/tests Makefile.am bug1.cc eqs_driver.cc
by Roberto Bagnara 30 Jan '02

30 Jan '02

CVSROOT: /cvs/purrs Module name: purrs Changes by: roberto(a)cs.unipr.it 2002-01-30 22:58:16 Modified files: tests : Makefile.am bug1.cc eqs_driver.cc Log message: Disable interval computations. Patches: http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/tests/Makefile.am.diff?cvsr… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/tests/bug1.cc.diff?cvsroot=… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/tests/eqs_driver.cc.diff?cv…

1 0

purrs infrastructure/README src/Makefile.am sr ...
by Roberto Bagnara 30 Jan '02

30 Jan '02

CVSROOT: /cvs/purrs Module name: purrs Changes by: roberto(a)cs.unipr.it 2002-01-30 19:33:49 Modified files: infrastructure : README src : Makefile.am alg_eq_solver.cc Log message: Do not build the approximation module. Some comments added/improved. Patches: http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/infrastructure/README.diff?… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/src/Makefile.am.diff?cvsroo… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/src/alg_eq_solver.cc.diff?c…

1 0

fi_lib++ (2nd attempt)
by Roberto Bagnara 30 Jan '02

30 Jan '02

Dear Sirs, first of all I would like to thank you for fi_lib++. I have just finished reading the documentation and it turns out that fi_lib++ has almost all the features of the library I was writing myself. I have a few questions: 1) What is the official web site for the library? 2) Is there a mailing list for fi_lib++ users and developers? 3) If not, to whom should I report bugs and improvements? To start with, I am having namespace problems (in my application I need to define FILIB_NAMESPACES). The problems I have found so far are solved by the patch attached below. 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 diff -rcp interval.orig/macro/include/Interval.h interval/macro/include/Interval.h *** interval.orig/macro/include/Interval.h Thu Nov 1 11:15:24 2001 --- interval/macro/include/Interval.h Fri Jan 25 18:52:55 2002 *************** *** 33,39 **** #endif #ifdef FILIB_NAMESPACES ! namesupersetace filib { #endif --- 33,39 ---- #endif #ifdef FILIB_NAMESPACES ! namespace filib { #endif diff -rcp interval.orig/macro/include/rounding.h interval/macro/include/rounding.h *** interval.orig/macro/include/rounding.h Wed Apr 25 15:48:56 2001 --- interval/macro/include/rounding.h Fri Jan 25 18:54:40 2002 *************** *** 22,33 **** #define ROUND_UP #define ROUND_NEAR static const double ZERO_SUCC=Double::succ(0.0); static const double ZERO_PRED=Double::pred(0.0); static const double ONE_SUCC=Double::succ(1.0); ! static const double ONE_PRED=Double::pred(1.0); ! #define low(a) ((a)>0)?((a)*ONE_PRED):((a)=0?ZERO_PRED:(a)*ONE_SUCC) #define high(a) ((a)>0)?((a)*ONE_SUCC):((a)=0?ZERO_SUCC:(a)*ONE_PRED) --- 22,40 ---- #define ROUND_UP #define ROUND_NEAR + #ifndef FILIB_NAMESPACES static const double ZERO_SUCC=Double::succ(0.0); static const double ZERO_PRED=Double::pred(0.0); static const double ONE_SUCC=Double::succ(1.0); ! static const double ONE_PRED=Double::pred(1.0); ! #else ! static const double ZERO_SUCC=filib::Double::succ(0.0); ! static const double ZERO_PRED=filib::Double::pred(0.0); ! ! static const double ONE_SUCC=filib::Double::succ(1.0); ! static const double ONE_PRED=filib::Double::pred(1.0); ! #endif #define low(a) ((a)>0)?((a)*ONE_PRED):((a)=0?ZERO_PRED:(a)*ONE_SUCC) #define high(a) ((a)>0)?((a)*ONE_SUCC):((a)=0?ZERO_SUCC:(a)*ONE_PRED)

2 1