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 Tue, 06 Jan 2009 03:53:36 GMT
Phil Steitz wrote:
> 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
Committed in r731822.  I noticed one more thing that I just commented 
for now in the code.  It might be faster and/or more accurate to use a 
solver to compute the nth root of the modulus.  Should we consider doing 
this and exposing the associated real function in MathUtils?  If yes, 
what solver with what parameters?

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