# §19.36 Methods of Computation

## §19.36(i) Duplication Method

Numerical differences between the variables of a symmetric integral can be reduced in magnitude by successive factors of 4 by repeated applications of the duplication theorem, as shown by (19.26.18). When the differences are moderately small, the iteration is stopped, the elementary symmetric functions of certain differences are calculated, and a polynomial consisting of a fixed number of terms of the sum in (19.19.7) is evaluated. For $R_{F}$ the polynomial of degree 7, for example, is

 19.36.1 $1-\tfrac{1}{10}E_{2}+\tfrac{1}{14}E_{3}+\tfrac{1}{24}E_{2}^{2}-\tfrac{3}{44}E_% {2}E_{3}-\tfrac{5}{208}E_{2}^{3}+\tfrac{3}{104}E_{3}^{2}+\tfrac{1}{16}E_{2}^{2% }E_{3},$ ⓘ Symbols: $E_{s}(\mathbf{z})$: symmetric function Referenced by: §19.19, §19.36(i), §19.36(i) Permalink: http://dlmf.nist.gov/19.36.E1 Encodings: TeX, pMML, png See also: Annotations for §19.36(i), §19.36 and Ch.19

where the elementary symmetric functions $E_{s}$ are defined by (19.19.4). If (19.36.1) is used instead of its first five terms, then the factor $(3r)^{-1/6}$ in Carlson (1995, (2.2)) is changed to $(3r)^{-1/8}$.

For both $R_{D}$ and $R_{J}$ the factor $(r/4)^{-1/6}$ in Carlson (1995, (2.18)) is changed to $(r/5)^{-1/8}$ when the following polynomial of degree 7 (the same for both) is used instead of its first seven terms:

 19.36.2 $1-\tfrac{3}{14}E_{2}+\tfrac{1}{6}E_{3}+\tfrac{9}{88}E_{2}^{2}-\tfrac{3}{22}E_{% 4}-\tfrac{9}{52}E_{2}E_{3}+\tfrac{3}{26}E_{5}-\tfrac{1}{16}E_{2}^{3}+\tfrac{3}% {40}E_{3}^{2}+\tfrac{3}{20}E_{2}E_{4}+\tfrac{45}{272}E_{2}^{2}E_{3}-\tfrac{9}{% 68}(E_{3}E_{4}+E_{2}E_{5}).$ ⓘ Symbols: $E_{s}(\mathbf{z})$: symmetric function Referenced by: §19.19 Permalink: http://dlmf.nist.gov/19.36.E2 Encodings: TeX, pMML, png See also: Annotations for §19.36(i), §19.36 and Ch.19

Polynomials of still higher degree can be obtained from (19.19.5) and (19.19.7).

The duplication method starts with computation of $\lambda_{0}=\sqrt{x}\sqrt{y}+\sqrt{y}\sqrt{z}+\sqrt{z}\sqrt{x}$. The reductions in §19.29(i) represent $x,y,z$ as squares, for example $x=U_{12}^{2}$ in (19.29.4). Because $U_{12}$ may be real and negative, or even complex, care is needed to ensure $\sqrt{x}=U_{12}$, and similarly for $y$ and $z$. This precaution is needed only for $\lambda_{0}$. Alternatively, the first duplication is done analytically as in Carlson and FitzSimons (2000), where further information can be found.

### Example

Three applications of (19.26.18) yield

 19.36.3 $R_{F}\left(1,2,4\right)=R_{F}\left(z_{1},z_{2},z_{3}\right),$ ⓘ Symbols: $R_{F}\left(\NVar{x},\NVar{y},\NVar{z}\right)$: symmetric elliptic integral of first kind Permalink: http://dlmf.nist.gov/19.36.E3 Encodings: TeX, pMML, png See also: Annotations for §19.36(i), §19.36(i), §19.36 and Ch.19

where, in the notation of (19.19.7) with $a=-\frac{1}{2}$ and $n=3$,

 19.36.4 \begin{aligned} z_{1}&=2.10985\;99098\;8,\\ z_{3}&=2.15673\;49098\;8,\\ Z_{1}&=0.00977\;77253\;5,\end{aligned}\quad\begin{aligned} z_{2}&=2.12548\;490% 98\;8,\\ A&=2.13069\;32432\;1,\\ Z_{2}&=0.00244\;44313\;4,\end{aligned}\\ {Z_{3}=-Z_{1}-Z_{2}=-0.01222\;21566\;9,}\\ {E_{2}=\Sci{-1.25480\;14}{-4},\quad E_{3}=\Sci{-2.9212}{-7}.} ⓘ Symbols: $E_{s}(\mathbf{z})$: symmetric function, $Z_{j}$: components and $A$ Permalink: http://dlmf.nist.gov/19.36.E4 Encodings: TeX, pMML, png See also: Annotations for §19.36(i), §19.36(i), §19.36 and Ch.19

