Jacobi am Implementation

Computing the Jacobi amplitude am(z|m) with Sala's series

Reference: K. L. Sala, "Transformations of the Jacobian Amplitude Function and Its Calculation via the Arithmetic-Geometric Mean," SIAM J. Math. Anal. 20(6), 1989, pp. 1514–1528.

bf-jacobi-am computes $\mathrm{am}(z|m)$ from Sala's series of Gudermannian functions (section 5), which converges for all complex $z$ and $m$ and follows Mathematica's branch cuts for am, and for small $m$ from the Fourier series in the nome. Section 6 describes the whole algorithm. This note derives the two formulas Sala's series uses:

1. $\mathrm{gd}(x + iy) = \mathrm{atan2}(\sinh x, \cos y) + i,\tanh^{-1}(\sin y,\mathrm{sech}, x)$
2. $\mathrm{gd}(A + w) - \mathrm{gd}(A - w) = 2\tan^{-1}\left(\frac{\sinh w}{\cosh A}\right)$

and states where each one, and the related identities $\mathrm{gd}(u) = 2\tan^{-1}(\tanh(u/2))$ and $\mathrm{gd}(u) = \tan^{-1}(\sinh u)$, is valid. Section 4 also gives a more accurate imaginary part as a possible improvement, and section 7 gives the error to expect in tests and other possible improvements.

Throughout, $\tan^{-1}$, $\tanh^{-1}$, $\sinh^{-1}$ and $\log$ denote principal values, and $\mathrm{atan2}(s, c)$ is the two-argument arctangent with values in $(-\pi, \pi]$.

1. Definition and branch cuts

Sala defines $\mathrm{gd}(z)$ as the integral of $\mathrm{sech}$ from $0$ to $z$. The integrand has poles at $z = (2t + 1),i\pi/2$, so gd is made single-valued by cutting the plane along the segments of the imaginary axis

$$\text{from } (4t + 1)\frac{i\pi}{2} \text{ to } (4t + 3)\frac{i\pi}{2}, \qquad t \in \mathbb{Z}.$$

Off these cuts gd is analytic, with $\mathrm{gd}'(z) = \mathrm{sech}, z$ and $\mathrm{gd}(0) = 0$. Write $z = x + iy$ with $x, y$ real.

2. $\mathrm{gd}(u) = 2\tan^{-1}(\tanh(u/2))$ for all complex $u$ off the cuts

Let $f(u) = 2\tan^{-1}(\tanh(u/2))$.

Derivative. Where $f$ is analytic,

$$f'(u) = \frac{2 \cdot \frac12 \mathrm{sech}^2(u/2)}{1 + \tanh^2(u/2)} = \frac{1}{\cosh^2(u/2) + \sinh^2(u/2)} = \frac{1}{\cosh u} = \mathrm{sech}\, u,$$

and $f(0) = 0$.

Branch cuts. Principal $\tan^{-1}(\zeta)$ is analytic except for $\zeta$ on the rays $\pm i[1, \infty)$. So $f$ fails to be analytic only where $\tanh(u/2) = is$ with $s$ real and $\lvert s\rvert \ge 1$, or where $\tanh(u/2)$ has a pole.

• $\tanh(u/2)$ is purely imaginary only for $u = iy$, where $\tanh(iy/2) = i\tan(y/2)$.
• $\lvert\tan(y/2)\rvert \ge 1$ exactly when $y \in [\pi/2, 3\pi/2] \bmod 2\pi$. This range also contains the poles of $\tanh(u/2)$, at $y = \pi \bmod 2\pi$.

So $f$ is analytic off the segments from $(4t + 1),i\pi/2$ to $(4t + 3),i\pi/2$, which are the same cuts as gd.

Conclusion. $f$ and gd are analytic on the same region, the plane minus these segments. That region is connected. They have the same derivative and agree at $u = 0$, so

$$\mathrm{gd}(u) = 2\tan^{-1}\left(\tanh\frac{u}{2}\right) \quad \text{for all } u \text{ off the cuts.}$$

Equivalently, $\tan(\mathrm{gd}(u)/2) = \tanh(u/2)$.

On the cuts. There $\tanh(u/2)$ lies on a cut of $\tan^{-1}$, so the value returned depends on how the complex atan in use treats an argument on its own cut.

This identity is what the difference formula in section 5 is built on. It could also be used to compute gd directly; section 4 compares its accuracy with the form that is implemented.

3. $\mathrm{gd}(u) = \tan^{-1}(\sinh u)$ only in the strip $\lvert\mathrm{Im}, u\rvert < \pi/2$

For real $u$ the identity $\mathrm{gd}(u) = \tan^{-1}(\sinh u)$ holds: both sides vanish at $0$ and have derivative $\mathrm{sech}, u$. It does not hold on the whole complex plane.

With $u = x + iy$,

$$\sinh u = \sinh x \cos y + i \cosh x \sin y.$$

Where $\cos y = 0$, this is $\pm i\cosh x$, and $\lvert\cosh x\rvert \ge 1$ for every $x$. So $\tan^{-1}(\sinh u)$ has cuts along the whole horizontal lines $\mathrm{Im}, u = \pi/2 + k\pi$, which gd does not have.

Inside the strip $\lvert\mathrm{Im}, u\rvert < \pi/2$, $\sinh u$ never lands on $\pm i[1, \infty)$: that needs $\cos y = 0$, or $x = 0$ and $\lvert\sin y\rvert \ge 1$. So $\tan^{-1}(\sinh u)$ is analytic there, and by the same argument as in section 2 it equals $\mathrm{gd}(u)$.

Outside the strip the two differ by a multiple of $\pi$. For example, at $u = 1 + 2i$:

$$\begin{aligned} \mathrm{gd}(u) &= 1.911125\ldots + 0.676553\ldots i \\\ \tan^{-1}(\sinh u) &= -1.230468\ldots + 0.676553\ldots i \end{aligned}$$

4. The real and imaginary parts of gd

Sala's eq. (2.11) gives, for $x \ne 0$,

