Gauss–Legendre algorithm: Difference between revisions

Revision as of 18:32, 10 January 2011

Latest revision as of 09:52, 15 June 2025

(77 intermediate revisions by 48 users not shown)

Line 1:

{{Short description|Quickly converging computation of π}}

The '''Gauss–Legendre algorithm''' is an [[algorithm]] to compute the digits of [[Pi|π{{pi}}]]. It is notable for being rapidly convergent, with only 25 iterations producing 45 million correct digits of π {{pi}}. However, theit drawbackhas issome thatdrawbacks (for example, it is memory[[Random-access_memory|computer memory]]-intensive) and ittherefore isall thereforerecord-breaking sometimescalculations notfor many years have used overother methods, almost always the [[Machin-likeChudnovsky algorithm]]. For details, see [[chronology of computation of π|Chronology of computation of formulas{{pi}}]].

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 root]]s. It repeatedly replaces two numbers by their [[arithmetic mean|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''';<ref>[[Richard Brent (scientist)|Brent, Richard]], ''Old and New Algorithms for pi'', Letters to the Editor, Notices of the AMS 60(1), p. 7</ref> it was independently discovered in 1975 by [[Richard Brent (scientist)|Richard Brent]] and [[Eugene Salamin (mathematician)|Eugene Salamin]]. 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]].

== Algorithm ==

# Initial value setting: <math display="block">a_0 = 1\qquad b_0 = \frac{1}{\sqrt{2}}\qquad p_0 = 1\qquad t_0 = \frac{1}{4}\qquad p_0 = 1\!.</math>

2.# Repeat the following instructions until the difference ofbetween <math>a_n\!a_{n+1}</math> and <math>b_n\!b_{n+1}</math> is within the desired accuracy: <math display="block"> \begin{align}

1. Initial value setting:

:<math> \begin{align} a_{n+1} & = \frac{a_n + b_n}{2}, \\

p_{n+1} & = 2p_n. \\

:<math>a_0 = 1\qquad b_0 = \frac{1}{\sqrt{2}}\qquad t_0 = \frac{1}{4}\qquad p_0 = 1\!</math>

b_{n+1} & = \sqrt{a_n b_n}, \\

\\

2. Repeat the following instructions until the difference of <math>a_n\!</math> and <math>b_n\!</math> is within the desired accuracy:

p_{n+1} & = 2p_n, \\

\\

:<math> \begin{align} a_{n+1} & = \frac{a_n + b_n}{2}, \\

b_t_{n+1} & = \sqrtt_n - p_n(a_{a_n b_nn+1}-a_{n},)^2. \\

t_{n+1} & = t_n - p_n(a_n - a_{n+1})^2, \\

p_{n+1} & = 2p_n.

\end{align}

</math>

# {{pi}} is then approximated as: <math display="block">\pi \approx \frac{(a_{n+1}+b_{n+1})^2}{4t_{n+1}}.</math>

3. π is then approximated as:

:<math>\pi \approx \frac{(a_n+b_n)^2}{4t_n}.\!</math>

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

:<math>3.140\dots\!</math>

:<math>3.14159264\dots\!</math>

:<math>3.1415926535897932382\dots\!</math>

:<math>3.14159265358979323846264338327950288419711\dots</math>

:<math>3.141592653589793238462643383279502884197169399375105820974944592307816406286208998625\dots</math>

The algorithm has second-order[[quadratic convergent natureconvergence]], which essentially means that the number of correct digits doubles with each step[[iteration]] of the algorithm.

== Mathematical background ==

=== Limits of the arithmetic-geometricarithmetic–geometric mean ===

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

:<math>\begin{align} a_{n+1} & = \frac{a_n+b_n}{2}, \\[6pt]

b_{n+1} & = \sqrt{a_n b_n},

\end{align}

</math>

which both converge to the same limit.<br />

If <math>a_0=1\!</math> and <math>b_0=\cos\varphi\!</math> then the limit is <math display="inline">{\pi \over 2K(\sin\varphi)}\!</math> where <math>K(k)\!</math> is the [[Elliptic integral#Complete elliptic integral of the first kind|complete elliptic integral of the first kind]]

:<math>K(k) = \int_0^{\pi/2} \frac{d\theta}{\sqrt{1-k^2 \sin^2\theta}}.\!</math>

If <math>c_0 = \sin\varphi\!</math>, <math>c_{i+1} = a_i - a_{i+1}\!</math>., then

:<math>\sum_{i=0}^\infty 2^{i-1} c_i^2 = 1 - {E(\sin\varphi)\over K(\sin\varphi)}\!</math>

where <math>E(k)\!</math> is the [[Elliptic integral#Complete elliptic integral of the second kind|complete elliptic integral of the second kind]]:

:<math>E(k) = \int_0^{\pi/2}\sqrt {1-k^2 \sin^2\theta}\,; d\theta.\!</math>

Gauss knew of boththese of thesetwo results.<ref name="brent">{{Citation

| last=Brent

| first=Richard

| author-link=Richard Brent (scientist)

| publication-date=

| date=

| year=1975

| title=Multiple-precision zero-finding methods and the complexity of elementary function evaluation

Line 79 ⟶ 75:

| oclc=

| accessdate=8 September 2007

| archive-date=23 July 2008

}}</ref>

| archive-url=https://web.archive.org/web/20080723170157/http://wwwmaths.anu.edu.au/~brent/pub/pub028.html

<ref name="salamin1">[[Eugene Salamin (mathematician)|Salamin, Eugene]]. ''Computation of pi'', Charles Stark Draper Laboratory ISS memo 74–19, 30 January, 1974, Cambridge, Massachusetts</ref>

| url-status=dead

}}</ref>

<ref name="salamin1">[[Eugene Salamin (mathematician)|Salamin, Eugene]]., ''Computation of pi'', Charles Stark Draper Laboratory ISS memo 74–19, 30 January, 1974, Cambridge, Massachusetts</ref>

<ref name="salamin2">{{Citation

| last=Salamin

Line 86 ⟶ 85:

| author-link=Eugene Salamin (mathematician)

| publication-date=

| date=1976

| year=1976

| title=Computation of pi Using Arithmetic-GeometricArithmetic–Geometric Mean

| periodical=Mathematics of Computation

| series=

Line 100 ⟶ 98:

| pages=565–570

| url=

| issn=0025--5718

| doi=10.2307/2005327

| jstor=2005327

| oclc=

| accessdate=

Line 107 ⟶ 106:

=== Legendre’s identity ===

Legendre proved the following identity:

:<math display="block">K(\sincos \varphitheta) E(\sin \theta ) + K(\sin \theta ) E(\sincos \varphitheta) - K(\sincos \varphitheta) K(\sin \theta) = {1\pi \over 2}\pi\!,</math><ref name="brent" />

for all <math>\theta</math>.<ref name="brent" />

=== Elementary proof with integral calculus ===

For <math>\varphi\!</math> and <math>\theta\!</math> such that <math>\varphi+\theta={1 \over 2}\pi\!</math> Legendre proved the identity:

:<math>K(\sin \varphi) E(\sin \theta ) + K(\sin \theta ) E(\sin \varphi) - K(\sin \varphi) K(\sin \theta) = {1 \over 2}\pi\!</math><ref name="brent" />

The Gauss-Legendre algorithm can be proven to give results converging to <math>\pi</math> using only integral calculus. This is done here<ref>{{citation|title=Recent Calculations of π: The Gauss-Salamin Algorithm|last1=Lord|first1=Nick|doi=10.2307/3619132|year=1992|journal=The Mathematical Gazette|volume=76|issue=476|pages=231–242|jstor=3619132|s2cid=125865215 }}</ref> and here.<ref>{{citation|title=Easy Proof of Three Recursive π-Algorithms|last1=Milla|first1=Lorenz|arxiv=1907.04110|year=2019}}</ref>

=== Gauss–Legendre method ===

The values <math>\varphi=\theta={\pi\over 4}\!</math> can be substituted into Legendre’s identity and the approximations to K, E can be found by terms in the sequences for the arithmetic geometric mean with <math>a_0=1\!</math> and <math>b_0=\sin{\pi \over 4}=\frac{1}{\sqrt{2}}\!</math>.<ref name="brent" />

== See also ==

* [[Numerical approximations of π|Numerical approximations of {{pi}}]]

== References ==

{{reflist}}

{{DEFAULTSORT:Gauss-Legendre algorithm}}

[[Category:Pi algorithms]]

[[es:Algoritmo de Gauss-Legendre]]

[[it:Algoritmo di Gauss-Legendre]]

[[he:אלגוריתם גאוס-לז'נדר]]

[[nl:Algoritme van Gauss-Legendre]]

[[ja:ガウス＝ルジャンドルのアルゴリズム]]

[[tr:Gauss-Legendre Algoritması]]

[[zh:高斯-勒让德算法]]