elliptic integrals with complex argument

[my-little-repository]

It seems that F is only defined on the real line. It would be nice if it was also defined on the complex domain, i.e. if it could take a complex argument phi. Actually it would be nice to have complex version of all the elliptic integrals and functions.

Surprisingly, few libraries/CAS provide such a feature. I rely on maxima at the moment for this kind of computation but it is slow.

[my-little-repository]

I had a look at the source code in slatec.jl and it seems that it is close to be generic. The only two problems are the hardcoded constants and the precision handling (eps). As a proof of concept I reimplemented inverse-jacobi-sn (this is the function I am interested in) from the maxima code, which is close to the code for DRF in slatec.jl. The result matches the maxima output.

  julia> inverse_jacobi_sn(big(1)+im, 0.5)
 5.02733945756855650035355109506569128759164736123056439213283279370016134997893e-01 + 1.066678964951364908926864502836271523211468266219013303457397076495607038282015im

Here is the code. As you see, it is close to yours and works both for complex and big arguments.

  inverse_jacobi_sn(u, m) = u * bf_rf(1-u*u, 1-m*u*u,1)

  bferrtol(x...) = mapreduce(x-> eps(real(float(x))), min, x)

  function bf_rf(x, y, z)
    xn = x; yn = y; zn = z; a = (xn+yn+zn)/3; an=a; power4=1; n=0
    epslon = max(abs(xn-a),abs(yn-a),abs(zn-a))/bferrtol(x,y,z)
    while power4*epslon > abs(a)
      lam = sqrt(xn*yn) + sqrt(yn*zn) + sqrt(zn*xn)
      power4 *= 0.25
      xn = (xn+lam)*0.25
      yn = (yn+lam)*0.25
      zn = (zn+lam)*0.25
      an = (an+lam)*0.25
      n += 1
    end
    xndev = (x-a)*power4/an; yndev = (y-a)*power4/an; zndev = -xndev-yndev
    ee2 = xndev*yndev - 6zndev*zndev; ee3 = xndev*yndev*zndev
    s = 1 -ee2/4 + ee3/14 + ee2*ee2/24 - ee2*ee3*3/44
    s/sqrt(an)
  end

[nolta]

Hmm. I'm afraid i have to close this, as maxima is GPL.

[ettersi]

See also https://github.com/JuliaMath/SpecialFunctions.jl/pull/79.