Re: help with how ML/Octave selects answer for multivalued functions.

[Nicholas Jankowski]

On Fri, Aug 21, 2015 at 5:13 PM, Nicholas Jankowski <address@hidden> wrote:
>
> Looking into how Octave does fractional powers to complex numbers.
> 

<SNIP>


First, apologizes for the multi-post last week. Second, apologies for the horrible typos to anyone who tried to follow the email. my cut-and-paste skills were lacking as I tried to do too many things at once. Just posting the clean up + outcome for posterity, since mailing lists live forever:

starting with a unit vector in C: z= a+bi = r*exp(i*theta)

z^(1/n) is multivalued for n>1, and the values can be found from:

z^(1/n) = r^(1/n)*exp(i*(theta+2*k*pi)/n),   (k = 0,1,...,n-1)

as an example: for r=1, theta = 1, n = 3:

>> theta = 1;n=3;k = 0:floor(n)-1;
>> exp(1i*(theta+2*k*pi)/n)
ans =
   0.94496 + 0.32719i  -0.75584 + 0.65476i  -0.18912 - 0.98195i

if Octave calculates the exponent directly, it chooses the principle, or k=0 value, which coincides with the 'algebraic non-complex simplification':

>> (exp(i*theta))^(1/n)
ans =  0.94496 + 0.32719i

>> exp(i*theta/n)
ans =  0.94496 + 0.32719i

OK... all well and good.

Now, for -pi<theta<pi, Octave sticks with the principle value. Exceeding pi, however:

>> theta =1.01*pi;n=3;k = 0:floor(n)-1;

>> exp(1i*(theta+2*k*pi)/n)
ans =
   0.49090 + 0.87121i  -0.99995 - 0.01047i   0.50904 - 0.86074i

>> (exp(i*theta)^(1/n)
ans =  0.50904 - 0.86074i

>> exp(i*theta/n)
ans =  0.49090 + 0.87121i

If we trust wikipedia, I believe that many codes use A^x = exp(x*log(A)) to calculate powers of complex numbers, which then uses the principle value of the complex logarithm.  In 'xpow.cc' I was pointed to a bit ago, z^x falls under "case 7", but I don't know what file std::pow pulls from (line 358) to verify the operation takes place in Octave. Evidence seems to suggest so.

The error comes from the log function truncating the complex argument to principle value (limiting to [-pi,pi) ), which then creates a different output when multiplied by x. This is unavoidable when complex arguments are stored in rectangular form. a third value is necessary to represent projection in the 3rd dimension (e.g., as a Riemann surface a la https://en.wikipedia.org/wiki/Complex_logarithm#The_associated_Riemann_surface )

in my case, since I have the exact polar form of the complex number with non-truncated argument, I can use this to calculate things directly now that I know how it all works. It doesn't streamline the main problem as nicely, but it works. For anyone interested, I did put an answer up to the original problem over at:
http://math.stackexchange.com/questions/1383443/matrix-with-fractional-exponent-not-getting-expected-output-in-matlab-octave

Cheers!

Nick J.