How exactly can we write an incomplete elliptic integral of the first kind as a sum of real and imaginary parts?

Question

All I want to do is numerically map the upper-half plane $\mathbb H:=\{z\in\mathbb C:\Im(z)>0\}$ to the unit square $[0,1)^2$. How this can be done is described in the Wikipedia article of the Schwarz-Christoffel mapping.$^1$ The mapping is given by $$f(z)=\int^z\frac1{\sqrt{\zeta(1-\zeta^2)}}\:{\rm d}\zeta.$$

Now the C++ standard only contains the function ellint_1 (see https://en.cppreference.com/w/cpp/experimental/special_functions), which computes an incomplete elliptic integral of the first kind with a real argument. So, in order to implement my mapping, it would be helpful to rewrite it as a sum of real and imaginary parts.

In the book of Abramowitz and Stegun, there is a formula which allows us to rewirte the incomplete elliptic integral of the first kind as a sum of real and imaginary parts:

However, I don't even understand the implicit equations (due to some print error maybe?). So, if anyone could clarify how we exactly need to apply this, I would be really thankful.

In the Handbook of Elliptic Integrals for Engineers and Scientists , there is a different formula:

However, since I have no idea what $k'$ is or how we can determine it, I again don't know how to implement this.

---

$^1$ But without a reference and I wasn't able to find a proof in the literature ... In the literature, I only find other formulas which are claimed to map $\mathbb H$ to a rectangle. Not only that, they somehow seem to map points on the real line (the border of $\mathbb H$) to the vertices of the rectange ... That's confusing to me, since shouldn't we map $-\infty$ to $(0,0)$, $\infty$ to $(0,1)$ and $\infty+{\rm i}\infty$ to $(1,1)$?

Answer 1

The answer is given in some posts linked to yours by MSE.

Another way is to use the Carlson $RF$ integral, which is easy to implement. Here is my Asymptote implementation. The Asymptote language is close to C++. The pair type is the type for complex numbers.

pair CarlsonRF(pair x, pair y, pair z) {
  int xzero = x == (0, 0) ? 1 : 0; 
  int yzero = y == (0, 0) ? 1 : 0; 
  int zzero = z == (0, 0) ? 1 : 0; 
  int nzeros = xzero + yzero + zzero;
  if(nzeros >= 2) {
    abort("At most one of `x`, `y`, `z` can be 0.");
  }
  real dx = realMax;
  real dy = realMax;
  real dz = realMax;
  real epsilon = 2.0 * realEpsilon;
  pair A;
  while(dx > epsilon || dy > epsilon || dz > epsilon) {
    pair lambda = sqrt(x)*sqrt(y) + sqrt(y)*sqrt(z) + sqrt(z)*sqrt(x);
    x = (x + lambda) / 4.0;
    y = (y + lambda) / 4.0;
    z = (z + lambda) / 4.0;
    A = (x + y + z) / 3.0;
    dx = length(1.0 - x/A);
    dy = length(1.0 - y/A);
    dz = length(1.0 - z/A);
  } 
  real E2 = dx*dy + dy*dz + dz*dx;
  real E3 = dy*dx*dz;
  return (1 - E2/10 + E3/14 + E2*E2/24 - 3*E2*E3/44 - 5*E2*E2*E2/208 +
    3*E3*E3/104 + E2*E2*E3/16) / sqrt(A);
}

private pair atanhz(pair z) {
  return 0.5 * log((1+z) / (1-z));
} 

private bool isinf(pair z) {
  return abs(xpart(z)) == inf || abs(ypart(z)) == inf;
} 

private real rfloor(real x) {
    return x - (x % 1);
}

private real rceil(real x) {
    return rfloor(x) + 1;
}

pair ellipticF(pair phi, pair m) {
  if(phi == 0 || m == inf || m == -inf){
    return (0, 0);
  }
  if(xpart(phi) == 0 && ypart(phi) == inf && ypart(m) == 0 && xpart(m) > 0 &&
      xpart(m) < 1){
    return sgn(ypart(phi)) *
      (ellipticF(pi/2, m) - ellipticF(pi/2, 1/m) / sqrt(m));
  }
  if(abs(xpart(phi)) == pi/2 && m == 1){
    return (nan, nan);
  }
  if(xpart(phi) >= -pi/2 && xpart(phi) <= pi/2){
    if(m == 1 && abs(xpart(phi)) < pi/2){
      return atanhz(sin(phi)); 
    }
    if(m == 0){
      return phi;
    }
    pair sine = sin(phi);
    if(isinf(sine)) {
      abort("`sin(phi)` is not finite.");
    }
    pair sine2 = sine*sine;
    pair cosine2 = 1 - sine2;
    pair oneminusmsine2 = 1 - m * sine2;
    return sine * CarlsonRF(cosine2, oneminusmsine2, (1.0, 0.0));
  }
  real k;
  if(xpart(phi) > pi/2) {
    k = rceil((xpart(phi)-pi/2) / pi);
    phi = phi - k * pi;
  } else {
    k = -rfloor((pi/2-xpart(phi)) / pi);
    phi = phi - k * pi;
  }
  return 2*k*ellipticF(pi/2, m) + ellipticF(phi, m);
}
```