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;^[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

Initial value setting: Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a_0 = 1\qquad b_0 = \frac{1}{\sqrt{2}}\qquad p_0 = 1\qquad t_0 = \frac{1}{4}.}
Repeat the following instructions until the difference between Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a_{n+1}} and $b_{n+1}$ is within the desired accuracy: Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \begin{align} 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{align} }
$π$ is then approximated as: Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\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):

3.140\dots

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 3.14159264\dots}

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 3.1415926535897932382\dots}

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 3.141592653589793238462643\dots}

3.14159265358979323846264338327\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 (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \begin{align} a_{n+1} & = \frac{a_n+b_n}{2}, \\[6pt] b_{n+1} & = \sqrt{a_n b_n}, \end{align} }

which both converge to the same limit.
If Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a_0=1} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle b_0=\cos\varphi} then the limit is Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\textstyle {\pi \over 2K(\sin\varphi)}} where Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle K(k)} is the complete elliptic integral of the first kind

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle K(k) = \int_0^{\pi/2} \frac{d\theta}{\sqrt{1-k^2 \sin^2\theta}}.}

If Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle c_0 = \sin\varphi} , $c_{i+1}=a_{i}-a_{i+1}$ , then

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sum_{i=0}^\infty 2^{i-1} c_i^2 = 1 - {E(\sin\varphi)\over K(\sin\varphi)}}

where Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle E(k)} is the complete elliptic integral of the second kind:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle E(k) = \int_0^{\pi/2}\sqrt {1-k^2 \sin^2\theta}\; d\theta}

Gauss knew of these two results.^[2] ^[3] ^[4]

Legendre’s identity

Legendre proved the following identity:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle K(\cos \theta) E(\sin \theta ) + K(\sin \theta ) E(\cos \theta) - K(\cos \theta) K(\sin \theta) = {\pi \over 2},}

for all Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \theta} .^[2]

Elementary proof with integral calculus

The Gauss-Legendre algorithm can be proven to give results converging to π using only integral calculus. This is done here^[5] and here.^[6]

References

^ Brent, Richard, Old and New Algorithms for pi, Letters to the Editor, Notices of the AMS 60(1), p. 7
^ ^2.0 ^2.1 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, archived from the original on 23 July 2008, retrieved 8 September 2007
^ Salamin, Eugene, Computation of pi, Charles Stark Draper Laboratory ISS memo 74–19, 30 January 1974, Cambridge, Massachusetts
^ Salamin, Eugene (1976), "Computation of pi Using Arithmetic–Geometric Mean", Mathematics of Computation, vol. 30, no. 135, pp. 565–570, doi:10.2307/2005327, ISSN 0025-5718, JSTOR 2005327
^ Lord, Nick (1992), "Recent Calculations of π: The Gauss-Salamin Algorithm", The Mathematical Gazette, 76 (476): 231–242, doi:10.2307/3619132, JSTOR 3619132, S2CID 125865215
^ Milla, Lorenz (2019), Easy Proof of Three Recursive π-Algorithms, arXiv:1907.04110

[1] Brent, Richard, Old and New Algorithms for pi, Letters to the Editor, Notices of the AMS 60(1), p. 7

[brent-2] 2.0 ^2.1 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, archived from the original on 23 July 2008, retrieved 8 September 2007

[salamin1-3] Salamin, Eugene, Computation of pi, Charles Stark Draper Laboratory ISS memo 74–19, 30 January 1974, Cambridge, Massachusetts

[salamin2-4] Salamin, Eugene (1976), "Computation of pi Using Arithmetic–Geometric Mean", Mathematics of Computation, vol. 30, no. 135, pp. 565–570, doi:10.2307/2005327, ISSN 0025-5718, JSTOR 2005327

[5] Lord, Nick (1992), "Recent Calculations of π: The Gauss-Salamin Algorithm", The Mathematical Gazette, 76 (476): 231–242, doi:10.2307/3619132, JSTOR 3619132, S2CID 125865215

[6] Milla, Lorenz (2019), Easy Proof of Three Recursive π-Algorithms, arXiv:1907.04110

[1]

[2]

[3]

[4]

[5]

[6]