Roots of unity

[Bob McElrath]

Roberto Bagnara [bagnara at 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 at 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 at students.wisc.edu) 
Univ. of Wisconsin at Madison, Department of Physics
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 240 bytes
Desc: not available
Url : http://www.cebix.net/pipermail/ginac-list/attachments/20020118/1419ab11/attachment.pgp