[PURRS-devel] Re: Roots of unity

[Richard B. Kreckel]

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.
-- 
Richard B. Kreckel
<Richard.Kreckel at Uni-Mainz.DE>
<http://wwwthep.physik.uni-mainz.de/~kreckel/>