[cfp-interest 3570] Re: [Cfp-interest 3332] special cases of cpow

[Paul Zimmermann]

thank you Jim for your answer, and sorry for the late reply.

It took us some time to think more about it.

What we decided for the power function in GNU MPC, when both x and y have
imaginary part +/-0, and the imaginary part of x^y should be zero:

* when replacing each imaginary part with tiny values of same sign and
  different magnitude, if we always get the same sign for Im(x^y), we
  return a zero of that sign for Im(x^y)

* if we get different signs for Im(x^y), we return +0 for Im(x^y)

This can be decided from the values and signs of the real parts of x and y.
We have documented this in the "algorithms.pdf" file coming along GNU MPC.

The same rule of course applies when the real part of x^y should be zero.
Thus in the following example we return (+0, 2), which matches your reasoning.

Anyway, if some standard one day specifies this case, we might follow it.

Best regards,
Paul

> From: Jim Thomas <jaswthomas at sbcglobal.net>
> Date: Tue, 3 Dec 2024 08:27:37 -0800
> Cc: CFP <cfp-interest at ucbtest.org>,
>  andreas.enge at inria.fr
> 
> Paul, 
> 
> Here’s a try at deducing an answer to your question using mathematical functions extended to match C and IEEE 754 special cases. As you note, Annex G defines cpow(z, w) to have a branch cut for z along the negative real axis. This fits with clog(z) which is defined to have a branch cut along the negative real axis, and the range of the strip [-pi i, +pi i] along the imaginary axis. Evaluating z ^ w = (-4 + 0 i) ^ (1/2 + 0 i) as e ^ (w ln(z)) ...
> 
> ln(-4 + 0i) = ln(4) + pi i
> 
> (1/2 + 0 i) (ln(4) + pi i) = ln(4)/2 + pi/2 i
> 
> Using the Annex F (and IEEE 754) specification that cospi(n + 1/2) returns +0 for integers n ...
> 
> e^(ln(4)/2 + pi/2 i) = e^(ln(4)/2) cis(pi/2) = 2 (+0 + i) = +0 + 2 i
> 
> 
> - Jim Thomas
> 
> 
> 
> > On Nov 21, 2024, at 6:42 AM, Paul Zimmermann <Paul.Zimmermann at inria.fr> wrote:
> > 
> >       Hi,
> > 
> > with Andreas Enge (in cc), co-author of GNU MPC, we wonder about the special cases of
> > the complex power, in particular when one of the arguments has a zero part (+0 or -0).
> > 
> > For the real power, IEEE 754 completely specifies the special cases of pow.
> > 
> > Is there a similar document for the complex power? N3301 only says "The cpow functions compute
> > the complex power function xy, with a branch cut for the first parameter along the negative real
> > axis." and "This allows cpow(z, c) to be implemented as cexp(c*clog(z)) without precluding
> > implementations that treat special cases more carefully".
> > 
> > For example, what should (-4,+0) ^ (0.5, +0) return? (+0, 2) or (-0,2)?
> > If we replace each +0 by a tiny value, we get different signs for the real part of the
> > result depending on the ratio of these tiny values:
> > 
> > sage: e=1e-6
> > sage: (-4+e*i)^(0.5+e^2*i)
> > 2.49997227485645e-7 + 1.99999999999373*I
> > sage: (-4+e^2*i)^(0.5+e*i)
> > -2.77257976169088e-6 + 1.99999371682264*I
> > 
> > Paul
> > 
> > 
> > _______________________________________________
> > Cfp-interest mailing list
> > Cfp-interest at oakapple.net
> > http://mailman.oakapple.net/mailman/listinfo/cfp-interest
>