$$\begin{aligned} \mathrm{gd}(x + iy) &= \alpha + i\beta, \\\ \alpha &= \mathrm{gd}(x) + \tan^{-1}(\mathrm{csch}\, x) - \tan^{-1}(\cos y \mathrm{csch}\, x), \\\ \beta &= \tanh^{-1}(\sin y \mathrm{sech}\, x). \end{aligned}$$

The implementation rewrites the real part as

$$\alpha = \mathrm{atan2}(\sinh x, \cos y)$$

and uses Sala's $\beta$ unchanged:

(gd (z)
  ;; Eqn 2.11, with the real part written as atan2(sinh x, cos y)
  (let ((x (realpart z))
        (y (imagpart z)))
    (complex (atan (sinh x) (cos y))
             (atanh (* (sin y) (sech x))))))

A possible improvement is to rewrite the imaginary part as

$$\beta = \sinh^{-1}\left(\frac{\sin y}{\sqrt{\sinh^2 x + \cos^2 y}}\right),$$

which avoids the cancellation in $\tanh^{-1}$ near the branch points.

Real part

Let $s = \sinh x$ and $c = \cos y$, with $x \ne 0$.

1. For real $x$, $\mathrm{gd}(x) = \tan^{-1}(\sinh x) = \tan^{-1} s$ (section 3, real case).

2. For real $s \ne 0$, $\tan^{-1} s + \tan^{-1}(1/s) = \mathrm{sign}(s),\pi/2$, and $\mathrm{sign}(s) = \mathrm{sign}(x)$. Since $\mathrm{csch}, x = 1/s$,

   \alpha = \mathrm{sign}(s)\frac{\pi}{2} - \tan^{-1}\frac{c}{s}.
   

3. For $s \ne 0$, the two-argument arctangent is

   \mathrm{atan2}(s, c) =
   \begin{cases}
   \frac{\pi}{2} - \tan^{-1}\frac{c}{s}, & s > 0 \quad (\text{values in } (0, \pi)), \\
   -\frac{\pi}{2} - \tan^{-1}\frac{c}{s}, & s < 0 \quad (\text{values in } (-\pi, 0)),
   \end{cases}
   

   which is $\mathrm{sign}(s),\pi/2 - \tan^{-1}(c/s)$ in both cases.

Therefore $\alpha = \mathrm{atan2}(\sinh x, \cos y)$.

At $x = 0$. Sala's formula divides by zero at $x = 0$ and he treats that case separately. The atan2 form needs no special case:

• $\cos y > 0$: $\mathrm{atan2}(0, \cos y) = 0$, matching Sala's $x = 0$ case.
• $\cos y < 0$: the point lies on a cut. atan2 gives $\pi$, which is the limit from $x > 0$. With IEEE signed zeros, $x = -0.0$ gives $-\pi$, the limit from $x < 0$. So the side of the cut is chosen explicitly, unlike in the $2\tan^{-1}(\tanh(u/2))$ form.

The atan2 form also avoids the cancellation in Sala's form for small $x$, where $\mathrm{gd}(x) + \tan^{-1}(\mathrm{csch}, x)$ and $\tan^{-1}(\cos y \mathrm{csch}, x)$ are both near $\pi/2$.

Imaginary part (possible improvement)

Sala's $\beta = \tanh^{-1}(\sin y \mathrm{sech}, x)$, as implemented, loses accuracy near the branch points $x = 0$, $\sin y = \pm 1$. There $\sin y \mathrm{sech}, x$ is close to $1$, so $\tanh^{-1}$ cancels, and in floating point the product can round to $1$ or above, where $\tanh^{-1}$ fails.

From $\tanh\beta = \sin y / \cosh x$,

$$\sinh\beta = \frac{\tanh\beta}{\sqrt{1 - \tanh^2\beta}} = \frac{\sin y}{\sqrt{\cosh^2 x - \sin^2 y}} = \frac{\sin y}{\sqrt{\sinh^2 x + \cos^2 y}},$$

using $\cosh^2 x - \sin^2 y = \sinh^2 x + 1 - \sin^2 y = \sinh^2 x + \cos^2 y$. So

$$\beta = \sinh^{-1}\left(\frac{\sin y}{\sqrt{\sinh^2 x + \cos^2 y}}\right).$$

The sum under the square root has no cancellation. It is zero only at $x = 0$, $\cos y = 0$, the branch points, where gd is singular.

Overflow

In double precision $\sinh^2 x$ overflows for $\lvert x\rvert \gtrsim 355$. Scaling by the larger of $\lvert\sinh x\rvert$ and $\lvert\cos y\rvert$, as hypot does, avoids that. $\sinh x$ itself, and $\cosh x$ in the implemented $\mathrm{sech}, x$, overflow for $\lvert x\rvert \gtrsim 710$; in that range $\alpha = \mathrm{sign}(x),\pi/2$ and $\beta = 0$ to working precision, so it would need its own branch. With bigfloats the exponent range is effectively unlimited and none of these limits arises. In bf-jacobi-am the reduction of section 6.1 keeps the arguments small enough that these limits are not reached.

The $\sinh^{-1}$ form with scaling:

(gd (z)
  (let* ((x (realpart z))
         (y (imagpart z))
         (sx (sinh x))
         (cy (cos y))
         (big (max (abs sx) (abs cy)))
         (small (min (abs sx) (abs cy))))
    (complex (atan sx cy)
             (asinh (/ (sin y)
                       (* big (sqrt (+ 1 (expt (/ small big) 2)))))))))

Accuracy comparison

Double precision, compared with a 50-digit reference evaluated at the same double inputs, so this measures the algorithm's error and not the conditioning of gd. 20,000 random points with $10^{-12} \le \lvert x\rvert \le 40$ and $\lvert y\rvert \le 10$, plus points near the branch points:

form	worst rel. error, real	worst rel. error, imag
atan2 + Sala's $\tanh^{-1}$ (implemented)	$4 \times 10^{-16}$	$3 \times 10^{-10}$; fails at $x = 10^{-10}$, $y = 1.5707963$
$2\tan^{-1}(\tanh(u/2))$	$1 \times 10^{-12}$	$1 \times 10^{-13}$
atan2 + $\sinh^{-1}$ form	$4 \times 10^{-16}$	$4 \times 10^{-16}$

