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.

Key words: computational algorithm; elementary symmetric function; elliptic integral; hyperelliptic integral; hypergeometric R-function; recurrence relations; series expansion.

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

F([phi], k) = [[integral].sup.[phi].sub.0] d[theta]/[square root of (1 - [k.sup.2] [sin.sup.2] [theta])]'

E([phi], k) = [[integral].sup.[phi].sub.0] [square root of (1 - [k.sup.2] [sin.sup.2] [theta] d[theta])],

and

[PI]([phi], [[alpha].sup.2], k) = [[integral].sup.[phi].sub.0] d[theta]/[square root of (1 - [k.sup.2] [sin.sup.2] [theta])] (1 - [[alpha].sup.2] [sin.sup.2] [theta]).

The complete forms of these integrals are obtained by setting [phi] = [pi]/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

[R.sub.F](x, y, z) = 1/2 [[integral].sup.[infinity].sub.0] dt/s(t),

[R.sub.D](x, y, z) = 3/2 [[integral].sup.[infinity].sub.0] dt/s(t)(t + z),

[R.sub.J](x, y, z, p) = 3/2 [[integral].sup.[infinity].sub.0] dt/s(t)(t + p)

where

s(t) = [square root of (t + x)] [square root of (t + y)] [square root of (t + z)].

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.

Daniel W. Lozier

NIST Mathematical and Computational Sciences Division

1. Simplified Formulas for Reducing Elliptic Integrals

A large class of elliptic integrals can be written in the form

I(m) = [[integral].sup.x.sub.y] [[PI].sup.h.sub.i=1][([a.sub.i] + [b.sub.i]t)].sup.-1/2] [[PI].sup.n.sub.j=1][([a.sub.j] + [b.sub.j]t).sup.m.sub.j] dt, (1.1)

where m = ([m.sub.1], ..., [m.sub.n] is an n-tuple of integers (positive, negative or zero), x > y, h = 3 or 4, n [greater than or equal to] 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.sub.i], + [b.sub.i]x and [a.sub.i] + [b.sub.i]y is assumed to lie in the cut plane C\(-[infinity], 0) for 1 [less than or equal to] [less than or equal to] h.

We write m = [summation over (n/j=1)] [m.sub.j][e.sub.j], where [e.sub.j] is an n-tuple with 1 in the jth 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 I(0), I([-e.sub.j]), 1[less than or equal to] j [less than or equal to] n, and (if h = 4) I([e.sub.k]), 1 [less than or equal to] k [less than or equal to] 4. A simple example of such a reduction is

[b.sub.i]I([e.sub.j] -- [e.sub.i]) = [d.sub.ji] I([--e.sub.i]) + [b.sub.j] I(0),

i,j [member of] {1, ..., n}, (1.2)

where [d.sub.ji] = [a.sub.j][b.sub.i] - [a.sub.i][b.sub.j]. This equation reduces all 36+72 integrals in Ref. (2), Eqs. (3.142) and (3.168) and also the 18 integrals in Ref. (2), Eq. (3.159) after taking [x.sup.2] as a new variable of integration. The basic integrals are expressed in terms of symmetric standard integrals [R.sub.F], [R.sub.D], and [R.sub.J] in Ref. (1), Sec. 4.

The general method first reduces I(m) by Ref. (1), Eq. (2.19) to integrals in which m has at most one nonzero component and then uses two recurrence relations Ref. (1) Eqs. (3.5) and (3.11) for further reduction to basic integrals. For example, the simplest recurrence relation is Ref. (1), Eq. (3.11):

[b.sup.q.sub.i]I(q[e.sub.j]) = [summation over (q/r=0)] (q/r) [b.sup.r.sub.j] [d.sup.q-r.sub.ji] I(r[e.sub.i]), q [member of] N. (1.3)

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

I(m) = I(m)/B, A(m) = A(m)/B,

B = [[PI].sup.n.sub.j=1][b.sup.mj.sub.j], [r.sub.ij] = [d.sub.ij]/[b.sub.i][b.sub.j] = [a.sub.i]/[b.sub.i] - [a.sub.j]/[b.sub.j]. (1.4)

Here A(m) is the algebraic function 

