Elliptic functions of complex argument

[Gregstrq]

As far as I understand, the elliptic functions in this package are defined only for real values of the first argument u. At the same time, it is good to be able to calculate them across the whole complex plane (all in all, the defining property of elliptic functions as doubly periodic meromorphic functions is intimately linked with the complex plane).

As a consequence, I have the following question: are there any plans to add the support of complex arguments? If not, could you advise some packages where such functionality is already implemented?

[ali-ramadhan]

I also need complex sn and cn but only Mathematica seems to support complex arguments.

Found these definitions in Abramowitz & Stegun which might be useful for the next person who stumbles upon this issue:

[ali-ramadhan]

Here's what I came up with:

import Elliptic.Jacobi: sn, cn, dn

"""
    sn(z::Complex, m::Real)

Compute the Jacobi elliptic function sn(z | m) following Abramowitz & Stegun (1964), Eq. 16.21.2.
"""
@inline function sn(z::Complex, m::Real)
    s  = sn(real(z), m)
    s₁ = sn(imag(z), 1-m)
    c  = cn(real(z), m)
    c₁ = cn(imag(z), 1-m)
    d  = dn(real(z), m)
    d₁ = dn(imag(z), 1-m)
    return (s * d₁ + im * c * d * s₁ * c₁) / (c₁^2 + m * s^2 * s₁^2)
end

"""
    cn(z::Complex, m::Real)

Compute the Jacobi elliptic function cn(z | m) following Abramowitz & Stegun (1964), Eq. 16.21.3.
"""
@inline function cn(z::Complex, m::Real)
    s  = sn(real(z), m)
    s₁ = sn(imag(z), 1-m)
    c  = cn(real(z), m)
    c₁ = cn(imag(z), 1-m)
    d  = dn(real(z), m)
    d₁ = dn(imag(z), 1-m)
    return (c * c₁ - im * s * d * s₁ * d₁) / (c₁^2 + m * s^2 * s₁^2)
end

"""
    dn(z::Complex, m::Real)

Compute the Jacobi elliptic function dn(z | m) following Abramowitz & Stegun (1964), Eq. 16.21.4.
"""
@inline function dn(z::Complex, m::Real)
    s  = sn(real(z), m)
    s₁ = sn(imag(z), 1-m)
    c  = cn(real(z), m)
    c₁ = cn(imag(z), 1-m)
    d  = dn(real(z), m)
    d₁ = dn(imag(z), 1-m)
    return (d * c₁ * d₁ - im * m * s * c * s₁) / (c₁^2 + m * s^2 * s₁^2)
end

Seems to work well for my conformal maps, but not sure which tests would be good to add before opening a PR.

[Gregstrq]

I have tried to search for the methods of computation of Jacobi Elliptic functions of complex argument, and I did not see any other methods besides the one you suggested.

Also. It seems like a good idea to add a function that computes all three Jacobi Elliptic functions -- sn, cn, and dn -- simultaneously.

[beddalumia]

I found a paper providing a nice set of algorithms to deal with incomplete and complete elliptic integrals with complex input (both the amplitude \phi and the parameter m=k^2).

-> Go to TL;DR

I did not implement nor check all the content, but it worked very very well for the complete elliptic integral of the first kind, K(m), which I needed to compute for generic m ∈ ℂ, in a many-body-physics code I'm developing.

Specifically I referred to the definition given in eq. 42 (in your notation: K(m) = RF(1-m,1,0)) and the ad hoc fast algorithm given in eqs. 43-45, which happens to be the slightest modification of the current implementation in base MATLAB (ellipke function, which gives K(m) and E(m) for 0 ≤ m ≤ 1).

So I just applied these diffs to the original MATLAB code for ellipke:

- function [k,e] = ellipke(m,tol)
+ function [k,e] = complex_ellipke(m,tol)
  %ELLIPKE Complete elliptic integral.
  %   [K,E] = ELLIPKE(M) returns the value of the complete elliptic
  %   integrals of the first and second kinds, evaluated for each
- %   element of M.  As currently implemented, M is limited to 0 <= M <= 1.
+ %   element of M. Domain: M ∈ ℂ (see [2] for details).

  [...]

  % References:
  %   [1] M. Abramowitz and I.A. Stegun, "Handbook of Mathematical
  %       Functions" Dover Publications", 1965, 17.6.
+ %   [2] B.C. Carlson, "Numerical computation of real or complex 
+ %       elliptic integrals" Numerical Algorithms volume 10, 
+ %       pages 13–26 (1995) [arXiv:math/9409227]

  [...]

- if ~isreal(m) || ~isreal(tol)
+ if ~isreal(tol)
     error(message('MATLAB:ellipke:ComplexInputs'))
  end
- if any(m(:) < 0) || any(m(:) > 1)
-    error(message('MATLAB:ellipke:MOutOfRange'));
- end

  [...]

- w1 = 2^i1*c1.^2;
+ w1 = 2^i1*norm(c1).^2;

  [...]

  end

As you can see, except for the routinary checks on the input, there is really just on line changed: a norm instead of a square. To my big surprise, this passed all the tests I could think of against ellipticK, the symbolic implementation provided by MATLAB's Symbolic Math Toolbox, with a speed gain of ≈300x with respect to it (basically the speed is that of ellipke, of course). As amazed as only child could be, I decided to include it in my code, dropping ellipticK, and so far everything seems to work like a charm, with blazing speed. (I'm considering writing to MathWorks, since to me at this point it seems they are really kind of "nerfing" ellipke with no reason, so that I will point them to the paper by Carlson, provide some evidence and see if they would like to modify the function in a future release...).

---

TL;DR

So, all this was to say:

1. I found an amazing paper, dealing with the topic of the issue

2. I had a wonderful experience with one of the algorithms it provides, namely I had to modify in the slightest detail an existing well-known algorithm

3. Maybe you could do the same with all the content, i.e. implement the provided generic algortihms for RF, RC, etc. and thus upgrade the code to the complex domain with a very moderate effort. 🙃

FInal Note

I've done some research: apparently also scipy's current implementations cannot deal with complex parameters. This is quite suprising to me, but well... we could consider to open a issue therein as well.