In both of the first two rows, the worst cases are near the branch points. The $\sinh^{-1}$ form also needs only real functions, where $2\tan^{-1}(\tanh(u/2))$ needs complex $\tanh$ and complex $\tan^{-1}$.

5. The difference gd(A + w) − gd(A − w)

The series

With $K = K(m)$, $K' = K(1 - m)$, and

$$c = \frac{\pi K}{K'}, \qquad w = \frac{\pi z}{2K'} = \frac{c\,z}{2K},$$

Sala's series for the amplitude is

$$\mathrm{am}(z|m) = \mathrm{gd}(w) + \sum_{n \ge 1} \left[\mathrm{gd}(cn + w) - \mathrm{gd}(cn - w)\right].$$

For $n \ge 1$, write $A = cn$. For small $\lvert w\rvert$ the two terms in each bracket nearly cancel, so we want a form that computes the difference directly.

The identity

$$\mathrm{gd}(A + w) - \mathrm{gd}(A - w) = 2\tan^{-1}\left(\frac{\sinh w}{\cosh A}\right).$$

Derivation. Let $D = \mathrm{gd}(A + w) - \mathrm{gd}(A - w)$, and let $a = (A + w)/2$, $b = (A - w)/2$, so that $a - b = w$ and $a + b = A$.

1. By section 2, $\tan(\mathrm{gd}(u)/2) = \tanh(u/2)$. So with $\alpha = \mathrm{gd}(A + w)/2$ and $\beta = \mathrm{gd}(A - w)/2$,

   \tan\alpha = \tanh a, \qquad \tan\beta = \tanh b, \qquad \frac{D}{2} = \alpha - \beta.
   

2. Tangent subtraction, an algebraic identity valid for complex values:

   \tan\frac{D}{2} = \frac{\tanh a - \tanh b}{1 + \tanh a \tanh b}.
   

3. Multiply the numerator and denominator by $\cosh a \cosh b$:

   \begin{aligned}
   \text{numerator:} &\quad \sinh a \cosh b - \cosh a \sinh b = \sinh(a - b) = \sinh w, \\
   \text{denominator:} &\quad \cosh a \cosh b + \sinh a \sinh b = \cosh(a + b) = \cosh A,
   \end{aligned}
   

   so

   \tan\frac{D}{2} = \frac{\sinh w}{\cosh A}.
   

4. Step 3 determines $D/2$ only modulo $\pi$. The branch is fixed in the next subsection.

Check at $A = 0$. The identity gives $\mathrm{gd}(w) - \mathrm{gd}(-w) = 2\mathrm{gd}(w)$ on the left and $2\tan^{-1}(\sinh w)$ on the right. These agree only for $\lvert\mathrm{Im}, w\rvert < \pi/2$ (section 3). This is not a contradiction: the branch condition below requires $\sinh\lvert w\rvert < \cosh 0 = 1$, i.e. $\lvert w\rvert < \sinh^{-1} 1 \approx 0.881 < \pi/2$. The series never uses $A = 0$; its $n = 0$ term is $\mathrm{gd}(w)$ itself.

Where the principal branch is correct

Let $\zeta = \sinh w / \cosh A$. The identity $D = 2\tan^{-1}\zeta$ holds with principal $\tan^{-1}$ when

$$\sinh\lvert w\rvert < \lvert\cosh A\rvert \quad\text{and}\quad \lvert\mathrm{Re}\, w\rvert < \mathrm{Re}\, A.$$

For real $A$ the second condition is implied by the first (shown below), but bf-jacobi-am is also called with complex $m$ and with real $m$ outside $[0, 1]$. Then $K' = K(1 - m)$, and so $c$ and $A = cn$, are complex, and both tests are needed. $\mathrm{Re}, A = n \mathrm{Re}, c &gt; 0$ always (section 6.2).

bf-jacobi-am uses the single test

$$\lvert w\rvert < \mathrm{Re}\, A,$$

which implies both: $\lvert\mathrm{Re}, w\rvert \le \lvert w\rvert &lt; \mathrm{Re}, A$, and $\sinh\lvert w\rvert &lt; \sinh(\mathrm{Re}, A) \le \lvert\cosh A\rvert$, since $\lvert\cosh A\rvert^2 = \sinh^2(\mathrm{Re}, A) + \cos^2(\mathrm{Im}, A)$. It needs no $\sinh$ of $\lvert w\rvert$, which can overflow. The region it gives up, $\mathrm{Re}, A \le \lvert w\rvert &lt; \sinh^{-1}\lvert\cosh A\rvert$, has $\lvert w\rvert$ no smaller than $\mathrm{Re}, A$, so there the gd difference loses at most about a digit.

Real $A$ and real $w$. gd takes values in $(-\pi/2, \pi/2)$, so $D \in (-\pi, \pi)$ and $D/2 \in (-\pi/2, \pi/2)$, the range of principal $\tan^{-1}$. The identity holds for all real $A$ and $w$, with no condition.

Complex $w$. Follow the path $tw$ for $0 \le t \le 1$. At $t = 0$ both sides are $0$. Both sides are analytic in $t$ along the path provided that

• (a) $\zeta(t) = \sinh(tw)/\cosh A$ stays off the cuts $\pm i[1, \infty)$ of $\tan^{-1}$, and
• (b) $A + tw$ and $A - tw$ stay off the cuts of gd.

If both hold, the two sides agree along the whole path, and in particular at $t = 1$.

Condition (a). For complex $v$, $\lvert\sinh v\rvert \le \sinh\lvert v\rvert$, since the power series of $\sinh$ has non-negative coefficients. So along the path

$$\lvert\zeta(t)\rvert \le \frac{\sinh(t\lvert w\rvert)}{\lvert\cosh A\rvert} \le \frac{\sinh\lvert w\rvert}{\lvert\cosh A\rvert} < 1,$$

and $\zeta(t)$ stays inside the unit disk, away from the cuts of $\tan^{-1}$.

