Digital Library of Mathematical Functions
About the Project
19 Elliptic IntegralsComputation

§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 RF the polynomial of degree 7, for example, is

19.36.1 1-110E2+114E3+124E22-344E2E3-5208E23+3104E32+116E22E3,

where the elementary symmetric functions Es 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 RD and RJ 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-314E2+16E3+988E22-322E4-952E2E3+326E5-116E23+340E32+320E2E4+45272E22E3-968(E3E4+E2E5).

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

The duplication method starts with computation of λ0=xy+yz+zx. The reductions in §19.29(i) represent x,y,z as squares, for example x=U122 in (19.29.4). Because U12 may be real and negative, or even complex, care is needed to ensure x=U12, and similarly for y and z. This precaution is needed only for λ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 RF(1,2,4)=RF(z1,z2,z3),

where, in the notation of (19.19.7) with a=-12 and n=3,

19.36.4 z1=2.10985 99098 8,z3=2.15673 49098 8,Z1=0.00977 77253 5,z2=2.12548 49098 8,A=2.13069 32432 1,Z2=0.00244 44313 4,Z3=-Z1-Z2=-0.01222 21566 9,E2=-1.25480 14×10⁻⁴, E3=-2.9212×10⁻⁷.

The first five terms of (19.36.1) suffice for

19.36.5 RF(1,2,4)=0.68508 58166.

All cases of RF, RC, RJ, and RD 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 RJ 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 RJ without computing RC more than once. Because of cancellations in (19.26.21) it is advisable to compute RG from RF and RD 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 0k21, cancellations may lead to loss of significant figures when k2 is close to 1 and ϕ>π/4, as shown by Reinsch and Raab (2000). The cancellations can be eliminated, however, by using (19.25.10).

Accurate values of F(ϕ,k)-E(ϕ,k) for k2 near 0 can be obtained from RD by (19.2.6) and (19.25.13).

§19.36(ii) Quadratic Transformations

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 RF(x,y,z) and RG(x,y,z) can be computed by successive transformations in which two of the three variables converge quadratically to a common value and the integrals reduce to RC, accompanied by two quadratically convergent series in the case of RG; compare Carlson (1965, §§5,6). (In Legendre’s notation the modulus k approaches 0 or 1.) Let

19.36.6 2an+1 =an+an2-cn2,
2cn+1 =an-an2-cn2
2tn+1 =tn+tn2+θcn2,

where n=0,1,2,, and

19.36.7 0 <c0
t0 0,
t02+θa02 0,
θ =±1.

Then (19.22.18) implies that

19.36.8 RF(tn2,tn2+θcn2,tn2+θan2)

is independent of n. As n, cn, an, and tn converge quadratically to limits 0, M, and T, respectively; hence

19.36.9 RF(t02,t02+θc02,t02+θa02)=RF(T2,T2,T2+θM2)=RC(T2+θM2,T2).

If t0=a0 and θ=-1, so that tn=an, 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 θ=1 (leading ultimately to a hyperbolic case of RC) or a descending Gauss transformation if θ=-1 (leading to a circular case of RC). If x, y, and z are permuted so that 0x<y<z, then the computation of RF(x,y,z) is fastest if we make c02a02/2 by choosing θ=1 when y<(x+z)/2 or θ=-1 when y(x+z)/2.


We compute RF(1,2,4) by setting θ=1, t0=c0=1, and a0=3. Then

19.36.10 c32 =6.65×10⁻¹²,
a32 =2.46209 30206 0
t32 =1.46971 53173 1


19.36.11 RF(1,2,4)=RC(T2+M2,T2)=0.68508 58166,

in agreement with (19.36.5). Here RC is computed either by the duplication algorithm in Carlson (1995) or via (19.2.19).

To (19.36.6) add

19.36.12 hn =tn2+θan2,
hn =hn-1tntn2+θcn2,


19.36.13 2RG(t02,t02+θc02,t02+θa02)=(t02+θm=02m-1cm2)RC(T2+θM2,T2)+h0+m=12m(hm-hm-1).

If the iteration of (19.36.6) and (19.36.12) is stopped when cs<rts (M and T being approximated by as and ts, and the infinite series being truncated), then the relative error in RF and RG is less than r if we neglect terms of order r2.

F(ϕ,k) can be evaluated by using (19.25.5). E(ϕ,k) can be evaluated by using (19.25.7), and RD by using (19.21.10), but cancellations may become significant. Thompson (1997, pp. 499, 504) uses descending Landen transformations for both F(ϕ,k) and E(ϕ,k). A summary for F(ϕ,k) is given in Gautschi (1975, §3). For computation of K(k) and E(k) 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 RJ or RD.

Descending Gauss transformations of Π(ϕ,α2,k) (see (19.8.20)) are used in Fettis (1965) to compute a large table (see §19.37(iii)). This method loses significant figures in ρ if α2 and k2 are nearly equal unless they are given exact values—as they can be for tables. If α2=k2, then the method fails, but the function can be expressed by (19.6.13) in terms of E(ϕ,k), 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 cel(kc,p,a,b) is computed by successive Bartky transformations (Bulirsch and Stoer (1968), Bulirsch (1969b)). The function el2(x,kc,a,b) 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 el3(x,kc,p) 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 kc2/2p< and -<p<kc2/2 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(k), E(k), and K(k)-E(k), 0k21, with four other methods. Also, see Todd (1975) for a special case of K(k). 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(k) and E(k) 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 Π(ϕ,α2,k) is described by Franke (1965) for α2<1.

When the values of complete integrals are known, addition theorems with ψ=π/219.11(ii)) ease the computation of functions such as F(ϕ,k) when 12π-ϕ is small and positive. Similarly, §19.26(ii) eases the computation of functions such as RF(x,y,z) 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).