Gauss-Legendre algorithm

The method is based on the individual work of Carl Friedrich Gauss (1777 - 1855) and Adrien-Marie Legendre (1752-1833) combined with modern algorithms for multiplication and square roots. It repeatedly replaces two numbers by their arithmetic and geometric mean, in order to approximate their arithmetic-geometric mean.

The version presented below is also known as the Salamin-Brent algorithm; it was independently discovered in 1976 by Eugene Salamin[?] and Richard Brent[?]. It was used to compute the first 206,158,430,000 decimal digits of Pi on September 18 to 20, 1999, and the results were checked with Borwein's algorithm.

1. Initial value setting;

a = 1, b = 1 / √2, t = 1/4, p = 1

2. Repeat the following instructins until the difference of a and b is within the desired accuracy:

x = (a+b) / 2

y = √(a*b)

t = t - p * (a-x)²

a = x

b = y

p = 2 * p

3. Pi is approximated with a, b and t as:

Pi ≈ (a+b)² / (4*t)

The algorithm has second order convergent nature, which essentially means that the number of correct digits doubles with each step of the algorithm.