\documentclass[11pt]{article}
\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{epstopdf}
\usepackage{hyperref}
\DeclareGraphicsRule{.tif}{png}{.png}{`convert #1 `dirname #1`/`basename #1 .tif`.png}

\textwidth = 6.5 in
\textheight = 9 in
\oddsidemargin = 0.0 in
\evensidemargin = 0.0 in
\topmargin = 0.0 in
\headheight = 0.0 in
\headsep = 0.0 in
\parskip = 0.2in
\parindent = 0.0in

\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{definition}{Definition}
\newtheorem{lemma}{Lemma}

\title{AGM calculation of Pi}
\author{Burton Rosenberg}
\begin{document}
\maketitle

The fastest method for calculating $\pi$ was discovered by Gauss.
He found a connection between a process called the {\em arithmetic-geometric mean}
and what we now call {\em elliptic integrals}.
The two integrals of interest are,
\[
I(a,b) = \int_0^{\pi/2} \frac{d\phi}{\sqrt{a^2 \cos^2\phi+b^2\sin^2\phi}}
\]
and a very similar integral,
\[
L(a,b) = \int_0^{\pi/2} \frac{\cos^2 \phi\,d\phi}{\sqrt{a^2 \cos^2\phi+b^2\sin^2\phi}}
\]
In 1748 Euler had discovered the relation,
\[
L(\sqrt{2},1) I(\sqrt{2},1) = \pi/4
\]
Gauss showed that these integrals can be evaluated by successively replacing
$a$ and $b$ by their arithmetic and geometric means,
\begin{eqnarray*}
   a_{i+1} &=& (a_i+b_i)/2\\
   b_{i+1} &=& \sqrt{a_ib_i}
\end{eqnarray*}
where $a_0=a$ and $b_0=b$.
The sequence $a_i$ and
$b_i$ converge extremely quickly to a common value, called the {\em arithmetic-geometric
mean of $a$ and $b$}, denoted ${\rm AGM}(a,b)$.  
Specifically,  what Gauss discovered that if $a_i$ and $b_i$ follow 
from $a$ and $b$ as a result of one
or several steps of the AGM procedure, the integral is unchanged, $I(a,b)=I(a_i,b_i)$.
Hence, in the limit, the common value ${\rm AGM}(a,b)$ can be factored out and the integrand
collapses to $1/\sqrt{\cos^2\theta+\sin^2\theta}=1$,
\[
I(\sqrt{2},1)=\frac{\pi}{2\,{\rm AGM}(\sqrt{2},1)}.
\]
Furthermore, he showed,
\[
c_0^2 L(a,b)  = (c_0^2-S)I(a,b),
\]
where $
S=\sum_{k=0}^\infty 2^{k-1}c_k^2
$,  $c_0=a^2-b^2$ and $c_i=a_i^2-b_i^2$.

Putting these results together,
\[
\pi = \frac{{\rm AGM}(\sqrt{2},1)^2}{1-S}.
\]
The use of this formula is that the AGM converges extremely quickly, and the values $c_i$
converge towards zero extremely quickly. With just several terms we already have a very
good approximation of $\pi$.


We take heavily from $\pi$ unleashed, Jorg Arndt and Christoph Haenel.
and Nick Lord.

\begin{definition}
\[\
I(a,b)=\int_0^{\pi/2}\frac{ d\phi}{\sqrt{ a^2\cos^2 \phi + b^2 \sin^2 \phi}}
\]
\end{definition}

\begin{definition}
\[
L(a,b) = \int_0^{\pi/2} \frac{\cos^2\phi\,d\phi}{\sqrt{a^2\cos^2\phi+b^2\sin^2\phi}}
\]
\end{definition}

\begin{lemma}
$L(a,b)+L(b,a)=I(a,b)$.
\end{lemma}

 Substitute $\phi=\pi/2-\phi'$ into $L(b,a)$,
\begin{eqnarray*}
L(b,a) &=&
  \int_{\pi/2}^{0} \frac{\cos^2(\pi/2-\phi')\,d(\pi/2-\phi')}
                  {\sqrt{b^2\cos^2(\pi/2-\phi')+a^2\sin^2(\pi/2-\phi)}}\\
&=& \int_0^{\pi/2} \frac{\sin^2 \phi'\,d\phi'}
                  {\sqrt{b^2\sin^2\phi'+a^2\cos^2\phi}}\\
\end{eqnarray*}
Hence $L(a,b)+L(b,a)$ has $\cos^2\phi+\sin^2\phi=1$ in the numerator of 
the integrand and thus reduces to $I(a,b)$.

\begin{lemma}
\[
I(a,b)=\frac{1}{2}\int_{-\infty}^\infty \frac{dx}{\sqrt{b^2+x^2}\sqrt{a^2+x^2}}
\]
\end{lemma}

Using the identity,
\[
\sec \phi = \sqrt{1+\tan \phi}
\]
we have,
\begin{eqnarray*}
\frac{ d\phi}{\sqrt{ a^2\cos^2 \phi + b^2 \sin^2 \phi}} 
&=&  \frac{d\phi}{\cos\phi \sqrt{a^2+b^2\tan^2\phi}}\\
&=&  \frac{b\sec^2\phi \,d\phi}{b\sec\phi \sqrt{a^2+b^2\tan^2\phi}}\\
&=&  \frac{b\sec^2\phi \,d\phi}{\sqrt{b^2+b^2\tan^2\phi} \sqrt{a^2+b^2\tan^2\phi}}
\end{eqnarray*}
Making the substitution $x = b \tan \phi$, with the differential,
\[
   dx =  b \sec^2 \phi \; d \phi
\]
we have,
\[
I(a,b)=\frac{1}{2}\int_{-\infty}^\infty \frac{dx}{\sqrt{b^2+x^2}\sqrt{a^2+x^2}}
\]

\begin{lemma}
\[
L(a,b) = \frac{1}{2}\int_{-\infty}^\infty\frac{b^2\,dx}{(b^2+x^2)\sqrt{a^2+x^2}\sqrt{b^2+x^2}}
\]
\end{lemma}

Proceeding as with $I(a,b)$,
\begin{eqnarray*}
\frac{ \cos^2\phi\,d\phi}{\sqrt{ a^2\cos^2 \phi + b^2 \sin^2 \phi}} 
&=&  \frac{d\phi}{\sec^2\phi\,\cos\phi \sqrt{a^2+b^2\tan^2\phi}}\\
&=&  \frac{b^3\sec^2\phi \,d\phi}{(b^2+b^2\tan\phi)\sqrt{b^2+b^2\tan^2\phi} \sqrt{a^2+b^2\tan^2\phi}}
\end{eqnarray*}
substituting $x=b\tan\phi$,
\[
L(a,b) = \frac{1}{2}\int_{-\infty}^\infty\frac{b^2\,dx}{(b^2+x^2)\sqrt{a^2+x^2}\sqrt{b^2+x^2}}
\]

\begin{lemma}
$I(a_1,b_1)=I(a,b)$ where $a_1=(a+b)/2$ and $b_1=\sqrt{ab}$.
\end{lemma}

Substituting $x=(1/2)(x'-ab/x')$,
\begin{eqnarray*}
\sqrt{a_1^2+x^2} &=& \sqrt{(a+b)^2/4 + (1/4)(x'^2-2ab+(ab/x')^2)}\\
&=& (1/2x') \sqrt{x'^4+(a^2+b^2)x'^2+(ab)^2}\\
&=&1/(2x')\sqrt{x'^2+a^2}\sqrt{x'^2+b^2}\\
~\\
\sqrt{b_1^2+x^2} &=& \sqrt{ab+(1/4)(x'^2-2ab+(ab/x')^2)}\\
   &=&1/(2x')(x'^2+ab)
\end{eqnarray*}
and,
\[
dx=(1/2)d(x'-ab/x') = 1/(2x'^2)( x'^2 +ab)dx'
\]
So,
\[
I(a_1,b_1)=\frac{1}{2}\int_{-\infty}^{\infty}\frac{dx}{\sqrt{b_1^2+x^2}\sqrt{a_1^2+x^2}}
=\frac{1}{2} \int_{0}^{\infty}\frac{2\,dx'}{\sqrt{b^2+x'^2}\sqrt{a^2+x'^2}}=I(a,b).
\]

\begin{lemma}
\[
L(b_1,a_1) =
\frac{a+b}{a-b}(L(b,a)-L(a,b))
\]
\end{lemma}

Substitute  $x=(1/2)(x'-ab/x)'$ into $L(b_1,a_1)$.
Note,
\begin{eqnarray*}
\frac{a_1^2}{a_1^2+((1/2)(x'-ab/x'))^2} 
&=& \frac{(a+b)^2x'^2}{ (a+b)^2x'^2+x'^4-2abx'^2+(ab)^2}\\
&=&\frac{(a+b)^2x'^2}{(a^2+x'^2)(b^2+x'^2)}\\
&= &\frac{a+b}{a-b}\left(\frac{a^2}{a^2+x'^2}-\frac{b^2}{b^2+x'^2} \right)
\end{eqnarray*}
Combining this with the identities above,
\begin{eqnarray*}
L(b_1,a_1) &=&
\frac{1}{2}\int_{-\infty}^\infty\frac{a_1^2}{a_1^2+x^2}\frac{dx}{\sqrt{b_1^2+x^2}\sqrt{a_1^2+x^2}}\\
&=&
\frac{a+b}{a-b}\int_{0}^{\infty} \left(\frac{a^2}{a^2+x'^2}-\frac{b^2}{b^2+x'^2} \right)\frac{dx'}{\sqrt{a^2+x'^2}\sqrt{b^2+x'^2}}\\
&=&\frac{a+b}{a-b}(L(b,a)-L(a,b))
\end{eqnarray*}

\begin{lemma}
\[
c_0^2L(b,a) = \left( \sum_{i=0}^{\infty}2^{i-1}c_i^2\right)I(a,b)
\]
\end{lemma}

Since,
\[
a_1^2-b_1^2= (a-b)^2/4
\]
and defining $c_0=a^2-b^2$ and $c_1^2=a_1^2-b_1^2$,
we can write the above relation as,
\[
4c_1^2 L(b_1,a_1) = c_0^2(L(b,a)-L(a,b))
\]
Using $I(a,b)=L(a,b)+L(b,a)$, 
\[
c_0^2L(b,a)-2c_1^2L(b_1,a_1)=(1/2)c_0^2I(a,b).
\]
Repeating the AGM to get $a_i, b_i$ and $c_i$, we get further relations,
\[
2^i(c_i^2L(b_i,a_i)-2c_{i+1}^2L(b_{i+1},a_{i+1}))=2^i(1/2)c_i^2I(a,b)
\]
Summing over $i$ and letting the sum telescope,
\[
c_0^2L(b,a)-2^nc_n^2L(b_n,a_n)=\left( \sum_{i=0}^{n-1}2^{i-1}c_i^2\right)I(a,b)
\]

 \end{document} 