handle branch cut using signed zero for imag(m) for ellipticK

[ahbarnett]

Dear stla,
Thanks for an timely library! Certainly helped me, since Matlab only covers (0,1) on real axis. I'm going to report some findings re limits on the branch cut of ellipticK, which may be in the weeds, but it had me confused in my application.

I noticed that for k^2=m>1 real, ie on the branch cut, the limit imag(m) = 0^- is taken (ie from below). However, for complex of the same value (with either sign of imaginary zero), the upper limit is taken. The latter seems standard to me, compared to say, IEEE-754 for sqrt of complex, which for negative real (and positive zero in imaginary) takes the limit from above. It is also usual to infer the side of the cut from the imaginary signed zero, but that doesn't happen:

julia> ellipticK(4.0)                                          # real input
0.8428751774062979 - 1.0782578237498217im
julia> ellipticK(4.0 +1e-15im)
0.8428751774062979 + 1.0782578237498217im
julia> ellipticK(4.0 -1e-15im)                         # ... matches this side of branch
0.8428751774062979 - 1.0782578237498217im
julia> m=Complex(4.0,0); ellipticK(m)          # complex input with +0 imag,  correct up side
0.8428751774062979 + 1.0782578237498217im
julia> m = Complex(-2.0)^2; ellipticK(m)      # hack to make m = 4-0im;  sign of imag ignored
0.8428751774062979 + 1.0782578237498217im

I'll point out that a simple AGM implementation of K with IEEE sqrt gets the sides of the cut I would consider standard, including heeding the sign of the imaginary zero.

Anyway this is a detail, apart from the difference between real m and complex m for the same number would be the main thing to fix. Maybe the above should be documented if it isn't going to be changed? Cheers, Alex

[stla]

Thanks! I'm not an expert in elliptic integrals and I don't understand what is going on here. I will implement ellipticK from the AGM.