- §19.36(i) Duplication Method
- §19.36(ii) Quadratic Transformations
- §19.36(iii) Via Theta Functions
- §19.36(iv) Other Methods

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-\frac{1}{10}{E}_{2}+\frac{1}{14}{E}_{3}+\frac{1}{24}{E}_{2}^{2}-\frac{3}{44}{E}_{2}{E}_{3}-\frac{5}{208}{E}_{2}^{3}+\frac{3}{104}{E}_{3}^{2}+\frac{1}{16}{E}_{2}^{2}{E}_{3},$$ | ||

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 ${\left(3r\right)}^{-1/6}$ in Carlson (1995, (2.2)) is changed to ${\left(3r\right)}^{-1/8}$.

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

19.36.2 | $$\begin{array}{l}1-\frac{3}{14}{E}_{2}+\frac{1}{6}{E}_{3}+\frac{9}{88}{E}_{2}^{2}-\frac{3}{22}{E}_{4}-\frac{9}{52}{E}_{2}{E}_{3}+\frac{3}{26}{E}_{5}-\frac{1}{16}{E}_{2}^{3}\\ \phantom{\rule{2em}{0ex}}+\frac{3}{40}{E}_{3}^{2}+\frac{3}{20}{E}_{2}{E}_{4}+\frac{45}{272}{E}_{2}^{2}{E}_{3}-\frac{9}{68}\left({E}_{3}{E}_{4}+{E}_{2}{E}_{5}\right).\end{array}$$ | ||

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.

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),$$ | ||

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

19.36.4 | $$\begin{array}{l}\begin{array}{cc}\hfill {z}_{1}& =\mathrm{2.10985\hspace{0.33em}99098\hspace{0.33em}8},\hfill \\ \hfill {z}_{3}& =\mathrm{2.15673\hspace{0.33em}49098\hspace{0.33em}8},\hfill \\ \hfill {Z}_{1}& =\mathrm{0.00977\hspace{0.33em}77253\hspace{0.33em}5},\hfill \end{array}\hspace{1em}\begin{array}{cc}\hfill {z}_{2}& =\mathrm{2.12548\hspace{0.33em}49098\hspace{0.33em}8},\hfill \\ \hfill A& =\mathrm{2.13069\hspace{0.33em}32432\hspace{0.33em}1},\hfill \\ \hfill {Z}_{2}& =\mathrm{0.00244\hspace{0.33em}44313\hspace{0.33em}4},\hfill \end{array}{Z}_{3}=-{Z}_{1}-{Z}_{2}=-\mathrm{0.01222\hspace{0.33em}21566\hspace{0.33em}9},\\ \phantom{\rule{0em}{0ex}}{E}_{2}=\mathit{-1.25480\hspace{0.33em}14\times 10\u207b\u2074},\hspace{1em}{E}_{3}=\mathrm{-2.9212\times 10\u207b\u2077.}\end{array}$$ | ||

The first five terms of (19.36.1) suffice for

19.36.5 | $${R}_{F}\left(1,2,4\right)=\mathrm{0.68508\hspace{0.33em}58166}\mathrm{\dots}.$$ | ||

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\le {k}^{2}\le 1$, cancellations may lead to loss of significant figures when ${k}^{2}$ is close to 1 and $\varphi >\pi /4$, as shown by Reinsch and Raab (2000). The cancellations can be eliminated, however, by using (19.25.10).

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 | $2{a}_{n+1}$ | $={a}_{n}+\sqrt{{a}_{n}^{2}-{c}_{n}^{2}},$ | ||

$2{c}_{n+1}$ | $={a}_{n}-\sqrt{{a}_{n}^{2}-{c}_{n}^{2}}={c}_{n}^{2}/\left(2{a}_{n+1}\right),$ | |||

$2{t}_{n+1}$ | $={t}_{n}+\sqrt{{t}_{n}^{2}+\theta {c}_{n}^{2}},$ | |||

where $n=0,1,2,\mathrm{\dots}$, and

19.36.7 | $0$ | $$ | ||

${t}_{0}$ | $\ge 0,$ | |||

${t}_{0}^{2}+\theta {a}_{0}^{2}$ | $\ge 0,$ | |||

$\theta $ | $=\pm 1.$ | |||

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)$$ | ||

is independent of $n$. As $n\to \mathrm{\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 $$, then the computation of ${R}_{F}\left(x,y,z\right)$ is fastest if we make ${c}_{0}^{2}\le {a}_{0}^{2}/2$ by choosing $\theta =1$ when $$ or $\theta =-1$ when $y\ge \left(x+z\right)/2$.

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 | ${c}_{3}^{2}$ | $=\mathit{6.65\times 10\u207b\xb9\xb2},$ | ||

${a}_{3}^{2}$ | $=\mathrm{2.46209\hspace{0.33em}30206\hspace{0.33em}0}={M}^{2},$ | |||

${t}_{3}^{2}$ | $=\mathrm{1.46971\hspace{0.33em}53173\hspace{0.33em}1}={T}^{2}.$ | |||

Hence

19.36.11 | $${R}_{F}\left(1,2,4\right)={R}_{C}\left({T}^{2}+{M}^{2},{T}^{2}\right)=\mathrm{0.68508\hspace{0.33em}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).

To (19.36.6) add

19.36.12 | ${h}_{n}$ | $=\sqrt{{t}_{n}^{2}+\theta {a}_{n}^{2}},$ | ||

$n=0,1,2,\mathrm{\dots}$, | ||||

${h}_{n}$ | $={h}_{n-1}{\displaystyle \frac{{t}_{n}}{\sqrt{{t}_{n}^{2}+\theta {c}_{n}^{2}}}},$ | |||

$n=1,2,\mathrm{\dots}$. | ||||

Then

19.36.13 | $$\begin{array}{l}2{R}_{G}\left({t}_{0}^{2},{t}_{0}^{2}+\theta {c}_{0}^{2},{t}_{0}^{2}+\theta {a}_{0}^{2}\right)\\ \phantom{\rule{2em}{0ex}}=\left({t}_{0}^{2}+\theta \sum _{m=0}^{\mathrm{\infty}}{2}^{m-1}{c}_{m}^{2}\right){R}_{C}\left({T}^{2}+\theta {M}^{2},{T}^{2}\right)+{h}_{0}+\sum _{m=1}^{\mathrm{\infty}}{2}^{m}\left({h}_{m}-{h}_{m-1}\right).\end{array}$$ | ||

If the iteration of (19.36.6) and (19.36.12) is stopped when $$ ($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(\varphi ,k\right)$ can be evaluated by using (19.25.5). $E\left(\varphi ,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(\varphi ,k\right)$ and $E\left(\varphi ,k\right)$. A summary for $F\left(\varphi ,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 $\mathrm{\Pi}\left(\varphi ,{\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(\varphi ,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 $\mathrm{cel}\left({k}_{c},p,a,b\right)$ is computed by successive Bartky transformations (Bulirsch and Stoer (1968), Bulirsch (1969b)). The function $\mathrm{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 $\mathrm{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 $$ and $$ require different treatment for numerical purposes, and again precautions are needed to avoid cancellations.

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\le {k}^{2}\le 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).

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 $\mathrm{\Pi}\left(\varphi ,{\alpha}^{2},k\right)$ is described by Franke (1965) for $$.

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(\varphi ,k\right)$ when $\frac{1}{2}\pi -\varphi $ 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\left(y,z\right)$. These special theorems are also useful for checking computer codes.