stla / EllipticFunctions.jl

Jacobi theta functions and related functions.
https://stla.github.io/EllipticFunctions.jl/dev/
MIT License
16 stars 4 forks source link

10x speedup possible for complete ellipticK by not calling ellipticF #34

Closed ahbarnett closed 9 months ago

ahbarnett commented 10 months ago

Dear stla, Thanks for a full-featured library that can handle the complex plane! To verify various branch cut limits, I coded a simple ellipticK, but found it a lot faster, at least on my laptop. I think the speed difference is because you call ellipticF with \phi=pi/2, and that routine has more generality. I just use Gauss' formula from the AGM (which converges extremely fast almost everywhere in the C plane); I think this is the standard method (wikipedia, etc). This is reverse of your repo which gets agm via calling K. Sorry, I haven't had time to make an arb-prec stopping test, safety exit, or docs. This is really dumb code (I checked it never exceeds 7 iterations):

function AGM(a0,b0)             # compute arithmetic-geometric mean
    a = Complex(a0); b = Complex(b0);   # needed even if a and b real, since ab<0 poss
    while abs(a-b)/abs(a) > 1e-15
        a,b = (a+b)/2, sqrt(a*b)
    end
    (a+b)/2
end
ellipticKA(m) = pi/(2*AGM(sqrt(1-Complex(m)),1))    # Gauss 1799 theorem

Timing (I checked rand is insignificant here):

julia> @time for i=1:1000000; ellipticK(2.0*rand()+1im*rand()); end
  2.316693 seconds
julia> @time for i=1:1000000; ellipticKA(2.0*rand()+1im*rand()); end
  0.232196 seconds

Accuracy is fine (at least in that chunk of C-plane):

julia> n=100000; maximum([abs(ellipticK(m)-ellipticKA(m)) for m in 2*rand(n)+1im*rand(n)])
2.085919684541963e-15

It also has more standard branch cuts (see Issue #33). Just a suggestion. I could try to PR it if you want, but you know your own library org and style better. Thanks and best wishes, Alex

stla commented 9 months ago

Thanks. I recently discovered this formula, but I didn't test it. You're right that seems better.

stla commented 9 months ago

I've done the implementation. Thanks.