# §22.20 Methods of Computation

## §22.20(i) Via Theta Functions

A powerful way of computing the twelve Jacobian elliptic functions for real or complex values of both the argument $z$ and the modulus $k$ is to use the definitions in terms of theta functions given in §22.2, obtaining the theta functions via methods described in §20.14.

## §22.20(ii) Arithmetic-Geometric Mean

Given real or complex numbers $a_{0},b_{0}$, with $b_{0}/a_{0}$ not real and negative, define

 22.20.1 $\displaystyle a_{n}$ $\displaystyle=\tfrac{1}{2}\left(a_{n-1}+b_{n-1}\right),$ $\displaystyle b_{n}$ $\displaystyle=\left(a_{n-1}b_{n-1}\right)^{1/2},$ $\displaystyle c_{n}$ $\displaystyle=\tfrac{1}{2}\left(a_{n-1}-b_{n-1}\right),$ ⓘ Symbols: $a_{n}$: numbers, $b_{n}$: numbers, $c_{n}$: numbers and $n$: positive Referenced by: §22.20(ii), §22.20(iv), §22.20(iv), §22.20(ii) Permalink: http://dlmf.nist.gov/22.20.E1 Encodings: TeX, TeX, TeX, pMML, pMML, pMML, png, png, png See also: Annotations for §22.20(ii), §22.20 and Ch.22

for $n\geq 1$, where the square root is chosen so that $\operatorname{ph}b_{n}=\tfrac{1}{2}(\operatorname{ph}a_{n-1}+\operatorname{ph}% b_{n-1})$, where $\operatorname{ph}a_{n-1}$ and $\operatorname{ph}b_{n-1}$ are chosen so that their difference is numerically less than $\pi$. Then as $n\to\infty$ sequences $\{a_{n}\}$, $\{b_{n}\}$ converge to a common limit $M=M(a_{0},b_{0})$, the arithmetic-geometric mean of $a_{0},b_{0}$. And since

 22.20.2 $\max\left(\left|a_{n}-M\right|,\left|b_{n}-M\right|,\left|c_{n}\right|\right)% \leq\text{(const.)}\times 2^{-2^{n}},$ ⓘ Symbols: $a_{n}$: numbers, $b_{n}$: numbers, $c_{n}$: numbers, $n$: positive and $M(a_{0},b_{0})$: limit Permalink: http://dlmf.nist.gov/22.20.E2 Encodings: TeX, pMML, png See also: Annotations for §22.20(ii), §22.20 and Ch.22

convergence is very rapid.

For $x$ real and $k\in(0,1)$, use (22.20.1) with $a_{0}=1$, $b_{0}=k^{\prime}\in(0,1)$, $c_{0}=k$, and continue until $c_{N}$ is zero to the required accuracy. Next, compute $\phi_{N},\phi_{N-1},\dots,\phi_{0}$, where

 22.20.3 $\phi_{N}=2^{N}a_{N}x,$ ⓘ Symbols: $x$: real, $a_{n}$: numbers and $\phi_{N}$: approximations Referenced by: §22.20(ii) Permalink: http://dlmf.nist.gov/22.20.E3 Encodings: TeX, pMML, png See also: Annotations for §22.20(ii), §22.20 and Ch.22
 22.20.4 $\phi_{n-1}=\frac{1}{2}\left(\phi_{n}+\operatorname{arcsin}\left(\frac{c_{n}}{a% _{n}}\sin\phi_{n}\right)\right),$

and the inverse sine has its principal value (§4.23(ii)). Then

 22.20.5 $\displaystyle\operatorname{sn}\left(x,k\right)$ $\displaystyle=\sin\phi_{0},$ $\displaystyle\operatorname{cn}\left(x,k\right)$ $\displaystyle=\cos\phi_{0},$ $\displaystyle\operatorname{dn}\left(x,k\right)$ $\displaystyle=\frac{\cos\phi_{0}}{\cos\left(\phi_{1}-\phi_{0}\right)},$ ⓘ Symbols: $\operatorname{cn}\left(\NVar{z},\NVar{k}\right)$: Jacobian elliptic function, $\operatorname{dn}\left(\NVar{z},\NVar{k}\right)$: Jacobian elliptic function, $\operatorname{sn}\left(\NVar{z},\NVar{k}\right)$: Jacobian elliptic function, $K\left(\NVar{k}\right)$: Legendre’s complete elliptic integral of the first kind, $\cos\NVar{z}$: cosine function, $\sin\NVar{z}$: sine function, $x$: real, $k$: modulus and $\phi_{N}$: approximations Referenced by: (22.20.5), §22.20(ii), 3rd item Permalink: http://dlmf.nist.gov/22.20.E5 Encodings: TeX, TeX, TeX, pMML, pMML, pMML, png, png, png Note (effective with 1.0.10): A note was added after (22.20.5) to deal with cases when computation of $\operatorname{dn}\left(x,k\right)$ with (22.20.5) becomes numerically unstable near $x=K$. Suggested 2014-10-20 by Hartmut Henkel See also: Annotations for §22.20(ii), §22.20 and Ch.22

