Phil Steitz wrote:
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?Luc Maisonobe wrote: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 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 andbehaves 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-INFINITYThanks, 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] 0sqrt(z)[1] Inf+0iArg(sqrt(z))[1] 0sqrt(z)*sqrt(z)[1] Inf+NaNiz <- complex(real=0, imaginary=Inf) Arg(z)[1] 1.570796sqrt(z)[1] Inf+InfiArg(sqrt(z))[1] 0.7853982sqrt(z)*sqrt(z)[1] NaN+Infiz <- complex(real=-Inf, imaginary=Inf) sqrt(z)[1] Inf+InfiArg(z)[1] 2.356194Arg(sqrt(z))[1] 0.7853982sqrt(z)*sqrt(z)[1] NaN+Infiz <- complex(real=-Inf, imaginary=0) Arg(z)[1] 3.141593sqrt(z)[1] 0+InfiArg(sqrt(z))[1] 1.570796sqrt(z)*sqrt(z)[1] -Inf+NaNiz <- complex(real=-Inf, imaginary=-Inf) Arg(z)[1] -2.356194sqrt(z)[1] Inf-InfiArg(sqrt(z))[1] -0.7853982sqrt(z)*sqrt(z)[1] NaN-Infi The problem with 3) is that it unfortunately *only* works for sqrt and4th roots. Doh! The problem is that pi/4 splits are the only args thatyou can "reasonably" represent with real and imaginary parts from {0, +-Inf}. The implementation of atan2 collapses all args for pairs withan 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 Idon'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 implybreaking 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 istoday. If the code was changed to work around these cases, the squareroots 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.Phil
Phil
LucPhil* 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. Thecomplex 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 andinfinite 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 willbehave 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 tothe 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