A(m) = [f.sub.m](x) - [f.sub.m](y), (1.5)

where [f.sub.m](t) is the integrand of Eq. (1.1). Note that I(0) = I(0).

For example, Eq. (1.2) becomes

I([e.sub.j] - [e.sub.i]) = [r.sub.ji] I(-[e.sub.i]) + I(0), (1.6)

and Eq. (1.3) becomes

I(q[e.sub.j]) = [summation over (q/r=0)] (q)(r) [r.sup.q-r.sub.ji] I(r[e.sub.i]), q [member of] N. (1.7)

The other recurrence relation, Ref. (1), Eq. (3.5), becomes

[summation over (h/r=0)](2q + r)[E.sub.h-r]([r.sub.1j], ..., [r.sub.hj])I((q - 1 + r)[e.sub.j])

= 2A(q[e.sub.j] + [summation over (h/i=1)][e.sub.i]), q [member of] Z, 1 [less than or equal to] j [less than or equal to] n (1.8)

where [E.sub.h-r] is the elementary symmetric function defined by

[[PI].sup.h.sub.i=1](1 + t[r.sub.ij]) = [summation over (h/k=0)][t.sup.k][E.sub.k]([r.sub.1j], ..., [r.sub.hj]), (1.9)

whence [E.sub.h] = 0 if 1 [less than or equal to] j [less than or equal to] h because [r.sub.jj] = 0.

The remaining formula, Ref. (1), Eq. (2.19), becomes

I(m) = [summation over (M/q=0)][C.sub.M-q](k)I(q[e.sub.k]) + [summation over (n/i=1)][[DELTA].sub.i][summation over (-[m.sub.i]/q=1)][C.sub.[m.sub.i]+q](i)I(-q[e.sub.i]), (1.10)

where M = [summation over (n/j=1)] [m.sub.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 [less than or equal to] k [less than or equal to] h. The coefficients are defined by

[[DELTA].sub.i] = [[PI].sup.n.sub.j=1,j[not equal to]i][r.sup.[m.sub.j].sub.ji], [C.sub.0](i) = 1, [C.sub.[+ or -]s](i) = [SIGMA] [[eta].sup.[[alpha].sub.l]].sub.[+ or -]l](i)...[[eta].sup.[[alpha].sub.s].sub.[+ or -]s](i)/[[alpha].sub.l]!...[[alpha].sub.s]!, [[eta].sub.[+ or -]p](i) = -1/p [[summation over (n/j=1,j[not equal to]i)][m.sub.j][r.sup.[+ or -]p.sub.ij], p [greater than or equal to] 1, (1.11)

where upper (lower) signs go together and the first sum extends over all nonnegative integers [[alpha].sub.l], ..., [[alpha].sub.s] such that [[alpha].sub.1] + 2[[alpha].sub.2] + ... + s[[alpha].sub.s] = s. A recurrence relation for the coefficients is

s [C.sub.[+ or -]s](i) = [summation over (s/p=1)]p[[eta].sub.[+ or -]p](i) [C.sub.[+ or -](s-p)](i), s [greater than or equal to] 1. (1.12)

2. 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.sub.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[pi]/4)[R.sub.L](y, z, p) = 3[R.sub.F](0, y, z) - p[R.sub.J](0, y, z, p), ([pi]/2)[R.sub.K](y, z) = [R.sub.F](0, y, z), (2.1)

can be written as

[R.sub.J](0, [g.sup.2.sub.n], [a.sup.2.sub.n], [p.sup.2.sub.n]) = [S.sub.n][R.sub.J](0, [g.sup.2.sub.n+1], [a.sup.2.sub.n+1], [p.sup.2.sub.n+1]) + (3/2[p.sup.2.sub.n])[R.sub.F](0, [g.sup.2.sub.n+1], [a.sup.2.sub.n+1]), (2.2)

where [a.sub.n], [g.sub.n], [p.sub.n] are positive, [S.sub.n] = ([p.sup.4.sub.n] - [a.sup.2.sub.n][g.sup.2.sub.n])/8[p.sup.4.sub.n], and

[a.sub.n+1] = [a.sub.n] + [g.sub.n]/2, [g.sub.n+1] = [square root of ([a.sub.n][g.sub.n])], [p.sub.n+1] = [p.sup.2.sub.n] + [a.sub.n][g.sub.n]/2[p.sub.n],