The first five terms of (19.36.1) suffice for

 19.36.5 $R_{F}\left(1,2,4\right)=0.68508\;58166\dots.$ ⓘ Symbols: $R_{F}\left(\NVar{x},\NVar{y},\NVar{z}\right)$: symmetric elliptic integral of first kind Referenced by: §19.36(ii) Permalink: http://dlmf.nist.gov/19.36.E5 Encodings: TeX, pMML, png See also: Annotations for §19.36(i), §19.36(i), §19.36 and Ch.19

All cases of $R_{F}$, $R_{C}$, $R_{J}$, and $R_{D}$ are computed by essentially the same procedure (after transforming Cauchy principal values by means of (19.20.14) and (19.2.20)). Complex values of the variables are allowed, with some restrictions in the case of $R_{J}$ that are sufficient but not always necessary. The computation is slowest for complete cases. For details see Carlson (1995, 2002) and Carlson and FitzSimons (2000). In the Appendix of the last reference it is shown how to compute $R_{J}$ without computing $R_{C}$ more than once. Because of cancellations in (19.26.21) it is advisable to compute $R_{G}$ from $R_{F}$ and $R_{D}$ by (19.21.10) or else to use §19.36(ii).

Legendre’s integrals can be computed from symmetric integrals by using the relations in §19.25(i). Note the remark following (19.25.11). If (19.25.9) is used when $0\leq k^{2}\leq 1$, cancellations may lead to loss of significant figures when $k^{2}$ is close to 1 and $\phi>\pi/4$, as shown by Reinsch and Raab (2000). The cancellations can be eliminated, however, by using (19.25.10).

Accurate values of $F\left(\phi,k\right)-E\left(\phi,k\right)$ for $k^{2}$ near 0 can be obtained from $R_{D}$ by (19.2.6) and (19.25.13).

Complete cases of Legendre’s integrals and symmetric integrals can be computed with quadratic convergence by the AGM method (including Bartky transformations), using the equations in §19.8(i) and §19.22(ii), respectively.

The incomplete integrals $R_{F}\left(x,y,z\right)$ and $R_{G}\left(x,y,z\right)$ can be computed by successive transformations in which two of the three variables converge quadratically to a common value and the integrals reduce to $R_{C}$, accompanied by two quadratically convergent series in the case of $R_{G}$; compare Carlson (1965, §§5,6). (In Legendre’s notation the modulus $k$ approaches 0 or 1.) Let

 19.36.6 $\displaystyle 2a_{n+1}$ $\displaystyle=a_{n}+\sqrt{a_{n}^{2}-c_{n}^{2}},$ $\displaystyle 2c_{n+1}$ $\displaystyle=a_{n}-\sqrt{a_{n}^{2}-c_{n}^{2}}=c_{n}^{2}/(2a_{n+1}),$ $\displaystyle 2t_{n+1}$ $\displaystyle=t_{n}+\sqrt{t_{n}^{2}+\theta c_{n}^{2}},$ ⓘ Symbols: $n$: nonnegative integer Referenced by: §19.36(ii), §19.36(ii) Permalink: http://dlmf.nist.gov/19.36.E6 Encodings: TeX, TeX, TeX, pMML, pMML, pMML, png, png, png See also: Annotations for §19.36(ii), §19.36 and Ch.19

where $n=0,1,2,\dots$, and

 19.36.7 $\displaystyle 0$ $\displaystyle $\displaystyle t_{0}$ $\displaystyle\geq 0,$ $\displaystyle t_{0}^{2}+\theta a_{0}^{2}$ $\displaystyle\geq 0,$ $\displaystyle\theta$ $\displaystyle=\pm 1.$ ⓘ Permalink: http://dlmf.nist.gov/19.36.E7 Encodings: TeX, TeX, TeX, TeX, pMML, pMML, pMML, pMML, png, png, png, png See also: Annotations for §19.36(ii), §19.36 and Ch.19

