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

where the elementary symmetric functions are defined by (19.19.4). If (19.36.1) is used instead of its first five terms, then the factor in Carlson (1995, (2.2)) is changed to .

For both and the factor in Carlson (1995, (2.18)) is changed to when the following polynomial of degree 7 (the same for both) is used instead of its first seven terms:

19.36.2

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

The duplication method starts with computation of . The reductions in §19.29(i) represent as squares, for example in (19.29.4). Because may be real and negative, or even complex, care is needed to ensure , and similarly for and . This precaution is needed only for . 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

where, in the notation of (19.19.7) with and ,

The first five terms of (19.36.1) suffice for

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

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

19.36.6

where , and

19.36.7

Then (19.22.18) implies that

is independent of . As , , , and converge quadratically to limits 0, , and , respectively; hence

If and , so that , then this procedure reduces to the AGM method for the complete integral.

The step from to is an ascending Landen transformation if (leading ultimately to a hyperbolic case of ) or a descending Gauss transformation if (leading to a circular case of ). If , , and are permuted so that , then the computation of is fastest if we make by choosing when or when .

### ¶ Example

We compute by setting , , and . Then

19.36.10

Hence

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

To (19.36.6) add

19.36.12
,
.

Then

If the iteration of (19.36.6) and (19.36.12) is stopped when ( and being approximated by and , and the infinite series being truncated), then the relative error in and is less than if we neglect terms of order .

can be evaluated by using (19.25.5). can be evaluated by using (19.25.7), and by using (19.21.10), but cancellations may become significant. Thompson (1997, pp. 499, 504) uses descending Landen transformations for both and . A summary for is given in Gautschi (1975, §3). For computation of and with complex see Fettis and Caslin (1969) and Morita (1978).

(19.22.20) reduces to if or , and (19.22.19) reduces to if or . Near these points there will be loss of significant figures in the computation of or .

Descending Gauss transformations of (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 and are nearly equal unless they are given exact values—as they can be for tables. If , then the method fails, but the function can be expressed by (19.6.13) in terms of , 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 is computed by successive Bartky transformations (Bulirsch and Stoer (1968), Bulirsch (1969b)). The function is computed by descending Landen transformations if is real, or by descending Gauss transformations if is complex (Bulirsch (1965b)). Remedies for cancellation when 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 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.

## §19.36(iii) Via Theta Functions

Lee (1990) compares the use of theta functions for computation of , , and , , with four other methods. Also, see Todd (1975) for a special case of . 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 and 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 is described by Franke (1965) for .

When the values of complete integrals are known, addition theorems with 19.11(ii)) ease the computation of functions such as when is small and positive. Similarly, §19.26(ii) eases the computation of functions such as when () is small compared with . 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).