Gauss–Legendre algorithm

The Gauss–Legendre algorithm is an algorithm to compute the digits of π. It is notable for being rapidly convergent, with only 25 iterations producing 45 million correct digits of π. However, it has some drawbacks (for example, it is computer memory-intensive) and therefore all record-breaking calculations for many years have used other methods, almost always the Chudnovsky algorithm. For details, see Chronology of computation of π.

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 Gauss–Euler, Brent–Salamin (or Salamin–Brent) algorithm; it was independently discovered in 1975 by Richard Brent and Eugene Salamin. It was used in 1999 to compute the first 200 billion decimal digits of π, with results checked using Borwein's algorithm.

Algorithm

1. Initial value setting: Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle a_{0}=1\qquad b_{0}={\frac {1}{\sqrt {2}}}\qquad p_{0}=1\qquad t_{0}={\frac {1}{4}}.}
2. Repeat the following instructions until the difference between $a_{n+1}$ and Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle b_{n+1}} is within the desired accuracy: Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle {\begin{aligned}a_{n+1}&={\frac {a_{n}+b_{n}}{2}},\\\\b_{n+1}&={\sqrt {a_{n}b_{n}}},\\\\p_{n+1}&=2p_{n},\\\\t_{n+1}&=t_{n}-p_{n}(a_{n+1}-a_{n})^{2}.\\\end{aligned}}}
3. π is then approximated as: $$\pi \approx {\frac {(a_{n+1}+b_{n+1})^{2}}{4t_{n+1}}}.$$

The first five iterations give (approximations given up to and including the first incorrect digit):

Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle 3.140\dots }

Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle 3.14159264\dots }

$3.1415926535897932382\dots$

Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle 3.14159265358979323846264338327950288419711\dots }

$3.141592653589793238462643383279502884197169399375105820974944592307816406286208998625\dots$

The algorithm has quadratic convergence, which essentially means that the number of correct digits doubles with each iteration of the algorithm.

Mathematical background

Limits of the arithmetic–geometric mean

The arithmetic–geometric mean of two numbers, a₀ and b₀, is found by calculating the limit of the sequences

Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle {\begin{aligned}a_{n+1}&={\frac {a_{n}+b_{n}}{2}},\\[6pt]b_{n+1}&={\sqrt {a_{n}b_{n}}},\end{aligned}}}

which both converge to the same limit.
If $a_{0}=1$ and Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle b_{0}=\cos \varphi } then the limit is Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\textstyle {\pi \over 2K(\sin \varphi )}} where $K(k)$ is the complete elliptic integral of the first kind

$K(k)=\int _{0}^{\pi /2}{\frac {d\theta }{\sqrt {1-k^{2}\sin ^{2}\theta }}}.$

If $c_{0}=\sin \varphi$, $c_{i+1}=a_{i}-a_{i+1}$, then

$\sum _{i=0}^{\infty }2^{i-1}c_{i}^{2}=1-{E(\sin \varphi ) \over K(\sin \varphi )}$

where $E(k)$ is the complete elliptic integral of the second kind:

Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle E(k)=\int _{0}^{\pi /2}{\sqrt {1-k^{2}\sin ^{2}\theta }}\;d\theta }

Gauss knew of these two results.

Legendre’s identity

Legendre proved the following identity:

Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle K(\cos \theta )E(\sin \theta )+K(\sin \theta )E(\cos \theta )-K(\cos \theta )K(\sin \theta )={\pi \over 2},}

for all $\theta$.

Elementary proof with integral calculus

The Gauss-Legendre algorithm can be proven to give results converging to $\pi$ using only integral calculus. This is done here and here.

See also

• Numerical approximations of π