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.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