The circumference of an ellipse

Share on:

Note: This blog post is mainly computational, with a hint of proof-oriented mathematics here and there. For a more in-depth analysis, read the excellent article "Gauss, Landen, Ramanujan, the Arithmetic-Geometric Mean, Ellipses, pi, and the Ladies Diary" by Gert Akmkvist and Bruce Berndt, in The American Mathematical Monthly, vol 95 no. 7 (August-September 1988), pages 585-608, and happily made available online for free at

This article was the 1989 winner of the MAA Halmos-Ford award "For outstanding expository papers in The American Mathematical Monthly" and you should certainly read it.


It's one of the delights or annoyances of mathematics, which ever way you look at it, that there's no simple formula for the cicrcumference of an ellipse comparable to \(2\pi r\) for a circle.

Indeed, for an ellipse with semi axes \(a\) and \(b\), the circumference can be expressed as the integral

\[ 4\int_0^{\pi/2}\sqrt{a^2\cos^2\theta + b^2\sin^2\theta}\,d\theta \]

which is one of a class of integrals called elliptic integrals, and which cannot be expressed using algebraic or standard transcendental functions.

However, it turns out that there are ways of very quickly computing the circumference on an ellipse to any desired accuracy, using methods which originate with Carl Friedrich Gauss.

The arithmetic-geometric mean

Given \(a>b>0\), we can define two sequences by:

\[ a_0 = a,\qquad a_{k+1}=(a_k+b_k)/2 \]


\[ b_0=b,\qquad b_{k+1}=\sqrt{a_kb_k}. \]

Thus the \(a\) values are the arithmetic means of the previous pair; the \(b\) values the geometric mean of the pair.


\[ b<\sqrt{ab}<\frac{a+b}{2}<a \]

and since \(a_{k+1}<a_k\) and \(b_{k+1}>b_k\), it follows that the \(a\) values are decreasing and bounded below, and the \(b\) values are increasing and bounded above, so they both converge. Also, if \(c_k=\sqrt{a_k^2-b_k^2}\) then (see Almqvist \& Berndt pp 587-588):

\begin{align*} c_{k+1}&=\sqrt{a_{k+1}^2-b_{k+1}^2}\\ &=\sqrt{\frac{1}{4}(a_k+b_k)^2-a_kb_k}\\ &=\frac{1}{2}(a_k-b_k)\\ &=\frac{c_k^2}{4a_{k+1}}\\ &<\frac{c_k^2}{4M(a,b)}. \end{align*}

This shows that not only do the sequences converge to the same limit, but that the sequences converge quadratically; each iteration being double the precision of the previous.

The common limit is called the arithmetic-geometric mean of \(a\) and \(b\) and will be denoted here as \(M(a,b)\).

To give an indication of this speed, in Python:

 1  from mpmath import mp
 2  mp.dps = 50
 3  a,b = mp.mpf('3'), mp.mpf('2')
 4  for i in range(10):
 5      a,b = (a+b)/2, mp.sqrt(a*b)
 6      print('{0:52s} {1:52s}'.format(str(a),str(b)))
 82.5                                                  2.4494897427831780981972840747058913919659474806567
 92.4747448713915890490986420373529456959829737403283  2.4746160019198827700554874647235766528956885806854
102.4746804366557359095770647510382611744393311605069  2.474680435816873015671798992747485357556612143505
112.4746804362363044626244318718928732659979716520059  2.4746804362363044625888873352899297125054403531594
122.4746804362363044626066596035914014892517060025827  2.4746804362363044626066596035914014892516421855508
132.4746804362363044626066596035914014892516740940667  2.4746804362363044626066596035914014892516740940667
142.4746804362363044626066596035914014892516740940667  2.4746804362363044626066596035914014892516740940667
152.4746804362363044626066596035914014892516740940667  2.4746804362363044626066596035914014892516740940667
162.4746804362363044626066596035914014892516740940667  2.4746804362363044626066596035914014892516740940667
172.4746804362363044626066596035914014892516740940667  2.4746804362363044626066596035914014892516740940667