Then (19.22.18) implies that

 19.36.8 $R_{F}\left(t_{n}^{2},t_{n}^{2}+\theta c_{n}^{2},t_{n}^{2}+\theta a_{n}^{2}\right)$ ⓘ Symbols: $R_{F}\left(\NVar{x},\NVar{y},\NVar{z}\right)$: symmetric elliptic integral of first kind and $n$: nonnegative integer Referenced by: §19.36(ii) Permalink: http://dlmf.nist.gov/19.36.E8 Encodings: TeX, pMML, png See also: Annotations for §19.36(ii), §19.36 and Ch.19

is independent of $n$. As $n\to\infty$, $c_{n}$, $a_{n}$, and $t_{n}$ converge quadratically to limits $0$, $M$, and $T$, respectively; hence

 19.36.9 $R_{F}\left(t_{0}^{2},t_{0}^{2}+\theta c_{0}^{2},t_{0}^{2}+\theta a_{0}^{2}% \right)=R_{F}\left(T^{2},T^{2},T^{2}+\theta M^{2}\right)=R_{C}\left(T^{2}+% \theta M^{2},T^{2}\right).$

If $t_{0}=a_{0}$ and $\theta=-1$, so that $t_{n}=a_{n}$, then this procedure reduces to the AGM method for the complete integral.

The step from $n$ to $n+1$ is an ascending Landen transformation if $\theta=1$ (leading ultimately to a hyperbolic case of $R_{C}$) or a descending Gauss transformation if $\theta=-1$ (leading to a circular case of $R_{C}$). If $x$, $y$, and $z$ are permuted so that $0\leq x, then the computation of $R_{F}\left(x,y,z\right)$ is fastest if we make $c_{0}^{2}\leq a_{0}^{2}/2$ by choosing $\theta=1$ when $y<(x+z)/2$ or $\theta=-1$ when $y\geq(x+z)/2$.

### Example

We compute $R_{F}\left(1,2,4\right)$ by setting $\theta=1$, $t_{0}=c_{0}=1$, and $a_{0}=\sqrt{3}$. Then

 19.36.10 $\displaystyle c_{3}^{2}$ $\displaystyle=\Sci{6.65}{-12},$ $\displaystyle a_{3}^{2}$ $\displaystyle=2.46209\;30206\;0=M^{2},$ $\displaystyle t_{3}^{2}$ $\displaystyle=1.46971\;53173\;1=T^{2}.$ ⓘ Symbols: $M$: limit and $T$: limit Permalink: http://dlmf.nist.gov/19.36.E10 Encodings: TeX, TeX, TeX, pMML, pMML, pMML, png, png, png See also: Annotations for §19.36(ii), §19.36(ii), §19.36 and Ch.19

Hence

 19.36.11 $R_{F}\left(1,2,4\right)=R_{C}\left(T^{2}+M^{2},T^{2}\right)=0.68508\;58166,$

in agreement with (19.36.5). Here $R_{C}$ is computed either by the duplication algorithm in Carlson (1995) or via (19.2.19).

 19.36.12 $\displaystyle h_{n}$ $\displaystyle=\sqrt{t_{n}^{2}+\theta a_{n}^{2}},$ $n=0,1,2,\dots$, $\displaystyle h_{n}$ $\displaystyle=h_{n-1}\frac{t_{n}}{\sqrt{t_{n}^{2}+\theta c_{n}^{2}}},$ $n=1,2,\dots$. ⓘ Symbols: $n$: nonnegative integer and $h_{n}$ Referenced by: §19.36(ii) Permalink: http://dlmf.nist.gov/19.36.E12 Encodings: TeX, TeX, pMML, pMML, png, png See also: Annotations for §19.36(ii), §19.36(ii), §19.36 and Ch.19

Then

 19.36.13 $2R_{G}\left(t_{0}^{2},t_{0}^{2}+\theta c_{0}^{2},t_{0}^{2}+\theta a_{0}^{2}% \right)=\left(t_{0}^{2}+\theta\sum_{m=0}^{\infty}2^{m-1}c_{m}^{2}\right)R_{C}% \left(T^{2}+\theta M^{2},T^{2}\right)+h_{0}+\sum_{m=1}^{\infty}2^{m}(h_{m}-h_{% m-1}).$

If the iteration of (19.36.6) and (19.36.12) is stopped when $c_{s}<\sqrt{r}t_{s}$ ($M$ and $T$ being approximated by $a_{s}$ and $t_{s}$, and the infinite series being truncated), then the relative error in $R_{F}$ and $R_{G}$ is less than $r$ if we neglect terms of order $r^{2}$.

