commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Phil Steitz <phil.ste...@gmail.com>
Subject Re: [math] Complex nth roots
Date Sat, 03 Jan 2009 15:44:07 GMT
Ted Dunning wrote:
> I just check on how sbcl handles this.  The result was mixed, but
> instructive as predicted.
>
> Positive and negative infinity behave as expected.
>
> ** (/ 1 0.0)
> #.SB-EXT:SINGLE-FLOAT-POSITIVE-INFINITY
>
> * (/ -1 0.0)
> *
> ***#.SB-EXT:SINGLE-FLOAT-NEGATIVE-INFINITY
>
> *
> Square roots of negative infinity work as for option 3.
>
> ** (sqrt (/ -1 .0))
> #C(0.0 #.SB-EXT:SINGLE-FLOAT-POSITIVE-INFINITY)
>
> *
> Fourth roots are reasonable as well.
> ** (sqrt (sqrt (/ -1 0.0)))
> #C(#.SB-EXT:SINGLE-FLOAT-POSITIVE-INFINITY
> #.SB-EXT:SINGLE-FLOAT-POSITIVE-INFINITY)
> *
> But, the exponentiation operator is inconsistent with these results and
> behaves like option (1) (without documentation as far as I can tell).
>
> ** (expt (/ -1 .0) 0.25)
> #.SB-EXT:SINGLE-FLOAT-POSITIVE-INFINITY
>
> * (expt (/ -1 .0) 0.5)
> #.SB-EXT:SINGLE-FLOAT-POSITIVE-INFINITY
>   

Thanks, Ted.   This is instructive.  R also handles square roots in a 
sort of 3)-like way, but inconsistently and without maintaining z = 
sqrt(z) * sqrt(z).  Arg works as expected, consistent with atan2 (what 
we use now); but Arg(sqrt(z)) = Arg(z)/2 fails in the lower half-plane.

 > z <- complex(real=Inf, imaginary=0)
 > Arg(z)
[1] 0
 > sqrt(z)
[1] Inf+0i
 > Arg(sqrt(z))
[1] 0
 > sqrt(z)*sqrt(z)
[1] Inf+NaNi
 > z <- complex(real=0, imaginary=Inf)
 > Arg(z)
[1] 1.570796
 > sqrt(z)
[1] Inf+Infi
 > Arg(sqrt(z))
[1] 0.7853982
 > sqrt(z)*sqrt(z)
[1] NaN+Infi
 > z <- complex(real=-Inf, imaginary=Inf)
 > sqrt(z)
[1] Inf+Infi
 > Arg(z)
[1] 2.356194
 > Arg(sqrt(z))
[1] 0.7853982
 > sqrt(z)*sqrt(z)
[1] NaN+Infi
 > z <- complex(real=-Inf, imaginary=0)
 > Arg(z)
[1] 3.141593
 > sqrt(z)
[1] 0+Infi
 > Arg(sqrt(z))
[1] 1.570796
 > sqrt(z)*sqrt(z)
[1] -Inf+NaNi
 > z <- complex(real=-Inf, imaginary=-Inf)
 > Arg(z)
[1] -2.356194
 > sqrt(z)
[1] Inf-Infi
 > Arg(sqrt(z))
[1] -0.7853982
 > sqrt(z)*sqrt(z)
[1] NaN-Infi

The problem with 3) is that it unfortunately *only* works for sqrt and 
4th roots.  Doh!  The problem is that pi/4 splits are the only args that 
you can "reasonably" represent with real and imaginary parts from {0, 
+-Inf}.  The implementation of atan2 collapses all args for pairs with 
an infinite element into these, which makes sense from a limit perspective.

So I think we are left with either 1) or 2) or finding a "standard" impl 
or spec to imitate.  Unfortunatly, the C99x standard is not public and I 
don't recall complex nth roots being included there in any case.   If 
anyone knows of a reasonable public standard to implement here, that 
would be great to at least look at.

An improvement to 1), lets call it 1.5),  that would make behavior a 
little less strange would be to avoid the NaNs from Inf * 0 and the 
"spurious" Infs from Inf * cos(phi) or Inf * sin(phi) where rounding 
error has moved phi off an analytical 0.  You can see examples of both 
of these in first two infinite test cases in ComplexTest.java as it is 
today.  If  the code was changed to work around these cases,  the square 
roots of Inf + 0i, for example, would be {Inf + 0i, -Inf + 0i}.
 I have not tested this, but it looks to me like the NaN in the first 
root is from Inf * 0 and the -Inf in the imaginary part of the second 
root is from Inf * sin(phi) with phi  from 2 * pi / 2 just slightly 
greater than pi.

Phil

> *
>
> On Fri, Jan 2, 2009 at 8:19 PM, Ted Dunning <ted.dunning@gmail.com> wrote:
>
>   
>> I think that (2) is the easiest for a user to understand.  Obviously, the
>> documentation aspect of (1) should be brought along.
>>
>> The behavior of Common Lisp is always instructive in these regards.  The
>> complex arithmetic was generally very well thought out.
>>
>> On Fri, Jan 2, 2009 at 7:14 PM, Phil Steitz <phil.steitz@gmail.com> wrote:
>>
>>     
>>>  ... I noticed another thing.   Behavior for complex numbers with NaN and
>>>       
>>>> infinite parts is, well, "interesting."  It is consistent with what we do
>>>> elsewhere in this class to just apply computational formulas and document
>>>> behavior for infinite and NaN; but the current implementation is hard to
>>>> document correctly and the results are a little strange for infinite values.
>>>>
>>>> I see three reasonable choices here, and am interested in others'
>>>> opinions.  Please select from the following or suggest something else.
>>>>
>>>> 1) Leave as is and just point out that the computational formula will
>>>> behave strangely for infinite values.
>>>>
>>>> 2) return {Complex.NaN} or {Complex.INF} respectively when the number has
>>>> a NaN or infinite part.
>>>>
>>>> 3) return {Complex.NaN} when either part is NaN; but for infinite values,
>>>> compute the argument using getArgument (atan2), generate the arguments for
>>>> the roots from this and select the real/im parts of the roots from {0, inf,
>>>> -inf} to match the argument (this will always be possible because atan2
>>>> always returns a multiple of pi/4 for infinite arguments).  For example,
the
>>>> 4th roots of real positive infinity would be {inf + 0i, 0 + infi, -inf +
0i,
>>>> 0 + -infi}
>>>>
>>>> 2) is probably the most defensible mathematically; but 3) is closer to
>>>> the spirit of C99x.  Unfortunately, since our implementation of multiply
is
>>>> 2)-like, 3) does not guarantee that nth roots actually solve the equation
>>>> r^n = z.
>>>>
>>>>         
>
>
>   


---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org


Mime
View raw message