Condition (b), for real $A \ge 0$. A point $A \pm tw$ on a cut of gd has real part $0$ and imaginary part at least $\pi/2$ in absolute value. So it needs $\lvert\mathrm{Re}(tw)\rvert = A$ and $\lvert\mathrm{Im}(tw)\rvert \ge \pi/2$, hence

$$\lvert w\rvert^2 \ge \lvert tw\rvert^2 \ge A^2 + \frac{\pi^2}{4}.$$

The condition $\sinh\lvert w\rvert &lt; \cosh A$ gives $\lvert w\rvert &lt; \sinh^{-1}(\cosh A)$. It is enough to show

$$\left[\sinh^{-1}(\cosh A)\right]^2 < A^2 + \frac{\pi^2}{4} \quad \text{for all } A \ge 0.$$

• The derivative of $\sinh^{-1}(\cosh A)$ is $\sinh A / \sqrt{\cosh^2 A + 1} &lt; 1$, so $\sinh^{-1}(\cosh A) \le A + \sinh^{-1} 1$, with $\sinh^{-1} 1 \approx 0.8814$. Then $\left[\sinh^{-1}(\cosh A)\right]^2 - A^2 \le 2 \cdot 0.8814,A + 0.7768$, which is below $\pi^2/4 \approx 2.4674$ for $A &lt; 0.95$.

• For $A \ge 0.95$, use $\sqrt{C^2 + 1} \le C + 1/(2C)$ with $C = \cosh A$:

  \sinh^{-1}(\cosh A) - A \le \ln\left(1 + e^{-2A} + \frac{1}{e^{2A} + 1}\right) \le 2e^{-2A} =: \delta,
  

  so $\left[\sinh^{-1}(\cosh A)\right]^2 - A^2 \le \delta(2A + \delta)$. At $A = 0.95$ this is about $0.66$, and it decreases for larger $A$.

So condition (b) holds whenever condition (a) does, for real $A \ge 0$. The largest value of $\left[\sinh^{-1}(\cosh A)\right]^2 - A^2$ is $\left[\sinh^{-1}(1)\right]^2 \approx 0.777$, at $A = 0$.

Condition (b), for complex $A$. The argument above assumes $A$ real, and for complex $A$ the first test alone is not enough. But if $\lvert\mathrm{Re}, w\rvert &lt; \mathrm{Re}, A$, then along the path

$$\mathrm{Re}(A \pm tw) \ge \mathrm{Re}\, A - t\lvert\mathrm{Re}\, w\rvert > 0,$$

so $A \pm tw$ never reaches the imaginary axis, where all the cuts of gd lie. That gives condition (b) for any $A$ with $\mathrm{Re}, A &gt; 0$.

Example where the second test is needed. Take $m = -6$, a real $m$. Then

$$K = 0.911490\ldots, \qquad K' = K(7) = 0.616803\ldots - 0.911490\ldots i, \qquad c = \frac{\pi K}{K'} = 1.458182\ldots + 2.154853\ldots i,$$

and take $z = 2K' \cdot 1.47/\pi = 0.577223\ldots - 0.853001\ldots i$, so $w = 1.47$. The first test passes, since $\sinh 1.47 = 2.04 &lt; \lvert\cosh c\rvert = 2.10$. But $\mathrm{Re}(c - w) = -0.012$, so $c - tw$ crosses the imaginary axis at $\mathrm{Im} \approx 2.15$, inside the cut of gd from $i\pi/2$ to $3i\pi/2$. The $n = 1$ term then differs by $2\pi$ between the two forms:

$$\begin{aligned} \text{all terms as gd differences:} &\quad \mathrm{am} = 5.77066731555 - 0.730409062547\,i \\\ \text{\$\tan^{-1}\$ form for } n = 1\text{:} &\quad \mathrm{am} = -0.51251799163 - 0.730409062547\,i \end{aligned}$$

Both satisfy $\sin(\mathrm{am}) = \mathrm{sn}(z|m)$ and $\cos(\mathrm{am}) = \mathrm{cn}(z|m)$ to 30 digits, so a check against sn and cn does not detect the error. The single test $\lvert w\rvert &lt; \mathrm{Re}, A$ also rejects this point, since $\lvert w\rvert = 1.47 &gt; \mathrm{Re}, c = 1.458$.

The gd-difference value is the one consistent with Mathematica's branch cuts. The line $\mathrm{Re}(c - w) = 0$ is $z = 2K - isK'$ with $s$ real, and the segment $1 \le s \le 3$ is Mathematica's cut from $2K - iK'$ to $2K - 3iK'$. This $z$ is $2K - iK'(1.3718 + 0.0075i)$, just beside that cut. The gd-difference series jumps by $2\pi$ across the cut, while the $\tan^{-1}$ form continues analytically through it and returns the value from the other side. At $w = 1.44$, on the near side, both forms give $-0.465903 - 0.740155i$.

Why the $\tan^{-1}$ form is needed

In double precision with $A = 3$ and $w = 10^{-8}$:

quantity	value	rel. error
$\mathrm{gd}(3 + 10^{-8})$	1.471304342110472
$\mathrm{gd}(3 - 10^{-8})$	1.4713043401239134
difference	1.9865584732770003e−9	$3.8 \times 10^{-8}$
$2\tan^{-1}(\sinh w / \cosh A)$	1.9865585483886647e−9	$2.6 \times 10^{-16}$
exact	1.986558548388664…e−9

Each gd value carries an absolute error of about $\varepsilon,\mathrm{gd}(A)$, while the true difference is about $2\lvert w\rvert \mathrm{sech}, A$. The subtraction therefore loses about

$$\log_{10}\frac{\mathrm{gd}(A)}{2\lvert w\rvert \mathrm{sech}\, A}$$