$F\left(\phi,k\right)$ can be evaluated by using (19.25.5). $E\left(\phi,k\right)$ can be evaluated by using (19.25.7), and $R_{D}$ by using (19.21.10), but cancellations may become significant. Thompson (1997, pp. 499, 504) uses descending Landen transformations for both $F\left(\phi,k\right)$ and $E\left(\phi,k\right)$. A summary for $F\left(\phi,k\right)$ is given in Gautschi (1975, §3). For computation of $K\left(k\right)$ and $E\left(k\right)$ with complex $k$ see Fettis and Caslin (1969) and Morita (1978).

(19.22.20) reduces to $0=0$ if $p=x$ or $p=y$, and (19.22.19) reduces to $0=0$ if $z=x$ or $z=y$. Near these points there will be loss of significant figures in the computation of $R_{J}$ or $R_{D}$.

Descending Gauss transformations of $\Pi\left(\phi,\alpha^{2},k\right)$ (see (19.8.20)) are used in Fettis (1965) to compute a large table (see §19.37(iii)). This method loses significant figures in $\rho$ if $\alpha^{2}$ and $k^{2}$ are nearly equal unless they are given exact values—as they can be for tables. If $\alpha^{2}=k^{2}$, then the method fails, but the function can be expressed by (19.6.13) in terms of $E\left(\phi,k\right)$, for which Neuman (1969b) uses ascending Landen transformations.

Computation of Legendre’s integrals of all three kinds by quadratic transformation is described by Cazenave (1969, pp. 128–159, 208–230).

Quadratic transformations can be applied to compute Bulirsch’s integrals (§19.2(iii)). The function $\operatorname{cel}\left(k_{c},p,a,b\right)$ is computed by successive Bartky transformations (Bulirsch and Stoer (1968), Bulirsch (1969b)). The function $\operatorname{el2}\left(x,k_{c},a,b\right)$ is computed by descending Landen transformations if $x$ is real, or by descending Gauss transformations if $x$ is complex (Bulirsch (1965b)). Remedies for cancellation when $x$ is real and near 0 are supplied in Midy (1975). See also Bulirsch (1969a) and Reinsch and Raab (2000).

Bulirsch (1969a, b) extend Bartky’s transformation to $\operatorname{el3}\left(x,k_{c},p\right)$ by expressing it in terms of the first incomplete integral, a complete integral of the third kind, and a more complicated integral to which Bartky’s method can be applied. The cases $k_{c}^{2}/2\leq p<\infty$ and $-\infty require different treatment for numerical purposes, and again precautions are needed to avoid cancellations.

## §19.36(iii) Via Theta Functions

Lee (1990) compares the use of theta functions for computation of $K\left(k\right)$, $E\left(k\right)$, and $K\left(k\right)-E\left(k\right)$, $0\leq k^{2}\leq 1$, with four other methods. Also, see Todd (1975) for a special case of $K\left(k\right)$. For computation of Legendre’s integral of the third kind, see Abramowitz and Stegun (1964, §§17.7 and 17.8, Examples 15, 17, 19, and 20). For integrals of the second and third kinds see Lawden (1989, §§3.4–3.7).

## §19.36(iv) Other Methods

Numerical quadrature is slower than most methods for the standard integrals but can be useful for elliptic integrals that have complicated representations in terms of standard integrals. See §3.5.

For series expansions of Legendre’s integrals see §19.5. Faster convergence of power series for $K\left(k\right)$ and $E\left(k\right)$ can be achieved by using (19.5.1) and (19.5.2) in the right-hand sides of (19.8.12). A three-part computational procedure for $\Pi\left(\phi,\alpha^{2},k\right)$ is described by Franke (1965) for $\alpha^{2}<1$.

When the values of complete integrals are known, addition theorems with $\psi=\pi/2$19.11(ii)) ease the computation of functions such as $F\left(\phi,k\right)$ when $\frac{1}{2}\pi-\phi$ is small and positive. Similarly, §19.26(ii) eases the computation of functions such as $R_{F}\left(x,y,z\right)$ when $x$ ($>0$) is small compared with $\min(y,z)$. These special theorems are also useful for checking computer codes.

For fast methods for computing the incomplete elliptic integral of the first kind see Karp and Sitnik (2007) and Fukushima (2010).