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
Hello there,
On Fri, 18 Jan 2002, Roberto Bagnara wrote:
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?
There are obviously a lot of choices that can be made in this field and they were the cause of much agony. The question is, what should be the guidelines for these choices? Instead of choosing `visual simplicity', which could be interpreted as a vote for (-1)^(1/3) => -1, we chose to accept `branch cut consistency', which is more rigorously definable.
The roots (even or odd) of negative numbers are always the first ones counting counter-clock-wise from the positive real axis.
Why is this? There are several guidelines. The usual convention is to define the n'th root of x==r*exp(I*p), r and p being real, like such:
x^(1/n) = exp(log(x)/n) == exp((log(r) + I*p)/n) == exp(log(r)/n) + exp(I*p/n)
The definition of the principal branch of the log function is such that the branch cut runs along the negative real axis, continuous with quadrant II, i.e. log(-1) == +I*Pi. Since the exp function does not feature a branch cut, this defines the branch cut of all the roots as well.
This definition is arbitrary from a mathematical point of few. However, for practical purposes, we need to fix the branch cuts somehow and we did not come up with our own definitions. For instance, the Common Lisp standard mandates them. I can warmly recommend Guy L. Steele's book "Common Lisp, the Language, 2nd Edition", which we should think of as the de-facto standard for Lisp. Section 12.5 specifies these things and is an excellent treatment. Just search Google for it, there are numerous copies of it on the web. (But skip the parts dealing with signed zeros, because we don't have these.)
Why, could one ask next, should we bother about Lisp? Because they did dare to specify these things and people have been following their recommendation. I checked a lot of critical cases and found that both MapleV and Mathematica follow the same convention with respect to branch cuts. Also, we studied the ISO C++ standard and while this standard unfortunately leaves some room for interpretation there is hope in sight: The revised C standard (C99) specifies complex numbers and the choice of all branch cuts is exactly identical to the one in Common Lisp. (One just has to read carefully and find the amusing section 7.3.3.2 saying "implementations shall map a cut so the function is continuous as the cut is approached coming around the finite endpoint of the cut in a counter clockwise direction.") The above consistency between Lisp, C, Maple and Mathematica (and GiNaC) also holds true for all inverse hyperbolic and inverse trigonometric functions like tanh, sinh, atanh, etc.
Unfortunately, other Computer Algebra Systems (notably MuPAD and Reduce) chose their own, mutually incompatible, definitions for branch cuts.
If somebody finds a discrepancy, then we should fix it -- I am not guaranteeing that we got every case correct (but we've tried quite hard). I am only saying that there are some external constraints and thinking that the choice is entirely ours would be foolish.
Well, I hope this explains your observation with regards to evalf().
Oh, and (-1)^(1/2) gets eval'ed while (-1)^(1/3) does not because the first one has an exact solution while the latter does not. That is all behind it.
Hope this helps -richy.
Roberto Bagnara [bagnara@cs.unipr.it] wrote:
Bob McElrath wrote:
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.
[...]
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"?
No, their choice is a reasonable one, but it makes wrong answers for this problem.
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:
[...]
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.
Thank you for your code, it is useful but I want to solve a more general problem.
I want to find the roots of a cubic equation _algebraically_. That is, I can't check if the discriminant is positive or negative because there may be undetermined constants involved. Then, given the _algebraic_ solution(s) to the cubic equation, I should be able to plug it in to the original polynomial and get zero. This behavior of evalf is giving me cube-roots which are all complex and none are actually roots.
Richard B. Kreckel [kreckel@zino.physik.uni-mainz.de] wrote:
Hello there,
[...]
There are obviously a lot of choices that can be made in this field and they were the cause of much agony. The question is, what should be the guidelines for these choices? Instead of choosing `visual simplicity', which could be interpreted as a vote for (-1)^(1/3) => -1, we chose to accept `branch cut consistency', which is more rigorously definable.
I am not questioning the choice of taking the roots of unity as you have. I was surprised, but I checked and Maple does the same thing. Rather I just want to find the roots of a cubic equation, and have those roots be correct. ;)
Examining the way Maple does it...it always seems to return to you something where the argument of the cube root is entirely positive. However when I ask it for the solution with undetermined coefficents:
myroots := solve(x^3+a2*x^2+a1*x-a0,x);
1/3 1/3 1/2 1/3 myroots := 1/6 %1 - 6 %2 - 1/3 a2, - 1/12 %1 + 3 %2 - 1/3 a2 + 1/2 I 3 (1/6 %1 + 6 %2),
1/3 1/2 1/3 - 1/12 %1 + 3 %2 - 1/3 a2 - 1/2 I 3 (1/6 %1 + 6 %2)
3 3 2 2 2 3 1/2 %1 := 36 a1 a2 + 108 a0 - 8 a2 + 12 (12 a1 - 3 a1 a2 + 54 a1 a2 a0 + 81 a0 - 12 a0 a2 )
2 1/3 a1 - 1/9 a2 %2 := ---------------- 1/3 %1
Looking at %1 there, it is the argument to the cube root. If it is negative that should demonstrate the roots of unity problem. The following choice is sufficent:
discr := 36*a1*a2+108*a0-8*a2^3+12*sqrt(12*a1^3-3*a1^2*a2^2+54*a1*a2*a0+81*a0^2-12*a0*a2^3);
3 3 2 2 2 3 1/2 discr := 36 a1 a2 + 108 a0 - 8 a2 + 12 (12 a1 - 3 a1 a2 + 54 a1 a2 a0 + 81 a0 - 12 a0 a2 )
evalf(subs(a2=11,a1=5,a0=7,discr));
-7912. + 3691.243693 I
However, if we look at the root itself:
evalf(subs({a2=11,a1=5,a0=7},myroots[1]));
.590815406
Therefore, Maple is not applying the branch cut definition of roots of unity to the roots of this polynomial. I wish to do a similar thing.
I think this should be possible. Any suggestions on how to go about doing it?
I was thinking along the lines of Maple's assume() functionality. Assuming assume() were introduced to ginac, what assumptions would you place on the subexpressions in the cube root to make the answers come out correctly?
Could the power class be modified to include information about which root-of-unity is the apropriate one? How does that affect nonrational exponents?
Cheers, -- Bob
Bob McElrath (rsmcelrath@students.wisc.edu) Univ. of Wisconsin at Madison, Department of Physics
Hi,
On Fri, 18 Jan 2002, Roberto Bagnara wrote: [...]
For the most up-to-date information see the PURRS site: http://www.cs.unipr.it/purrs/ . */
I just wanted to have a look and found that this site is password-protected. Is this intentional?
Regards -richy.
Richard B. Kreckel wrote:
Hi,
On Fri, 18 Jan 2002, Roberto Bagnara wrote: [...]
For the most up-to-date information see the PURRS site: http://www.cs.unipr.it/purrs/ . */
I just wanted to have a look and found that this site is password-protected. Is this intentional?
Hi Richard,
yes, the protection is meant to be there. The reason is simply that the PURRS web site is, for the time being, almost identical to the PPL web site (http://www.cs.unipr.it/ppl/) with the string "PPL" replaced by the string "PURRS". In our sources the web site is already mentioned, but release time is months away.
We will open the web site well before the initial release of
the code, but we must provide some contents before that ;-) Ciao
Roberto
participants (3)
-
Bob McElrath -
Richard B. Kreckel -
Roberto Bagnara