You see that we have reached the limit of precision in six steps. To better demonstrate the speed of convergence, use greater precision and display the difference \(a_k-b_k\):

 1  mp.dps = 2000
 2  a,b = mp.mpf('3'), mp.mpf('2')
 3  for i in range(10):
 4      a,b = (a+b)/2, mp.sqrt(a*b)
 5      mp.nprint(a-b,10)

You can see that the precision does indeed roughly double at each step.

Elliptic integrals and the AGM

This integral:

\[ K(k) = \int_0^{\pi/2}\frac{1}{\sqrt{1-k^2\sin^2\theta}}d\theta \]

is called a complete elliptic integral of the first kind. If

\[ k^2=1-\frac{b^2}{a^2} \]


\[ \int_0^{\pi/2}\frac{1}{\sqrt{1-k^2\sin^2\theta}} d\theta = a\int_0^{\pi/2} \frac{1}{\sqrt{a^2\cos^2\theta +b^2\sin^2\theta}} d\theta \]

and the integral on the right is denoted \(I(a,b)\). Gauss observed (and proved) that

\[ I(a,b)=\frac{\pi}{2M(a,b)}. \]

This is in fact equivalent to the assertion that

\[ K(x) = \frac{\pi}{2M(1-x,1+x)}. \]

Gauss started with assuming that \(M(1-x,1+x)\) was an even function, so could be expressed as a power series in even powers of \(x\). He then compared the power series representing \(K(x)\) with the integral computed by expanding the integrand as a power series using the binomial theorem and integrating term by term, to obtain

\[ \frac{1}{M(1-x,1+x)}=\frac{2}{\pi}K(x)=1+\left(\frac{1}{2}\right)^{2}x^2+ \left(\frac{1\cdot 3}{2\cdot 4}\right)^{2}x^4+ \left(\frac{1\cdot 3\cdot 5}{2\cdot 4\cdot 6}\right)^{2}x^6+\cdots \]

and this power series can be written as

\[ \sum_{n=0}^{\infty}\frac{(2n)!^2}{2^{4n}(n!)^4}x^{2n}. \]

Another proof, apparently originally due to Newman, uses first the substition \(t=a\tan\theta\) to rewrite the integral:

\[ I(a,b)=\int_0^{\pi/2}\frac{1}{\sqrt{a^2\cos^2\theta+b^2\sin^2\theta}} d\theta = \int^{\infty}_{-\infty}\frac{1}{\sqrt{(a^2+t^2)(b^2+t^2)}} dt \]

and make the substitution

\[ u = \frac{1}{2}\left(t-\frac{ab}{t}\right) \]

which produces (after "some care", according to the Borwein brothers in their book Pi and the AGM):

\[ I(a,b) = \int^{\infty}_{-\infty}\frac{1}{\sqrt{(\left(\frac{a+b}{2}\right)^2+t^2)(ab+t^2)}} dt =I(\frac{a+b}{2},\sqrt{ab}) = I(a_1,b_1). \]

Continuing this process we see that \(I(a,b) = I(a_k,b_k)\) for any \(k\), and taking the limit, that \(I(a,b)=I(M,M)\). Finally

\[ I(M,M) = \int_0^{\pi/2}\frac{1}{\sqrt{M^2\cos^2\theta+M^2\sin^2\theta}} d\theta = \int_0^{\pi/2}\frac{1}{M} d\theta = \frac{\pi}{2M}. \]

So we now know that the complete elliptic integral of the first kind is computable by means of the AGM. But what about the complete elliptic integral of the second kind?

Complete elliptic integrals of the second kind

These are defined as

\begin{align*} E(k)&=\int_0^{\pi/2}\sqrt{1-k^2\sin^2\theta} d\theta\\ &=\int_0^1\sqrt{\frac{1-k^2t^2}{1-t^2}}\,dt. \end{align*}

Note: An alternative convention is to write:

\begin{align*} E(m)&=\int_0^{\pi/2}\sqrt{1-m\sin^2\theta} d\theta\\ &=\int_0^1\sqrt{\frac{1-mt^2}{1-t^2}}\,dt. \end{align*}

This elliptic integral is the one we want, since the perimeter of an ellipse given in cartesion form as

\[ \frac{x^2}{a^2}+\frac{y^2}{b^2}=1 \]

is equal to

\[ 4aE\left(\sqrt{1-\frac{b^2}{a^2}}\right) \]

