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 17:33:49 GMT
Luc Maisonobe wrote:
> Phil Steitz a écrit :
>   
>> 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.
>>     
>
> Solution 2 seems a good compromise to me. It is reasonably logical and
> could probably be implemented simply. Solution 1 looks clumsy and we
> could probably not withstand it on the long run, which would imply
> breaking compatibility some day to go back to a more definitive solution.
>
>   
>> 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.
>>     
>
> Is this simpler or more difficult to implement than 2 ? If simpler, then
> this could be a reasonable solution too.
>   
2) is easiest, since all you have to do is check isInfinite().  The 1.5) 
idea is more complicated and unfortunately the "avoid spurious infs from 
args near analytical zeros" bit requires an arbitrary choice of what 
"near" means and also still does not guarantee that the roots solve the 
definitional equation.  So my vote is for 2).  If there are no 
objections, I will implement that with the other changes.

Phil
> Luc
>
>   
>> 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
>>
>>     
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
>   


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


Mime
View raw message