From dev-return-112022-apmail-commons-dev-archive=commons.apache.org@commons.apache.org Tue Jan 06 03:54:12 2009
Return-Path:
Delivered-To: apmail-commons-dev-archive@www.apache.org
Received: (qmail 70673 invoked from network); 6 Jan 2009 03:54:12 -0000
Received: from hermes.apache.org (HELO mail.apache.org) (140.211.11.2)
by minotaur.apache.org with SMTP; 6 Jan 2009 03:54:12 -0000
Received: (qmail 70923 invoked by uid 500); 6 Jan 2009 03:54:11 -0000
Delivered-To: apmail-commons-dev-archive@commons.apache.org
Received: (qmail 70827 invoked by uid 500); 6 Jan 2009 03:54:10 -0000
Mailing-List: contact dev-help@commons.apache.org; run by ezmlm
Precedence: bulk
List-Help:
List-Unsubscribe:
List-Post:
List-Id:
Reply-To: "Commons Developers List"
Delivered-To: mailing list dev@commons.apache.org
Received: (qmail 70816 invoked by uid 99); 6 Jan 2009 03:54:10 -0000
Received: from nike.apache.org (HELO nike.apache.org) (192.87.106.230)
by apache.org (qpsmtpd/0.29) with ESMTP; Mon, 05 Jan 2009 19:54:10 -0800
X-ASF-Spam-Status: No, hits=-0.0 required=10.0
tests=SPF_PASS
X-Spam-Check-By: apache.org
Received-SPF: pass (nike.apache.org: domain of phil.steitz@gmail.com designates 74.125.92.145 as permitted sender)
Received: from [74.125.92.145] (HELO qw-out-1920.google.com) (74.125.92.145)
by apache.org (qpsmtpd/0.29) with ESMTP; Tue, 06 Jan 2009 03:54:01 +0000
Received: by qw-out-1920.google.com with SMTP id 4so3448757qwk.60
for ; Mon, 05 Jan 2009 19:53:39 -0800 (PST)
DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed;
d=gmail.com; s=gamma;
h=domainkey-signature:received:received:message-id:date:from
:user-agent:mime-version:to:subject:references:in-reply-to
:content-type:content-transfer-encoding;
bh=dfAr5KGJqykllKs0iWwBvRAUgpmTszVV2D+IjgcbgRk=;
b=XoDHXQqUMjrfhRDuYtZf5sdqWzzlGNnyZaS95ghXPzsk6l4K9h5i4gi15xbXxMoJjY
8LWTYImJj5g0ss8l0KyvWHltkz9mQX9oE+n/O9tZ9ZL1DyCUHhHz7scYDF3kqMPanQG+
uWxtwj2FPZ+htvkt+h3r3L9vWPMvDCyqxnee4=
DomainKey-Signature: a=rsa-sha1; c=nofws;
d=gmail.com; s=gamma;
h=message-id:date:from:user-agent:mime-version:to:subject:references
:in-reply-to:content-type:content-transfer-encoding;
b=F3vTD/ngui5o0d84PPQznRyvVWpc3hDgnC+wkoK6C6OgsyHX9mTm3IXSOJwKCrHoce
rU/zmSD1ohlPLvhVi8W4Ja/JXCjA/9ZP1vazQHFEdS6oKc0Jm2Al/LOih8cjeCVLIFcp
Jb4CraxT4U+jNZl+3r9CTgtFRsiTkvYXgL8Fg=
Received: by 10.215.12.13 with SMTP id p13mr3583428qai.82.1231214019802;
Mon, 05 Jan 2009 19:53:39 -0800 (PST)
Received: from phil-steitzs-macbook-pro.local (pool-98-114-72-207.phlapa.east.verizon.net [98.114.72.207])
by mx.google.com with ESMTPS id 4sm1245879yxj.18.2009.01.05.19.53.38
(version=SSLv3 cipher=RC4-MD5);
Mon, 05 Jan 2009 19:53:38 -0800 (PST)
Message-ID: <4962D5C0.50705@gmail.com>
Date: Mon, 05 Jan 2009 22:53:36 -0500
From: Phil Steitz
User-Agent: Thunderbird 2.0.0.19 (Macintosh/20081209)
MIME-Version: 1.0
To: Commons Developers List
Subject: Re: [math] Complex nth roots
References: <495D5168.4060702@gmail.com> <495EB9DC.90005@gmail.com> <495ED82B.40006@gmail.com> <495F87C7.8020002@gmail.com> <495F9C2E.1040706@free.fr> <495FA17D.4080505@gmail.com>
In-Reply-To: <495FA17D.4080505@gmail.com>
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Virus-Checked: Checked by ClamAV on apache.org
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
>>>> 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
>>>>> 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