(using the first definition). This can be written as

\[ 4aE(e) \]

where \(e\) is the eccentricity of the ellipse.

The computation of elliptic integrals goes back as far as Euler, taking in Gauss and then Legendre along the way. Adrien-Marie Legendre (1752-1833) is better known now than he was in his lifetime; his contributions to elliptic integrals are now seen as fundamental, and in paving the way for the greater work of Abel and Jacobi.

In particular, Legendre showed that

\[ \frac{K(k)-E(k)}{K(k)}=\frac{1}{2}(c_0^2+2c_1^2+4c_2^2+\cdots + 2^nc_n^2+\cdots ) \]


\[ c_n^2=a_n^2-b_n^2. \]

The above equation can be rewritten to provide a fast-converging series for \(E(k)\):

\begin{align*} E(k)&=K(k)(1-\frac{1}{2}(c_0^2+2c_1^2+4c_2^2+\cdots + 2^nc_n^2+\cdots))\\ &=\frac{\pi}{2M(a,b)}(1-\frac{1}{2}(c_0^2+2c_1^2+4c_2^2+\cdots + 2^nc_n^2+\cdots)). \end{align*}

In the sequences for the AGM, let \(M(a,b)\) be approximated by \(a_n\). This produces a sequence that converges very quickly to \(E(k)\):

\begin{align*} e_0&=\frac{\pi}{2a}(1-\frac{1}{2}c_0^2)\\ e_1&=\frac{\pi}{2a_1}(1-\frac{1}{2}(c_0^2+2c_1^2))\\ e_2&=\frac{\pi}{2a_2}(1-\frac{1}{2}(c_0^2+2c_1^2+4c_2^2)) \end{align*}

This can easily be managed recursively, given \(a\) and \(b\), by setting

\[ a_0=a,\qquad,b_0=b,\qquad p_0=1,\qquad s_0=a^2-b^2 \]

and then iterating by

\begin{align*} a_{k+1}&=(a_k+b_k)/2\\ b_{k+1}&=\sqrt{a_kb_k}\\ p_{k+1}&=2p_k\\ s_{k+1}&=s_k+p_k(a_{k+1}^2-b_{k+1}^2) \end{align*}

Then the values

\[ e_k = \frac{2\pi}{a_k}(a_0^2-\frac{1}{2}s_k) \]

approach the perimeter of the ellipse.

We can demonstrate this in Python for \(a,b=3,2\), again using the multiprecision library mpmath:

 1  mp.dps = 100
 2  a = mp.mpf(3)
 3  b = mp.mpf(2)
 4  p = 1
 5  s = a**2-b**2
 6  a0 = a
 7  e = 2*mp.pi/a*(a0*a0-1/2*s)
 8  for i in range(11):
 9      a, b, p = (a+b)/2, np.sqrt(a*b), p*2
10      s += p*(a**2-b**2)
11      print(2*mp.pi/a*(a0*a0-1/2*s))

We see that we have reached the limits of this precision very quickly. And as above, we can look instead at the differences between succesive approximations:

 2a = mp.mpf(3)
 3b = mp.mpf(2)
 4p = 1
 5s = a**2-b**2
 6a0 = a
 7e = 2*mp.pi/a*(a0*a0-1/2*s)
 8for i in range(11):
 9    a, b, p = (a+b)/2, np.sqrt(a*b), p*2
10    s += p*(a**2-b**2)
11    e1 = 2*mp.pi/a*(a0*a0-1/2*s)
12    mp.nprint(e1-e,10)
13    e = e1

The convergence is seen to be quadratic.

Using Excel

It might seem utterly regressive to use Excel - or indeed, any spreadsheet - for computations such as this, but in fact Excel can be used very easily for recursive computations, as long as you're happy to use only IEEE double precision.

In Cells A1-E2, enter:

1 3 2 1 =A1^2-B1^2
2 =(A1+B1)/2 =SQRT(A1*B1) =2*C1 =D1+C2*(A2^2-B2^2) =2*PI()/A2*($A$1^2-1/2*D2)

and then copy A2-E2 down a few rows. Because of overflow errors, if you show 14 decimal places in column E, they'll never quite settle down. But in fact you reach the limits of double precision by row 5.