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