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, the drawback is that it is computer memory-intensive and therefore sometimes Machin-like formulas are used instead.

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;[1] it was independently discovered in 1975 by Richard Brent and Eugene Salamin. It was used to compute the first 206,158,430,000 decimal digits of π on September 18 to 20, 1999, and the results were checked with Borwein's algorithm.

Algorithm

1. Initial value setting:

${\displaystyle a_{0}=1\qquad b_{0}={\frac {1}{\sqrt {2}}}\qquad t_{0}={\frac {1}{4}}\qquad p_{0}=1.}$

2. Repeat the following instructions until the difference of ${\displaystyle a_{n}}$ and ${\displaystyle b_{n}}$ is within the desired accuracy:

${\displaystyle {\begin{aligned}a_{n+1}&={\frac {a_{n}+b_{n}}{2}},\\\\b_{n+1}&={\sqrt {a_{n}b_{n}}},\\\\t_{n+1}&=t_{n}-p_{n}(a_{n}-a_{n+1})^{2},\\\\p_{n+1}&=2p_{n}.\\\end{aligned}}}$

3. π is then approximated as:

${\displaystyle \pi \approx {\frac {(a_{n+1}+b_{n+1})^{2}}{4t_{n+1}}}.}$

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

${\displaystyle 3.140\dots }$

${\displaystyle 3.14159264\dots }$

${\displaystyle 3.1415926535897932382\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

${\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 ${\displaystyle a_{0}=1}$ and ${\displaystyle b_{0}=\cos \varphi }$ then the limit is ${\textstyle {\pi \over 2K(\sin \varphi )}}$ where ${\displaystyle K(k)}$ is the complete elliptic integral of the first kind

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

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

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

where ${\displaystyle E(k)}$ is the complete elliptic integral of the second kind:

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

Gauss knew of both of these results.[2] [3] [4]

Legendre’s identity

For ${\displaystyle \varphi }$ and ${\displaystyle \theta }$ such that ${\textstyle \varphi +\theta ={1 \over 2}\pi }$ Legendre proved the identity:

${\displaystyle K(\sin \varphi )E(\sin \theta )+K(\sin \theta )E(\sin \varphi )-K(\sin \varphi )K(\sin \theta )={1 \over 2}\pi .}$[2]

Equivalently,

${\displaystyle \forall \varphi :K(\sin \varphi )[E(\cos \varphi )-K(\cos \varphi )]+K(\cos \varphi )E(\sin \varphi )={\frac {\pi }{2}}}$

Elementary proof with integral calculus

The Gauss-Legendre algorithm can be proven without elliptic modular functions. This is done here[5] and here[6] using only integral calculus.

See also

• Numerical approximations of π

References

1. Brent, Richard, Old and New Algorithms for pi, Letters to the Editor, Notices of the AMS 60(1), p. 7
2. Brent, Richard (1975), Traub, J F (ed.), "Multiple-precision zero-finding methods and the complexity of elementary function evaluation", Analytic Computational Complexity, New York: Academic Press, pp. 151–176, retrieved 8 September 2007
3. Salamin, Eugene, Computation of pi, Charles Stark Draper Laboratory ISS memo 74–19, 30 January 1974, Cambridge, Massachusetts
4. Salamin, Eugene (1976), "Computation of pi Using Arithmetic–Geometric Mean", Mathematics of Computation, 30 (135), pp. 565–570, doi:10.2307/2005327, ISSN 0025-5718
5. Lord, Nick (1992), "Recent Calculations of π: The Gauss-Salamin Algorithm", The Mathematical Gazette, 76 (476): 231–242, doi:10.2307/3619132, JSTOR 3619132
6. Milla, Lorenz (2019), Easy Proof of Three Recursive π-Algorithms, arXiv:1907.04110