Mailing-List: contact dev-help@commons.apache.org; run by ezmlm
Precedence: bulk
Reply-To: "Commons Developers List" <dev@commons.apache.org>
Received-SPF: neutral (athena.apache.org: local policy)
Message-ID: <495F9C2E.1040706@free.fr>
Date: Sat, 03 Jan 2009 18:11:10 +0100
From: Luc Maisonobe <Luc.Maisonobe@free.fr>
User-Agent: Thunderbird 2.0.0.18 (X11/20081125)
MIME-Version: 1.0
To: Commons Developers List <dev@commons.apache.org>
Subject: Re: [math] Complex nth roots
References: <495D5168.4060702@gmail.com>
 <495EB9DC.90005@gmail.com>	 <495ED82B.40006@gmail.com>
	 <c7d45fc70901022019v3305cc3dg3b2f9c03703c95ed@mail.gmail.com>
 <c7d45fc70901022045g25f25d74j724e48894c523ac4@mail.gmail.com>
 <495F87C7.8020002@gmail.com>
In-Reply-To: <495F87C7.8020002@gmail.com>
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

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.

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