digits. Every factor of 10 decrease in $\lvert w\rvert$ costs one more digit. Once $\lvert w\rvert \lesssim \varepsilon,\mathrm{gd}(A)\cosh(A)/2$, the difference is entirely noise. Since $w = \pi z/(2K')$, this is the case of small $\lvert z\rvert$.

Where the $\tan^{-1}$ form fails

Outside the region where the tests hold, principal $\tan^{-1}$ can return the wrong branch. For $0 &lt; m &lt; 1$, $\zeta = \sinh w / \cosh(cn)$ reaches $\pm i[1, \infty)$ on the lines $\mathrm{Im}, w = \pi/2 + k\pi$ with $\lvert\mathrm{Re}, w\rvert \ge cn$, that is on

$$\mathrm{Im}\, z = (2k + 1)K', \qquad \lvert\mathrm{Re}\, z\rvert \ge 2nK.$$

am has no cuts there. At $m = 0.5$, using the $\tan^{-1}$ form for every term gave relative errors of $1.24$ at $z = 5 + 2.5iK'$ and $1.57$ at $z = 9 + 1.5iK'$. In that region $\lvert w\rvert$ is not small, so the gd difference loses little precision and is used instead.

6. The algorithm in bf-jacobi-am

With $K = K(m)$, $K' = K(1 - m)$, $c = \pi K/K'$ and $w = cz/(2K) = \pi z/(2K')$:

1. $m = 1$. $\mathrm{am}(z|1) = \mathrm{gd}(z)$ (DLMF 22.16.5). $K(1)$ is infinite, so the series cannot be used.
2. $1 - m$ rounds to $1$ ($\lvert m\rvert$ below about $1.1 \times 10^{-16}$ for doubles, or about $2^{-\text{fpprec}}$ for bigfloats, with fpprec in bits). $K' = K(1)$ cannot be computed. Use the nome series (6.3) with $q = m/16$, which is the nome to working precision.
3. Otherwise, reduce $z$ by a multiple of $2K$ (6.1), then use the nome series if it converges faster than Sala's series (6.4), and Sala's series otherwise.

6.1 Reduction by 2K

Sala's series is the bilateral sum $\sum_{k \in \mathbb{Z}} \mathrm{gd}(w + kc)$, with the $k$ and $-k$ terms paired. Replacing $z$ by $z + 2K$ replaces $w$ by $w + c$, which reindexes the sum. The truncated sum $\sum_{k=-N}^{N}$ gains $\mathrm{gd}(w + (N+1)c) - \mathrm{gd}(w - Nc)$, which tends to $\mathrm{gd}(+\infty) - \mathrm{gd}(-\infty) = \pi$ because $\mathrm{Re}, c &gt; 0$. So

$$\mathrm{am}(z + 2K | m) = \mathrm{am}(z | m) + \pi$$

for every $m$ other than $0$ and $1$, including complex $m$, and with the same branch cuts as the series. For the nome series the same shift changes $t$ by $2\pi$ and $t/2$ by $\pi$, so the identity holds there too.

bf-jacobi-am computes $\mathrm{am}(z - 2rK) + r\pi$ with

$$r = \mathrm{truncate}\left(\frac{\mathrm{Re}\, w}{\mathrm{Re}\, c}\right),$$

so the reduced $w$ has $\lvert\mathrm{Re}, w\rvert &lt; \mathrm{Re}, c$. Then:

• $\sinh$ and $\cosh$ never overflow in double precision. Without the reduction, $\mathrm{gd}(A + w)$ needs $\sinh$ at up to about $2\lvert\mathrm{Re}, w\rvert$, which overflows for $z \gtrsim 420$ at $m = 0.5$.
• The multiple of $\pi$ is added once, instead of accumulating over $r$ terms each close to $\pi$.
• Fewer terms are needed. Before the series' terms start to decay, $n$ must reach about $\lvert w\rvert/\lvert c\rvert$, which is at most about $1$ after the reduction. For $z = 100$, $m = 0.5$ in double precision this is 13 terms instead of 38.

Truncation, not rounding. Truncating toward zero keeps the reduced $\mathrm{Re}, w$ the same sign as the original, so the reduced am and $r\pi$ have the same sign and their sum does not cancel. Rounding can flip the sign. For example, at $z = 1.5$, $m = 1.5 + i$, $\mathrm{Re}, w / \mathrm{Re}, c = 0.58$: rounding gives $r = 1$ and $\mathrm{am} = \pi + (\approx -2.21)$, and the result was 10 ulps off; truncating gives $r = 0$, 3 ulps off, the same as without reduction.

Effect on accuracy. The reduced argument $z/(2K) - r$ is exact by Sterbenz's lemma, so the reduction adds no rounding error beyond the existing error in $z/(2K)$. That error is about $\lvert z/(2K)\rvert,\varepsilon$ absolute, which makes am wrong by about $\lvert\mathrm{am}\rvert,\varepsilon$, a few ulps. Section 7.1 has examples.

6.2 Why Re c > 0

Let $\tau = iK'/K$, so that the nome is $q = e^{i\pi\tau}$. Then $c = \pi K/K' = i\pi/\tau$ and

$$\mathrm{Re}\, c = \frac{\pi \mathrm{Im}\, \tau}{\lvert\tau\rvert^2}.$$

So $\mathrm{Re}, c &gt; 0$ exactly when $\mathrm{Im}, \tau &gt; 0$, i.e. $\lvert q\rvert &lt; 1$.

• Complex $m$ off the cuts $(-\infty, 0]$ and $[1, \infty)$: $\mathrm{Im}, \tau &gt; 0$ is the standard property of the map from $m$ to $\tau$.
• Real $m &lt; 0$: $K$ is real and positive, and $\mathrm{Im}, \tau = \mathrm{Re}, K' / K$ with $\mathrm{Re}, K(x) &gt; 0$ for real $x &gt; 1$ on either side of its cut.
• Real $m &gt; 1$: $K'$ is real and positive, and $\mathrm{Im}, \tau = K' \mathrm{Re}, K / \lvert K\rvert^2$ with $\mathrm{Re}, K(m) = K(1/m)/\sqrt{m} &gt; 0$.
• $0 &lt; m &lt; 1$: $c$ is real and positive.

So the reduction, which divides by $\mathrm{Re}, c$, and the convergence of Sala's series are safe for all $m$ except $0$ and $1$. $\mathrm{Re}, c \to 0$ as $m \to 0$, because $K'$ grows like $\frac12 \ln(16/m)$; then Sala's series converges slowly, which is what the nome series is for.

6.3 The nome series

For $\lvert\mathrm{Im}, z\rvert &lt; K'$ (DLMF 22.16.9; functions.wolfram.com 09.24.06.0002.01):

$$\mathrm{am}(z|m) = \frac{t}{2} + 2\sum_{n \ge 1} \frac{q^n}{n(1 + q^{2n})} \sin(nt),$$

with $t = \pi z/K$ and $q = \exp(-\pi K'/K)$, computed from the $K$ and $K'$ already needed for the reduction. This $q$ uses the same $K'$ as Sala's series, so the two agree for complex $m$. For small $m$, $K' = K(1 - m)$ is computed from a rounded $1 - m$, but $q \approx m/16$, and an absolute error of about $\varepsilon$ in the effective $m$ changes $q$ by only about $\varepsilon/16$.

Since $\lvert\sin(nt)\rvert \le e^{n\lvert\mathrm{Im}, t\rvert}$, the $n$-th term is bounded by $2\rho^n/n$ with

$$\rho = \lvert q\rvert\, e^{\lvert\mathrm{Im}\, t\rvert} = \lvert q\rvert \exp\left(\pi \left\lvert\mathrm{Im}\, \frac{z}{K}\right\rvert\right).$$

The series converges for $\rho &lt; 1$, which is the strip $\lvert\mathrm{Im}, z\rvert &lt; K'$, where am has no cuts. The loop stops when $2\rho^n/n \le \text{tol}\cdot\lvert\text{sum}\rvert$, not when a term is small: $\sin(nt)$ can be zero for one $n$ while later terms are not. At $z = K/2$, for example, $t = \pi/2$ and $\sin 2t = 0$, so a test on the term stops at $n = 2$ and misses the $n = 3$ term.

At $z = 0$ the sum is $0$ and the bound test never succeeds in bigfloat arithmetic, where $\rho^n$ never underflows, so $t = 0$ returns $0$ directly. This is reachable through the reduction when $z$ is an exact multiple of $2K$.

6.4 Choosing between the series

After the reduction, Sala's terms decay like $e^{-n\mathrm{Re}, c}$, independent of $\mathrm{Im}, z$, and the nome series terms like $\rho^n$. bf-jacobi-am uses the nome series when $\rho &lt; e^{-\mathrm{Re}, c}$. Since $\lvert q\rvert = \exp(-\pi \mathrm{Re}(K'/K))$ and $\mathrm{Re}, c = \pi \mathrm{Re}(K/K')$, taking logs and dividing by $\pi$ gives

$$\left\lvert\mathrm{Im}\, \frac{z}{K}\right\rvert + \mathrm{Re}\, \frac{K}{K'} < \mathrm{Re}\, \frac{K'}{K}.$$

• For real $z$ and real $0 &lt; m &lt; 1$ this is $K &lt; K'$, i.e. $m &lt; 1/2$.
• A nonzero $\mathrm{Im}, z$ moves the switch toward smaller $m$, since it uses up the strip where the nome series converges.
• $\rho &lt; 1$ is implied, so the nome series is never used where it diverges.
• $\mathrm{Im}(z/K)$ is unchanged by the reduction, so the test can use either $z$.

Terms needed at $z = 3 + 0.5i$, 30 digits:

$m$	nome series	Sala's series (approx.)
$10^{-3}$	8	62
$10^{-10}$	3	167
$10^{-20}$	2	($K(1 - m)$ fails in double)
$10^{-300}$	1

7. Accuracy: expected errors and possible improvements

Sections 7.4 and 7.5 and the $\sinh^{-1}$ form of section 4 are not implemented.

7.1 Expected error for large |z|

Without the reduction of 6.1, for real $z$ with $A = cn &lt; w$ each term $\mathrm{gd}(A + w) - \mathrm{gd}(A - w)$ is close to $\pi$, and the sum grows to about $\lvert\mathrm{am}\rvert$. Each addition rounds at half an ulp of the running sum. At $z = 100$, $m = 0.5$ the correctly rounded value is 84.70311272411382 (exact 84.70311272411382440…); the unreduced series gave 84.70311272411385, 2 ulps high.

With the reduction, Maxima under SBCL in double precision, against 60-digit references at the same inputs ($m = 0.5$):

$z$	result	error
100	84.70311272411384	+0.84 ulp
1000	847.1374038213263	−0.56 ulp
3000	2541.6555265646302	−0.47 ulp

The unreduced series overflows for $z \gtrsim 420$.

7.2 Small m: before and after the nome series

Maxima under SBCL, with the previous bf-jacobi-am (Sala's series only) and the current one. References are 50-digit values at the same inputs; times are per call.

Double precision:

call	correct value	before	after
jacobi_am(0.5, 1e-20)	0.5	error: elliptic_kc(1) is undefined	0.5
jacobi_am(10.0, 1e-3)	9.997613139826175…	9.997613139826168 (4 ulps), 84 µs	9.997613139826175 (correctly rounded), 56 µs
jacobi_am(3.0, 1e-10)	2.99999999992150731…	2.999999999921508 (1.6 ulps), 144 µs	2.9999999999215072 (correctly rounded), 66 µs
jacobi_am(0.5, 1e-15)	0.49999999999999998… (rounds to 0.5)	0.5, 132 µs	0.5000000000000001 (1 ulp), 48 µs

Bigfloat with fpprec = 32:

call	correct value to 32 digits	before	after
jacobi_am(3b0, 1b-30)	2.9999999999999999999999999999992	2.9999999999999999999999999999996 (4 in the last digit), 40.9 ms	2.9999999999999999999999999999992, 1.1 ms
jacobi_am(10b0, 1b-3)	9.9976131398261754702124306920947	9.9976131398261754702124306920951 (4 in the last digit), 6.4 ms	9.9976131398261754702124306920948 (1 in the last digit), 1.2 ms
jacobi_am(3b0, 1b-40)	3.0	error: elliptic_kc(1) is undefined	3.0

The gains come from three places:

• Tiny $m$ no longer fails. When $1 - m$ rounds to $1$, the previous code called $K(1)$.
• Fewer terms. For small $m$, Sala's terms decay like $e^{-n\mathrm{Re}, c}$ with $\mathrm{Re}, c \approx 0.38$ at $m = 10^{-10}$, while the nome series decays like $q^n$ with $q \approx m/16$.
• Less accumulated rounding at larger $z$. At $z = 10$ the previous code summed about four terms of size $\pi$; the current code reduces $z$ and adds $4\pi$ once.

It is not correctly rounded everywhere: at jacobi_am(0.5, 1e-15) the rounding in computing $t = \pi z/K$ leaves the result 1 ulp high.

7.3 Expected error in tests of sin(am) against sn

Even with am correct to a few ulps of the computation of $w$, two factors enlarge $\lvert\sin(\mathrm{am}) - \mathrm{sn}\rvert$:

• The series computes $w = \pi z/(2K')$, whose relative error is about $\varepsilon$. Since $d\mathrm{am}/dz = \mathrm{dn}(z)$, this gives an error in am of about $\lvert z\rvert,\lvert\mathrm{dn}(z)\rvert,\varepsilon$. Any method that scales $z$ by a period has this error.
• $\sin(\mathrm{am})$ changes with am at rate $\lvert\cos\mathrm{am}\rvert = \lvert\mathrm{cn}(z)\rvert$.

So

$$\lvert\sin(\mathrm{am}) - \mathrm{sn}\rvert \lesssim \lvert\mathrm{cn}\rvert\,\lvert\mathrm{dn}\rvert\,\lvert z\rvert\,\varepsilon + \lvert\mathrm{am}\rvert\,\lvert\mathrm{cn}\rvert\,\varepsilon,$$

the second term from the rounding of am itself. Both $\lvert\mathrm{cn}\rvert$ and $\lvert\mathrm{dn}\rvert$ become large near the poles of sn.

Example: $z = 2ki$, $m = 0.75$, at about 106 bits (fpprec 32, $\varepsilon \approx 1.2 \times 10^{-32}$). The poles of sn on the imaginary axis are at $i(2s + 1)K' = 1.69i,\ 5.06i,\ 8.43i, \ldots$

$k$	$\lvert\mathrm{sn}\rvert$	am	am error	$\lvert\sin(\mathrm{am}) - \mathrm{sn}\rvert$
1	3.57	$\pi + 1.984i$	$1.3 \times 10^{-31}$	$4.7 \times 10^{-31}$
4	2.55	$1.665i$	$2.3 \times 10^{-31}$	$6.2 \times 10^{-31}$
6	5.71	$-2.444i$	$4.1 \times 10^{-31}$	$2.4 \times 10^{-30}$
9	1.94	$\pi - 1.416i$	$9.6 \times 10^{-31}$	$2.1 \times 10^{-30}$
5	0.11	$\pi + 0.115i$	$2.7 \times 10^{-32}$	$2.7 \times 10^{-32}$

At $k = 6$, $\lvert\mathrm{cn}\rvert,\lvert\mathrm{dn}\rvert,\lvert z\rvert,\varepsilon \approx 4 \times 10^{-30}$. The rtest_elliptic tolerance for this test was raised from 5.5b−32 to 3.2b−30 for this reason. A relative test, closeto((sin(am) - sn)/sn, ...), would need less slack because the test values reach $\lvert\mathrm{sn}\rvert \approx 5.7$.

7.4 Real m > 1: reciprocal-modulus transformation

For $m &gt; 1$, $K(m)$ is complex, so $c$ is complex and a real $z$ gives a result with a spurious imaginary part of order $\varepsilon$. For example, $\mathrm{am}(16 | 2.5) = 24.878832226728505 - 1.95 \times 10^{-15},i$ in double precision, while the true imaginary part is about $3 \times 10^{-41}$.

A&S 16.11 maps $m &gt; 1$ to $\mu = 1/m &lt; 1$. With $v = \sqrt{m},u$,

$$\mathrm{sn}(u|m) = \sqrt{\mu}\,\mathrm{sn}(v|\mu), \qquad \mathrm{cn}(u|m) = \mathrm{dn}(v|\mu),$$

so $e^{i\mathrm{am}(u|m)} = \mathrm{cn}(u|m) + i\mathrm{sn}(u|m) = \mathrm{dn}(v|\mu) + i\sqrt{\mu},\mathrm{sn}(v|\mu)$, and

$$\mathrm{am}(u|m) = -i\log\left(\mathrm{dn}(v|\mu) + i\sqrt{\mu}\,\mathrm{sn}(v|\mu)\right) + 2\pi N, \qquad N = \mathrm{round}\left(\frac{\mathrm{Re}\, v}{4K(\mu)}\right).$$

The $2\pi N$ term. The log gives am only modulo $2\pi$, and its principal value is the branch that is continuous along the real axis. For $m &gt; 1$, Mathematica's cuts for am cross the real axis at $u = 2(2s + 1)K(\mu)/\sqrt{m}$, i.e. at $v = 2(2s + 1)K(\mu)$, and the series jumps by $2\pi$ at each crossing. $N$ counts the crossings between $0$ and $\mathrm{Re}, v$.

Example: $u = 16$, $m = 2.5$. The principal log gives $-0.253909001989843$, which is also the integral of dn along the real axis. $N = 4$, and $-0.253909001989843 + 8\pi = 24.8788322267285$, the series value.

Check. 400 random points with real $m \in [1.05, 20]$ and $u \in [-15, 15]$, half with $\mathrm{Im}, u \in [-4, 4]$, at 25 digits: the principal log alone differed from the series at 354 points, always by a multiple of $2\pi$; with the $2\pi N$ term all 400 agreed.

Advantages. $K(\mu)$ and $K'(\mu)$ are real, so real $u$ gives an exactly real result, and the existing sn and dn code replaces the series.

Implementation notes.

• Real $u$: $\mathrm{am} = \mathrm{atan2}(\sqrt{\mu},\mathrm{sn}(v|\mu), \mathrm{dn}(v|\mu)) + 2\pi N$. Since $\mathrm{dn}(v|\mu) &gt; 0$, the atan2 gives the continuous part directly.
• Complex $u$: compute the real part of the log as an atan2 and the imaginary part, $-\log\lvert\mathrm{dn} + i\sqrt{\mu},\mathrm{sn}\rvert$, with log1p, so that small am keeps its relative accuracy.
• Complex $m$: $K(\mu)$ is complex, and whether $\mathrm{round}(\mathrm{Re}, v / 4K(\mu))$ is still the right crossing count has not been checked.

7.5 Real m < 0: negative-parameter transformation

A&S 16.10 maps $m &lt; 0$ to $\mu = -m/(1 - m) \in (0, 1)$. With $v = u\sqrt{1 - m}$,

$$\mathrm{sn}(u|m) = \frac{\mathrm{sd}(v|\mu)}{\sqrt{1 - m}}, \qquad \mathrm{cn}(u|m) = \mathrm{cd}(v|\mu),$$

so

$$e^{i\mathrm{am}(u|m)} = \frac{\mathrm{cn}(v|\mu) + i\mathrm{sn}(v|\mu)/\sqrt{1 - m}}{\mathrm{dn}(v|\mu)}.$$

For $m &lt; 0$, $K' = K(1 - m)$ is complex, which caused the branch issue in the $m = -6$ example of section 5. This transformation keeps everything real for real $u$. The $2\pi$ correction for complex $u$, the analogue of $N$ in 7.4, has not been worked out or tested.

8. Numerical checks

All at 30–40 digits with mpmath, using Sala's definitions above:

• $\mathrm{gd}(u)$ vs $2\tan^{-1}(\tanh(u/2))$: 20,000 random points with $\lvert\mathrm{Re}, u\rvert \le 4$, $\lvert\mathrm{Im}, u\rvert \le 12$; largest difference $3 \times 10^{-27}$.
• $\mathrm{gd}(u)$ vs $\tan^{-1}(\sinh u)$ on the same points: largest difference $\pi$ (wrong branch outside $\lvert\mathrm{Im}, u\rvert &lt; \pi/2$).
• Sala's $\alpha$ vs $\mathrm{atan2}(\sinh x, \cos y)$: 20,000 random points with $\lvert x\rvert \le 6$, $\lvert y\rvert \le 20$; largest difference $5 \times 10^{-41}$.
• $\mathrm{gd}(A + w) - \mathrm{gd}(A - w)$ vs $2\tan^{-1}(\sinh w / \cosh A)$: 20,000 random points with real $A \in [0.05, 6]$ and complex $w$ satisfying $\sinh\lvert w\rvert &lt; \cosh A$; largest difference $1.6 \times 10^{-30}$.
• The same comparison with complex $A \in [0, 3] \times [-6, 6],i$ and complex $w$ satisfying only $\sinh\lvert w\rvert &lt; \lvert\cosh A\rvert$: 3,862 of 100,000 points differed, all by $\pm 2\pi$. With $\lvert\mathrm{Re}, w\rvert &lt; \mathrm{Re}, A$ added: none.
• $\mathrm{am}(z|m)$ with both tests vs all terms as gd differences, and $\sin(\mathrm{am})$, $\cos(\mathrm{am})$ vs sn, cn: 1,500 random points with $z \in [-3, 3] \times [-3, 3],i$ and $m$ complex in $[-10, 10] \times [-10, 10],i$, real in $[-20, 0]$, or real in $[1, 20]$; largest relative difference $10^{-29}$, largest sn/cn error $7 \times 10^{-28}$.
• The reduction by $2K$ with the single test $\lvert w\rvert &lt; \mathrm{Re}, A$, vs the unreduced series with both tests: 1,500 random points with $\lvert\mathrm{Re}, z\rvert \le 60$, $\lvert\mathrm{Im}, z\rvert \le 4$ and $m$ complex, real negative, real above 1, or in $(0, 1)$; largest relative difference $5 \times 10^{-29}$.
• The nome series vs Sala's series, with $q = \exp(-\pi K'/K)$: 909 points with $m$ real in $\pm[10^{-30}, 0.1]$ or complex with $\lvert m\rvert \lesssim 0.3$, $\lvert\mathrm{Re}, z\rvert \le 30$, $\lvert\mathrm{Im}, z\rvert \le 3$ and $\rho &lt; 0.5$; largest relative difference $7 \times 10^{-28}$.
• Double-precision accuracy of the three forms of gd: see section 4.
• Errors at large $\lvert z\rvert$, at small $m$, in sin(am) vs sn tests, and the reciprocal-modulus check for $m &gt; 1$: see section 7.

Of bf-jacobi-am itself, in Maxima under SBCL 2.2.9:

• Double precision, 300 random points, against 40-digit references. Largest $\lvert\text{error}\rvert / (\varepsilon \max(1, \lvert\mathrm{am}\rvert))$ by kind of $m$: 3.9 for $m \in (0, 1)$; 2.6 for $m &lt; 0$; 8.2 for $m &gt; 1$ (including the spurious imaginary part of section 7.4); 2.9 for complex $m$; 2.1 for $\lvert m\rvert \in [10^{-25}, 10^{-2}]$; 3.1 for $m \in (0.05, 0.95)$ with $\lvert z\rvert$ up to 2000.
• Bigfloat with fpprec = 40, 40 random points with $\lvert\mathrm{Re}, z\rvert \le 300$ and $m$ as above or down to $10^{-60}$: largest relative error $8 \times 10^{-40}$.
• rtest_elliptic, before and after the change: problems 229 and 278 fail both times, since their tolerances were set under CMUCL. Problem 280 failed before and passes after. Problem 249, the $z = 2ki$, $m = 1.75$ bigfloat sin(am) vs sn test, passed before and fails after: its largest difference went from below its 6.3091b−31 tolerance to 6.33b−31.
• run_testsuite() after the change: the only other failure is rtest5 problem 80, describe(sin), because this build has no documentation.