Three Improvements in Reduction and Computation of Elliptic Integrals

Three improvements in reduction and computation of elliptic integrals are made. 1. Reduction formulas, used to express many elliptic integrals in terms of a few standard integrals, are simplified by modifying the definition of intermediate “basic integrals.” 2. A faster than quadratically convergent series is given for numerical computation of the complete symmetric elliptic integral of the third kind. 3. A series expansion of an elliptic or hyperelliptic integral in elementary symmetric functions is given, illustrated with numerical coefficients for terms through degree seven for the symmetric elliptic integral of the first kind. Its usefulness for elliptic integrals, in particular, is important.


Foreword
Elliptic integrals have many applications, for example in mathematics and physics: • arclengths of plane curves (ellipse, hyperbola, Bernoulli's lemniscate) • surface area of an ellipsoid in 3-dimensional space • electric and magnetic fields associated with ellipsoids • periodicity of anharmonic oscillators • mutual inductance of coaxial circles • age of the universe in a Friedmann model These applications are mentioned in the chapter on elliptic integrals, written by B. C. Carlson, that will appear in the NIST Digital Library of Mathematical Functions .
The DLMF is scheduled to begin service in 2003 from a NIST Web site. A hardcover book will be published also. These resources will provide a complete guide to the higher mathematical functions for use by experienced scientific professionals. The book will provide mathematical formulas, references to proofs, references to extensions and generalizations, graphs, brief descriptions of computational methods, a survey of useful published tables, and sample applications. The Web site will include, in addition, interactive visualizations of 3-dimensional surfaces, a mathematics-aware search engine, a downloading capability for equations, live links to Web sites that provide mathematical software, and a limited facility for generating tables on demand.
The DLMF is modeled after the Handbook of Mathematical Functions , published in 1964 by the National Bureau of Standards with M. Abramowitz and I. A. Stegun as editors. This handbook has been enormously successful: it has sold more than 500,000 copies, its sales remain high, and it is very frequently cited in journal articles in physics and many other fields. But new properties of the higher functions have been developed, and new functions have risen in importance in applications, since the publication of the Abramowitz and Stegun handbook. More than half of the old handbook was devoted to tables, now made obsolete by the revolutionary improvements since 1964 in computers and software. The need for a modern reference is being filled by more than 30 expert authors who are working under contract to NIST and supervised by four NIST editors. The writing is being edited carefully to assure consistent style and level of content.
Elliptic integrals have long been associated with the name of Legendre. Legendre's incomplete elliptic integrals are The complete forms of these integrals are obtained by setting = /2. Over a period of more than 35 years, Carlson has published a series of papers that provide valuable new mathematical and computational foundations for the subject in terms of the symmetric elliptic integrals , The complete forms are obtained by setting x = 0. In comparison with Legendre's integrals, Carlson's integrals simplify the reduction of general elliptic integrals to standard forms and open the way to efficient computations by application of a duplication theorem.
One of the purposes of the DLMF project is to stimulate research into the theory, computation and application of the higher mathematical functions. The paper which follows is an example. It is a further development of material that appears in Carlson's DLMF chapter on elliptic integrals.

Simplified Formulas for Reducing Elliptic Integrals
A large class of elliptic integrals can be written in the form where m = (m 1 , . . ., m n ) is an n -tuple of integers (positive, negative or zero), x > y , h = 3 or 4, n Ն h , and the different linear factors are not proportional. The a 's and b 's may be complex (with the b 's not equal to zero), but the integral is assumed to be well defined, possibly as a Cauchy principal value. In particular the line segment with endpoints a i + b i x and a i + b i y is assumed to lie in the cut plane ( C\(Ϫϱ, 0) for 1 Յ i Յ h . We write m = ͚ n j=1 m j e j , where e j is an n -tuple with 1 in the j th position and 0's elsewhere, and we define 0 = (0, . . ., 0). Reference [1] gives a general method of reducing I (m) to a linear combination of "basic integrals," defined as The b 's appear also in the other two formulas and therefore in all reduction formulas, sometimes in great profusion if m is considerably more complicated than in Eq.
(1.2). It has been found that the b 's disappear from all three formulas, and therefore from all reduction formulas, if we define Here A (m) is the algebraic function The other recurrence relation, Ref. [1], Eq. (3.5), becomes where E hϪr is the elementary symmetric function defined by where M = ͚ n j=1 m j and each sum is empty if its upper limit is less than its lower limit. The first term on the right is independent of k , which is usually best chosen so that 1 Յ k Յ h . The coefficients are defined by

Algorithms for Complete Elliptic Integrals of the Third Kind
Complications formerly encountered in numerical computation of Legendre's complete elliptic integral of the third kind were avoided by defining and tabulating Heuman's Lambda function (for circular cases) and a modification of Jacobi's Zeta function (for hyperbolic cases). For example, the method of Ref. [3] was later superseded by Bartky's transformation and its application by Bulirsch [4] to his complete integral cel (k c , p , a , b ). Bartky's transformation for the complete case of the symmetric integral of the third kind, obtained from Ref. [5], Eq. (4.14) with the help of (3/4)R L (y , z , p ) = 3R F (0, y , z ) Ϫ pR J (0, y , z , p ), can be written as where a n , g n , p n are positive, s n = (p 4 n Ϫ a 2 n g 2 n )/8p 4 n , and a n+1 = a n + g n 2 , g n+1 = ͙a n g n , p n+1 = p 2 n + a n g n 2p n , n ʦ ‫.ގ‬ (2.3) As n → ϱ, a n and g n converge quadratically to Gauss's arithmetic-geometric mean, M (a 0 , g 0 ), and R F (0, g 2 n , a 2 n ) = /2M (a 0 , g 0 ), n ʦ ‫,ގ‬ Since g n+1 converges quadratically to M (a 0 , g 0 ), so does p n , and s n converges quadratically to 0. Iteration of Eq.
(2.2) gives where and therefore Q n+1 Q n = s n p 2 n p 2 n+1 = 1 2 p 2 n Ϫ a n g n p 2 n + a n g n . (2.8) Letting n → ϱ in Eq. (2.6), we find, for positive a 0 , g 0 , p 0 , n Ϫ a n g n p 2 n + a n g n , (2.9) where ⑀ n converges to 0 quadratically and Q n converges to 0 faster than quadratically. If p 0 = a 0 , Eq. (2.9) becomes Q n ⑀ n , ⑀ n = a n Ϫ g n a n + g n , (2.10) where R D is a complete integral of the second kind, symmetric in only its first two arguments. If 0 < a 0 Յ g 0 then Ϫ1 < ⑀ 0 Յ 0, but ⑀ n Ն 0 and Q n Յ 0 for n Ն 1.

Expansion in Elementary Symmetric Functions
The duplication method of computing the symmetric elliptic integrals R F and R J (including their degenerate cases R C and R D ) consists in iterating their duplication theorems until their variables are nearly equal and then expanding in a series of elementary symmetric functions of the small differences between the variables. In the absence of a duplication theorem the method is useful for a hyperelliptic integral only if the variables are nearly equal. The series is truncated to a polynomial of fixed degree; the higher the degree, the fewer duplications are needed for a desired accuracy of the result but the larger the number of terms to be calculated. No tests have been made to determine an optimal compromise, which would depend in large part on the speed of extracting square roots in the duplication theorem and therefore on the equipment used. In Ref. [8] polynomials of degree five were chosen for simplicity, but later it seemed worthwhile to increase the degree for the comparatively slow computation of R J , and Ref. [9], Eq. (A.11) gives the terms of degree six and seven for R J . The change in speed would be significant only if a very large number of computations were performed, since the result of a single computation is returned with no delay apparent to the eye. We shall give here the corresponding terms for R F and the general form of the infinite series of which the polynomials are truncations.
Any R -function R Ϫa (b 1 , . . ., b k ; z 1 , . . ., z k ) (see Ref. [6]) in which all the b -parameters are positive integral multiples of a single number ␤ can be rewritten with repeated variables and all b 's equal to ␤ . An example is and R D (x , y , z ) is the case with p = z . Therefore we consider . ., ␤ ; z 1 /A , . . ., z n /A ), (3.2) where B is the beta function and we have used the homogeneity of R to divide the variables by their arithmetic average, The relative difference between A and z i is Each application of a duplication theorem decreases by a factor of four the differences between the variables, therefore the difference A Ϫ z i , and ultimately |Z | as A approaches a limit. It is easy to determine whether another duplication is needed before using the truncated series to achieve the desired accuracy. If a = ␤ = 1/2 and n = 3, the series [Eq. (3.8)] with E 1 = 0 can be rearranged as