n [member of] N. (2.3)

As n [right arrow] [infinity], [a.sub.n] and [g.sub.n] converge quadratically to Gauss's arithmetic-geometric mean, M([a.sub.0], [g.sub.0]), and

[R.sub.F](0, [g.sup.2.sub.n], [a.sup.2.sub.n]) = [pi]/2M([a.sub.0], [g.sub.0]), n [member of] N, (2.4)

by Ref. (6), Eqs. (6.10-8) and Eq. (2.1). It follows from

(2.3) that

[p.sub.n+1] - [g.sub.n+1] = [([p.sub.n] - [g.sub.n+1]).sup.2]/2[p.sub.n]. (2.5)

Since [g.sub.n+1] converges quadratically to M ([a.sub.0], [g.sub.0]), so does [p.sub.n] and [S.sub.n] converges quadratically to 0. Iteration of Eq. (2.2) gives

[R.sub.J](0, [g.sup.2.sub.0], [a.sup.2.sub.0], [p.sup.2.sub.0]) = [Q.sub.n] [p.sup.2.sub.n]/[p.sup.2.sub.0] [R.sub.J](0, [g.sup.2.sub.n], [a.sup.2.sub.n], [p.sup.2.sub.n])

+ 3[pi]/4[p.sup.2.sub.0]M([a.sub.0], [g.sub.0]) [summation over (n-1/m=0)] [Q.sub.m], (2.6)

where

[Q.sub.0] = 1, [Q.sub.m]/[p.sup.2.sub.0] = [S.sub.0][S.sub.1]...[S.sub.m-1]/[p.sup.2.sub.m], m [greater than or equal to] 1, (2.7)

and therefore

[Q.sub.n+1]/[Q.sub.n] = [S.sub.n][p.sup.2.sub.n]/[p.sup.2.sub.n+1] = 1/2 [p.sup.2.sub.n] - [a.sub.n][g.sub.n]/[p.sup.2.sub.n] + [a.sub.n][g.sub.n]. (2.8)

Letting n [right arrow] [infinity] in Eq. (2.6), we find, for positive [a.sub.0], [g.sub.0], [p.sub.0],

[R.sub.J](0, [g.sup.2.sub.0], [a.sup.2.sub.0], [p.sup.2.sub.0]) = 3[pi]/4[p.sup.2.sub.0]M([a.sub.0], [g.sub.0])[summation over ([infinity]/n=0)] [Q.sub.n],

[Q.sub.0] = 1, [Q.sub.n+1] = 1/2 [Q.sub.n] [[epsilon].sub.n], [[epsilon].sub.n] = [p.sup.2.sub.n] - [a.sub.n][g.sub.n]/[p.sup.2.sub.n] + [a.sub.n][g.sub.n], (2.9)

where [[epsilon].sub.n] converges to 0 quadratically and [Q.sub.n] converges to 0 faster than quadratically.

If [p.sub.0] = [a.sub.0], Eq. (2.9) becomes

[R.sub.D](0, [g.sup.2.sub.0], [a.sup.2.sub.0] = 3[pi]/4[a.sup.2.sub.0]M([a.sub.0], [g.sub.0]) [summation over ([infinity]/n=0)] [Q.sub.n],

[Q.sub.0] = 1, [Q.sub.n+1] = 1/2 [Q.sub.n] [[epsilon].sub.n], [[epsilon].sub.n] = [a.sub.n] - [g.sub.n]/[a.sub.n] + [g.sub.n], (2.10)

where [R.sub.D] is a complete integral of the second kind, symmetric in only its first two arguments. If 0 < [a.sub.0] [less than or equal to] [g.sub.0] then -1 < [[epsilon].sub.0] [less than or equal to] 0, but [[epsilon].sub.n] [greater than or equal to] 0 and [Q.sub.n] [less than or equal to] 0 for n [greater than or equal to] 1.

If the last variable of [R.sub.J] is negative, the Cauchy principal value is given by

([q.sup.2.sub.0] + [a.sup.2.sub.0])[R.sub.J](0, [g.sup.2.sub.0], [a.sup.2.sub.0], - [q.sup.2.sub.0]) = ([p.sup.2.sub.0] - [a.sup.2.sub.0])[R.sub.J](0, [g.sup.2.sub.0], [a.sup.2.sub.0], [p.sup.2.sub.0])

- 3[pi]/2M([a.sub.0], [g.sub.0]), [p.sup.2.sub.0] = [a.sup.2.sub.0] ([q.sup.2.sub.0] + [g.sup.2.sub.0])/[q.sup.2.sub.0] + [a.sup.2.sub.0]), (2.11)

