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.dunn...@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.ste...@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-unsubscr...@commons.apache.org
For additional commands, e-mail: dev-h...@commons.apache.org
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org
For additional commands, e-mail: dev-h...@commons.apache.org
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org
For additional commands, e-mail: dev-h...@commons.apache.org
