Algorithms for the arithmetic-geometric mean

With complex variables, it is convenient to work with the univariate function M(z)=agm(1,z). The general case is given by agm(a,b)=aM(1,b/a).

Functional equation

If the real part of z initially is not completely nonnegative, we apply the functional equation M(z)=(z+1)M(u)/2 where u=z/(z+1).

Note that u has nonnegative real part, absent rounding error. It is not a problem for correctness if rounding makes the interval contain negative points, as this just inflates the final result.

For the derivative, the functional equation becomes M(z)=[M(u)(z1)M(u)/((1+z)z)]/2.

AGM iteration

Once z is in the right half plane, we can apply the AGM iteration (2an+1=an+bn,bn+12=anbn) directly. The correct square root is given by ab, which is computed as ab,iab,iab,ab respectively if both a and b have positive real part, nonnegative imaginary part, nonpositive imaginary part, or otherwise.

The iteration should be terminated when an and bn are close enough. For positive real variables, we can simply take lower and upper bounds to get a correct enclosure at this point. For complex variables, it is shown in [Dup2006], p. 87 that, for z with nonnegative real part, |M(z)an||anbn|, giving a convenient error bound.

Rather than running the AGM iteration until an and bn agree to p bits, it is slightly more efficient to iterate until they agree to about p/10 bits and finish with a series expansion. With z=(ab)/(a+b), we have

agm(a,b)=(a+b)π4K(z2),

valid at least when |z|<1 and a,b have nonnegative real part, and

π4K(z2)=1218z25128z411512z646932768z8+

where the tail is bounded by k=10|z|k/64.

First derivative

Assuming that z is exact and that |arg(z)|3π/4, we compute (M(z),M(z)) simultaneously using a finite difference.

The basic inequality we need is |M(z)|max(1,|z|), which is an immediate consequence of the AGM iteration.

By Cauchy’s integral formula, |M(k)(z)/k!|CDk where C=max(1,|z|+r) and D=1/r, for any 0<r<|z| (we choose r to be of the order |z|/4). Taylor expansion now gives

|M(z+h)M(z)hM(z)|CD2h1Dh|M(z+h)M(zh)2hM(z)|CD3h21Dh|M(z+h)+M(zh)2M(z)|CD2h21Dh

assuming that h is chosen so that it satisfies hD<1.

The forward finite difference would require two function evaluations at doubled precision. We use the central difference as it only requires 1.5 times the precision.

When z is not exact, we evaluate at the midpoint as above and bound the propagated error using derivatives. Again by Cauchy’s integral formula, we have

|M(z+ε)|max(1,|z|+|ε|+r)r|M(z+ε)|2max(1,|z|+|ε|+r)r2

assuming that the circle centered on z with radius |ε|+r does not cross the negative half axis. We choose r of order |z|/2 and verify that all assumptions hold.

Higher derivatives

The function W(z)=1/M(z) is D-finite. The coefficients of W(z+x)=k=0ckxk satisfy

2z(z21)c2=(3z21)c1+zc0,
(k+2)(k+3)z(z21)ck+3=(k+2)2(3z21)ck+2+(3k(k+3)+7)zck+1+(k+1)2ck

in general, and

(k+2)2ck+2=(3k(k+3)+7)ck+1+(k+1)2ck

when z=1.