Calculating pi with AGM and mpmath

This post gives an algorithm based on the arithmetic-geometric mean that rapidly converges to pi. I’ll use it to illustrate multiple precision arithmetic using Python’s mpmath module.

Given two non-negative numbers a and b, their arithmetic mean is (a + b)/2 and their geometric mean is √(ab). Suppose you start with two non-negative numbers and take their arithmetic mean and geometric mean. Then take the arithmetic and geometric mean of the results. Keep doing this repeatedly, and the result converges to what is called the arithmetic-geometric mean (AGM).

There is a formula for calculating pi based on the AGM that goes back to Gauss.

pi = frac{ 4M^2left(1, frac{1}{sqrt{2}}right)}{1 - sum_{n=1}^infty 2^{n+1} (a_n^2 - b_n^2)}

Here M(a, b) is the AGM of a and b, and an and bn are the values after n steps of the iteration.

The process defining the AGM converges quickly, and so the formula is practical for computing pi. I found it in a paper from 1988, and at that time the formula had been used to compute pi to over 29 million decimal places. In the code below, we compute pi to 997 decimal places in only 10 iterations.

from mpmath import mp, sqrt, mpf, pi, log10, fabs

# carry out calculations to 1000 decimal places
decimal_places = 1000
mp.dps = decimal_places
epsilon = 1/mpf(10**decimal_places)

# set a = 1 and b = 1/sqrt(2) as multi-precision numbers
a = mpf(1)
b = 1/sqrt(mpf(2))

diff = a - b
series = mpf(0)

n = 0
while diff > epsilon:
    n += 1
    arith = (a + b)/2
    geom = sqrt(a*b)
    a, b = arith, geom
    series += 2**(n+1) * (a*a - b*b)
    diff = a - b

# a and b have converged to the AGM
my_pi = 4*a*a/(1 - series)

error = fabs(pi - my_pi)
decimal_places = int(-log10(error))

print "Number of steps used: %d" % n
print "Number of correct decimal places: %d" % decimal_places

The code can be used to calculate pi out further. The number of correct decimal places roughly doubles with each iteration. For example, computing pi to 10,000 places takes only 3 more iterations.

Related posts:

10 thoughts on “Calculating pi with AGM and mpmath

  1. Reference: “Gauss, Landen, Ramanujan, the Arithmetic-Geometric Mean, Ellipses, pi, and the Ladies Diary” by Gert Almevist and Bruce Berndt. Available here.

  2. sherifffruitfly

    what about the computation’s stability characteristics? in particular, if the summation is close to 1, or a_i close to b_i ?

  3. Steven: All calculations are carried out to the same precision, in this case 1000 decimal places. I suppose you’re asking how the square root of 2 accuracy propagates throughout the rest of the calculation since it is computed first. I don’t have a direct answer to that, but the final result of carrying all calculations out to 1000 places is that 3 decimal places are lost.

    I don’t know just where those decimal places are lost. It would be interesting to look at this computation in more detail and see which steps contribute the most to the error.

    It seems that about the same number of decimal places are lost as the precision varies. For example, when I reran the calculation with 20,000 decimal precision, I still lose 3 decimal place, i.e. the result is correct to 19,997 places. I suppose the precision loss depends on the number of iterations, not the precision. Using 20,000 decimal places only required 14 iterations, only 4 more than using 1,000 decimal places.

    When I go up to 100,000 decimal places, it takes 17 iterations, and I lose 5 decimal places of precision. This is more evidence that the number of decimal places lost increases very slowly as the working precision increases.

  4. On a hunch that the expression (a*a - b*b) was contributing to loss of precision, I eventually discovered that the following version of the loop gives full precision:

    while diff > epsilon:
    series += 2**n * diff**2
    n += 1
    arith = (a + b)/2
    geom = sqrt(a*b)
    a, b = arith, geom
    diff = a - b

    This comes from substituting (a + b)/2 for a and sqrt(a*b) for b in the expression (a*a - b*b), so that each term in the series is based on the previous values of a and b.

  5. Loop is done while diff (=a-b) > epsilon.
    Final error is pi – 4*a*a/(1 – series), not diff (=a-b).

  6. UltimateQuantumGuy

    Having come to this website by chance after attempting a calculation of my own using the Borwein quartic convergence AGM algorithm, I thought it might be of interest for you to know that in my implementation, I found a striking resemblance between your error in the decimals and my own – having noted almost exactly the same results with my version.

Leave a Reply

Your email address will not be published. Required fields are marked *