where we have used Eq. (2.4) and chosen [x.sub.i] = [a.sup.2.sub.0] in Ref. [7], Eq. (4.6). Substitution of Eq. (2.9) gives

[R.sub.J](0, [g.sup.2.sub.0], [a.sup.2.sub.0], - [q.sup.2.sub.0]) = - 3[pi]/4M([a.sub.0], [g.sub.0])([q.sup.2.sub.0] + [a.sup.2.sub.0])

(2 + [a.sup.2.sub.0] - [g.sup.2.sub.0]/[q.sup.2.sub.0] + [g.sup.2.sub.0] [summation over ([infinity]/n=0)] [Q.sub.n]). (2.12)

Equation (2.3) and the second line of Eq. (2.9) still apply, with [p.sub.0] given by Eq. (2.11).

For the complete case of Legendre's third integral,

[PI]([[alpha].sup.2], k) = ([[alpha].sup.2]/3)[R.sub.J](0, [k'.sup.2], 1, 1 - [[alpha].sup.2]) + [R.sub.F](0, [k'.sup.2], 1), (2.13)

where [k'.sup.2] = 1 - [k.sup.2], Eq. (2.9) implies

[PI]([[alpha].sup.2], k) = [pi]/4M(1, k') (2 + [[alpha].sup.2]/1 - [[alpha].sup.2] [summation over ([infinity]/n=0)] [Q.sub.n]),

- [infinity] < [k.sup.2] < 1, - [infinity] < [[alpha].sup.2] < 1, (2.14)

Where Eq. (2.3) applies with [a.sub.0] = 1, [g.sub.0] = k', and [p.sup.2.sub.0] = 1 - [[alpha].sup.2].

If [[alpha].sup.2] > 1, Eq. (2.12) gives the Cauchy principal value,

[PI]([[alpha].sup.2], k) = -[pi][k.sup.2]/4M(1, k')([[alpha].sup.2] - [k.sup.2]) [summation over ([infinity]/n=0)] [Q.sub.n],

- [infinity] < [k.sup.2] < 1, - 1 < [[alpha].sup.2] < [infinity], (2.15)

where Eq. (2.3) and the second line of Eq. (2.9) apply with [a.sub.0] = 1, [g.sub.0] = k', and [p.sup.2.sub.0] = 1 - [k.sup.2]/[[alpha].sup.2].

3. Expansion in Elementary Symmetric Functions

The duplication method of computing the symmetric elliptic integrals [R.sub.F] and [R.sub.J] (including their degenerate cases [R.sub.C] and [R.sub.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 comparative ly slow computation of [R.sub.J], and Ref. (9), Eq. (A.11) gives the terms of degree six and seven for [R.sub.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.sub.F] and the general form of the infinite series of which the polynomials are truncations.

Any R-function [R.sub.-a]([b.sub.1], ..., [b.sub.k]; [z.sub.1], ..., [z.sub.k]) (see Ref. (6)) in which all the b-parameters are positive integral multiples of a single number [beta] can be rewritten with repeated variables and all b's equal to [beta]. An example is

[R.sub.J](x,y,z,p) = [R.sub.-3/2](1/2,1/2,1/2, 1; x,y,z,p)

= [R.sub.-3/2](1/2,1/2,1/2,1/2,1/2,x,y,z,p,p), (3.1)

and [R.sub.D](x,y,z) is the case with p = z. Therefore we consider

[R.sub.-a]([beta], ..., [beta]; [z.sub.1], ..., [z.sub.n])

= 1/B(a, n[beta] - a) [[integral].sup.[infinity].sub.0] [t.sup.n[beta]-a-1] [[PI].sup.n.sub.i=1][(t + [z.sub.i]).sup.-[beta]]dt

= [A.sup.-a][R.sub.-a]([beta], ..., [beta]; [z.sub.1]/A, ..., [z.sub.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,

A = 1/n [summation over(n/i=1)] [z.sub.i]. (3.3)

The relative difference between A and [z.sub.i] is

[Z.sub.i] = A - [z.sub.i]/A = 1 - [z.sub.i]/A, (3.4)

whence

[R.sub.-a]([beta], ..., [beta]; [z.sub.1], ..., [z.sub.n])

= [A.sup.-a][R.sub.-a]([beta], ..., [beta]; 1 - [Z.sub.1], ..., 1 - [Z.sub.n]). (3.5)

Because the function is symmetric in the Z's, it can be expanded in elementary symmetric functions [E.sub.m] = [E.sub.m]([Z.sub.1], ..., [Z.sub.n]) defined by

[[PI].sup.n.sub.i=1](1 + t[Z.sub.i]) = [summation over(n/r=0)] [t.sup.r] [E.sub.r]. (3.6)

Applying Ref. (10), Eqs. (A.5) and (A.12) to the right side of Eq. (3.5), we find

[R.sub.-a]([beta], ..., [beta]; [z.sub.1], ..., [z.sub.n])

= [A.sup.-a][summation over([infinity]/N=0)] (a)N/[(n[beta]).sub.N] [T.sub.N]([beta], ..., [beta]; [Z.sub.1], ..., [Z.sub.n]), (3.7)

= [A.sup.-a][summation over([infinity]/N=0)] (a)N/[(n[beta]).sub.N][summation over([(-1).sup.N+M][([beta]).sub.M])] [E.sup.[m.sub.1].sub.1]...[E.sup.[m.sub.n].sub.n]/[m.sub.1]!...[m.sub .n]!, (3.8)

where [(a).sub.N] is Pochhammer's shifted factorial M = [summation over (n/i=1)] [m.sub.i], and the inner sum (representing [T.sub.N]) extends over all nonnegative integers [m.sub.1], ..., [m.sub.n], such that [m.sub.1] + 2[m.sub.2] + ... + [nm.sub.n] = N. Reference (10), Eq. (A.6) provides the recurrence relation

[NT.sub.N] + [summation over(n/r=2)] [(-1).sup.r](r[beta] + N - r)[E.sub.r][T.sub.N-r] = 0, (3.9)

where [T.sub.0] = 1 and [T.sub.N-r] = 0 if r > N, whence [T.sub.1] = 0. The term with r = 1 is missing because Eqs. (3.3) and (3.4) imply [E.sub.1] = 0, which greatly simplifies [T.sub.N].

Let \Z\= [max.sub.i] \[Z.sub.i]\ and [lambda] = max{\a\, 1}. The series [Eq. (3.7)] converges absolutely if [beta] > 0 and \Z\ < 1, and the truncation error [r.sub.k] resulting from neglect of terms of degree N [greater than or equal to] K can be shown to satisfy

\[r.sub.k]\ [less than or equal to] [(\a\).sub.k][\Z\.sup.k]/K![(1 - \Z\).sup.[lambda]] (3.10)

Each application of a duplication theorem decreases by a factor of four the differences between the variables, therefore the difference A - [z.sub.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 = [beta] = 1/2 and n = 3, the series [Eq. (3.8)] with [E.sub.1] = 0 can be rearranged as

[R.sub.F]([z.sub.1],[z.sub.2],[z.sub.3]) = [A.sup.-1/2][summation over([infinity]/r=0)] [summation over([infinity]/s=0)] [(-1).sup.r][(1/2).sub.r+s]/4r + 6s + 1 [E.sup.r.sub.2][E.sup.s.sub.3]/r!s!. (3.11)

This double series is convenient for obtaining numerical coefficients but requires \[E.sub.2]\ + \[E.sub.3]\ < 1 for absolute convergence. Keeping together all terms of the same degree in the Z's, as in Eq. (3.8), we find

[R.sub.F]([z.sub.1], [z.sub.2], [z.sub.3]) = [A.sup.-1/2](1 - 1/10 [E.sub.2] + 1/14 [E.sub.3] + 1/24 [E.sup.2.sub.2]

- 3/44 [E.sub.2][E.sub.3] - 5/208 [E.sup.3.sub.2] + 3/104 [E.sup.2.sub.3] + 1/16 [E.sup.2.sub.2][E.sub.3] + [r.sub.8]), (3.12)

where

\[r.sub.8]\ < 0.2[\Z\.sup.8]/1 - \Z\, \Z\ = [max.sub.i] \A - [Z.sub.i]/A < 1,

A = ([z.sub.1] + [z.sub.2] + [z.sub.3]/3. (3.13)

One more duplication would (ultimately) have decreased \[r.sub.8]by a factor of [4.sup.8] = 65 536.

For convenient reference we restate the corresponding result for [R.sub.J] and [R.sub.D] [see Eq. (3.1) and the discussions preceding it]:

[R.sub.-3/2](1/2, ..., 1/2; [z.zub.1], ..., [z.sub.5])

= [A.sup.-3/2] (1 - 3/14 [E.sub.2] + 1/6 [E.sub.3] + 9/88 [E.sup.2.sub.2] - 3/22 [E.sub.4]

- 9/52 [E.sub.2][E.sub.3] + 3/26 [E.sub.5] - 1/6 [E.sup.3.sub.2] + 3/40 [E.sup.2.sub.3] + 3/20 [E.sub.2][E.sub4]

+ 45/272 [E.sup.2.sub.2][E.sub.3] - 9/68 [E.sub.3][E.sub.4] - 9/68 [E.sub.2][E.sub.5] + [r.sub.8]), (3.14)

where

\[r.sub.8]\ < 3.4\[\Z\.sup.8]/[(1 - \Z\).sup.3/2]), \Z\ = [max.sub.i] \A - [z.sub.i]/A\ < 1,

A = ([z.sub.1] + ... + [z.sub.5])/5. (3.15)

The duplication theorems for these two functions are more complicated than that for [R.sub.F]; see Ref. (10) and Ref. (9), Appendix.

Acknowledgment

This manuscript was authored by a contractor of the U.S. Government (Ames Laboratory, Department of Energy), under contract No. W-7405-ENG-82. Accordingly, the U.S. Government retains a nonexclusive, royalty-free license to publish or reproduce the published form of this contribution, or allow others to do so for U.S. Government purposes. This work was supported also, in part, by the National Science Foundation under Agreement No. 9980036. Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the author and do not necessarily reflect the views of the National Science Foundation.

Accepted: August 6, 2002

4. References

(1.) B. C. Carlson, Toward symbolic integration of elliptic integrals, J. Symb. Comput. 28, 739-753 (1999).

(2.) I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series, and Products (Sixth Ed.), Academic Press, San Diego (2000).

(3.) M. Abramowitz and I. A. Stegun (eds.), Handbook of Mathematical Functions, U.S. Government Printing Office, Washington, D.C. (1964) Sec. 17.8 (Examples 16 and 18). Reprinted by Dover, New York (1965).

(4.) R. Bulirsch, An extension of the Bartky-transformation to incomplete elliptic integrals of the third kind, Numer. Math. 13, 266-284 (1969) Sec. 1.2.

(5.) B. C. Carlson, Quadratic transformations of Appell functions, SIAM J. Math. Anal. 7, 291-304 (1976).

(6.) B. C. Carlson, Special Functions of Applied Mathematics, Academic Press, New York (1977).

(7.) D. G. Zill and B. C. Carlson, Symmetric elliptic integrals of the third kind, Math. Comp. 24, 199-214 (1970).

(8.) B. C. Carlson, Numerical computation of real or complex elliptic integrals, Numer. Algorithms 10, 13-26 (1995).

(9.) B. C. Carlson and J. FitzSimons, Reduction theorems for elliptic integrands with the square root of two quadratic factors, J. Comput. Appl. Math. 118, 71-85 (2000).

(10.) B. C. Carlson, Computing elliptic integrals by duplication, Numer. Math. 33, 1-16 (1979).

About the author: B. C. Carlson was formerly a professor of physics and is currently a professor emeritus of mathematics at Iowa State University and an associate at the Ames Laboratory of the U.S. Department of Energy.