and the subsidiary functions can be found using (22.2.10). This formula for $\operatorname{dn}$ becomes unstable near $x=K$. If only the value of $\operatorname{dn}\left(x,k\right)$ at $x=K$ is required then the exact value is in the table 22.5.1. If both $k$ and $x$ are real then $\operatorname{dn}$ is strictly positive and $\operatorname{dn}\left(x,k\right)=\sqrt{1-k^{2}{\operatorname{sn}^{2}}\left(x,% k\right)}$ which follows from (22.6.1). If either $k$ or $x$ is complex then (22.2.6) gives the definition of $\operatorname{dn}\left(x,k\right)$ as a quotient of theta functions.

### Example

To compute $\operatorname{sn}$, $\operatorname{cn}$, $\operatorname{dn}$ to 10D when $x=0.8$, $k=0.65$.

Four iterations of (22.20.1) lead to $c_{4}=\Sci{6.5}{-12}$. From (22.20.3) and (22.20.4) we obtain $\phi_{1}=1.40213\;91827$ and $\phi_{0}=0.76850\;92170$. Then from (22.20.5), $\operatorname{sn}\left(0.8,0.65\right)=0.69506\;42165$, $\operatorname{cn}\left(0.8,0.65\right)=0.71894\;76580$, $\operatorname{dn}\left(0.8,0.65\right)=0.89212\;34349$.

## §22.20(iii) Landen Transformations

By application of the transformations given in §§22.7(i) and 22.7(ii), $k$ or $k^{\prime}$ can always be made sufficently small to enable the approximations given in §22.10(ii) to be applied. The rate of convergence is similar to that for the arithmetic-geometric mean.

### Example

To compute $\operatorname{dn}\left(x,k\right)$ to 6D for $x=0.2$, $k^{2}=0.19$, $k^{\prime}=0.9$.

From (22.7.1), $k_{1}=\tfrac{1}{19}$ and $x/(1+k_{1})=0.19$. From the first two terms in (22.10.6) we find $\operatorname{dn}\left(0.19,\tfrac{1}{19}\right)=0.999951$. Then by using (22.7.4) we have $\operatorname{dn}\left(0.2,\sqrt{0.19}\right)=0.996253$.

If needed, the corresponding values of $\operatorname{sn}$ and $\operatorname{cn}$ can be found subsequently by applying (22.10.4) and (22.7.2), followed by (22.10.5) and (22.7.3).

## §22.20(iv) Lattice Calculations

If either $\tau$ or $q=e^{i\pi\tau}$ is given, then we use $k={\theta_{2}^{2}}\left(0,q\right)/{\theta_{3}^{2}}\left(0,q\right)$, $k^{\prime}={\theta_{4}^{2}}\left(0,q\right)/{\theta_{3}^{2}}\left(0,q\right)$, $K=\frac{1}{2}\pi{\theta_{3}^{2}}\left(0,q\right)$, and $K^{\prime}=-i\tau K$, obtaining the values of the theta functions as in §20.14.

If $k,k^{\prime}$ are given with $k^{2}+{k^{\prime}}^{2}=1$ and $\Im k^{\prime}/\Im k<0$, then $K,K^{\prime}$ can be found from

 22.20.6 $\displaystyle K$ $\displaystyle=\frac{\pi}{2M(1,k^{\prime})},$ $\displaystyle K^{\prime}$ $\displaystyle=\frac{\pi}{2M(1,k)},$

using the arithmetic-geometric mean.

### Example 1

If $k=k^{\prime}=1/\sqrt{2}$, then three iterations of (22.20.1) give $M=0.84721\;30848$, and from (22.20.6) $K=\pi/(2M)=1.85407\;46773$ — in agreement with the value of $\left(\Gamma\left(\tfrac{1}{4}\right)\right)^{2}/\left(4\sqrt{\pi}\right)$; compare (23.17.3) and (23.22.2).

### Example 2

If $k^{\prime}=1-i$, then four iterations of (22.20.1) give $K=1.23969\;74481+i0.56499\;30988$.

## §22.20(v) Inverse Functions

See Wachspress (2000).

## §22.20(vi) Related Functions

$\operatorname{am}\left(x,k\right)$ can be computed from its definition (22.16.1) or from its Fourier series (22.16.9). Alternatively, Sala (1989) shows how to apply the arithmetic-geometric mean to compute $\operatorname{am}\left(x,k\right)$.

Jacobi’s epsilon function can be computed from its representation (22.16.30) in terms of theta functions and complete elliptic integrals; compare §20.14. Jacobi’s zeta function can then be found by use of (22.16.32).

## §22.20(vii) Further References

For additional information on methods of computation for the Jacobi and related functions, see the introductory sections in the following books: Lawden (1989), Curtis (1964b), Milne-Thomson (1950), and Spenceley and Spenceley (1947).