Elliptic Perimeters

Here is a brief reminder of the anatomy of an ellipse.

• An ellipse is the set of points $𝑃$ whose distances ${𝑟}_1$ and ${𝑟}_2$ to two fixed points, the foci, add up to a constant.

• The longest chord is the major axis, and half of its length is denoted by $𝑎$.
• The chord through the centre at right angles to it is the minor axis, and half of its length is denoted by $𝑏$.

• The distance from center to either foci is denoted by $𝑐$.

  Why is the constant $2𝑎$?

  Choose $𝑃$ at the right-hand end of the major axis; $𝑎-𝑐$ and $𝑎+𝑐$ are the distances to the foci, which add to $2𝑎$.
  Where does ${𝑎}^2={𝑏}^2+{𝑐}^2$ come from?

  Choose $𝑃$ at the top, at $\left(0,𝑏\right)$. The distances to either focus are equal by symmetry and still add to $2𝑎$, hence each is $𝑎$. It’s hopefully clear that ${𝑎}^2={𝑏}^2+{𝑐}^2$ due to Pythagoras!
• A circle is the case $𝑐=0$, so its eccentricity $𝑒=\frac{𝑐}{𝑎}$ is $0$; otherwise $0<𝑒<1$, and as $𝑒\to 1$ the ellipse flattens towards its major axis.
• We say a function / expression is elementary if it can be expressed as the sum, product, quotient or composition of the trigonometric, inverse-trigonometric, exponential, logarithmic, or n-th root functions.
• An algebraic function is a function $𝑦=𝑓\left(𝑥\right)$ that satisfies a finite polynomial in $𝑥$ and $𝑦$. For example, $𝑦={𝑥}^2$, $𝑦=\frac{\sqrt{𝑥}}{𝑥+3}$, $𝑦=𝑥\sqrt[3]{{𝑥}^2+4}$, etc.
• Liouville’s theorem tells us that if a function has an elementary integral, the integrand (the expression you’re integrating) must be expressed as the derivative of an algebraic function, plus the sum of a bunch of terms in the form $\frac{{𝑢}^{\prime }\left(𝑥\right)}{𝑢\left(𝑥\right)}$ (called a logarithmic derivative), where each $𝑢$ is an algebraic function.
• A circle’s area, $𝐴=𝜋{𝑟}^2$, and perimeter, $𝑃=2𝜋𝑟$, are elementary.

• An ellipse’s area, $𝐴=𝜋𝑎𝑏$ (note the symmetry with a circle, where $𝑟=𝑎=𝑏$), is elementary. However, its perimeter is not elementary. It is given by $𝐿=4𝑎𝐸\left(𝑒\right)$ where $𝐸$ is the “complete elliptic integral of the second kind”.

  Why is the perimeter’s integral not elementary?

  The integral

  Let’s parametrise to get $𝑥=𝑎sin𝜃$ and $𝑦=𝑏cos𝜃$. Using the arc-length formula on quadrant one:

  $${\int }_0^{\frac{𝜋}{2}}\sqrt{{𝑎}^2{cos}^2𝜃+{𝑏}^2{sin}^2𝜃}d𝜃$$

  Since

  $${𝑎}^2{cos}^2𝜃+{𝑏}^2{sin}^2𝜃={𝑎}^2\left(1-{𝑒}^2{sin}^2𝜃\right)$$

  we have the following for the total perimeter:

  $$𝐿=4𝑎{\int }_0^{\frac{𝜋}{2}}\sqrt{1-{𝑒}^2{sin}^2𝜃}d𝜃$$

  And as we will soon show, this integral has no elementary antiderivative, so we give it a name, $𝐸$. The same was done for $\text{erf}\left(𝑥\right)=\frac{2}{\sqrt{𝜋}}{\int }_0^{𝑥}{𝑒}^{-{𝑡}^2}d𝑡$. So we have:

  $$𝐿=4𝑎𝐸\left(𝑒\right)$$

  A brief introduction to complex analysis

  This proof uses some complex analysis. In a nutshell: consider a function like $\frac{7}{{𝑥}^2}+\frac{4}{𝑥}$. We call this the Laurent series expansion form; it is the sum of a bunch of terms in the form ${𝑎}_{𝑛}{\left(𝑧-{𝑧}_0\right)}^{𝑛}$. This function shoots off to $\infty$ as $𝑥⟶0$. At $𝑥=0$, we have a pole. We have different ways of classifying poles - one such way is via numbers called residues. The example above has residue 4, and we calculate it as the coefficient of the ${𝑥}^{-1}$ term. Another measure is the order, which is simply the highest negative power. The above example has order 2. There are some useful theorems / definitions concerning these:

  • Differentiation makes poles one order stronger. This is somewhat intuitive if we consider an example like $\frac{𝑑}{d𝑥}\frac{1}{𝑥}=-\frac{1}{{𝑥}^2}$ - the order increased from 1 to 2.
  • A holomorphic function is one that has a complex derivative at every point in a region. In layperson’s terms, when you zoom in close enough near any point, they can be expressed as a convergent power series (an infinite polynomial). Note that their differentiability means they can’t have poles.
  • Derivatives of the functions we’re considering (called meromorphic functions; these are holomorphic everywhere in a region except at their poles) always have residue 0. This should hopefully be clear: there’s no way for the derivative of a function containing finitely many powers of $𝑥$ to have a $\frac{1}{𝑥}$ term.

  • The residues of logarithmic derivatives are always integers. An example:

    $$𝑢\left(𝑥\right)={𝑥}^3,\frac{{𝑢}^{\prime }\left(𝑥\right)}{𝑢\left(𝑥\right)}=\frac{3{𝑥}^2}{{𝑥}^3}=\frac{3}{𝑥}$$

    which has residue $3$.

  That’s not quite the whole story, and there is a slight non-triviality. We are considering the integral of a differential $𝑓\left(𝑥\right)d𝑥$. When we change coordinates as we often do in calculus, both function and differential must be transformed too. This can cause poles to disappear. Consider for example $𝜔=\frac{1}{\sqrt{𝑥}}d𝑥$, and substituting $𝑥={𝑡}^2\text{so that}d𝑥=2𝑡d𝑡$, so we obtain $𝜔=2d𝑡$ which trivially has no pole.

  Rewriting the integral in a more useful form

  Now let’s make the substitution $𝑥=sin\left(𝜃\right)$, so that $d𝜃=\frac{d𝑥}{\sqrt{1-{𝑥}^2}}$:

  $${\int }_0^1\frac{\sqrt{1-{𝑒}^2{𝑥}^2}}{\sqrt{1-{𝑥}^2}}d𝑥$$

  We want to show that it’s not possible to represent it as a Liouville decomposition. First we multiply by the numerator over itself to get:

  $${\int }_0^1\frac{1-{𝑒}^2{𝑥}^2}{\sqrt{\left(1-{𝑥}^2\right)\left(1-{𝑒}^2{𝑥}^2\right)}}d𝑥$$

  Let ${𝑦}^2=\left(1-{𝑥}^2\right)\left(1-{𝑒}^2{𝑥}^2\right)$, so our differential is

  $$𝜔=\frac{𝑦}{1-{𝑥}^2}d𝑥$$

  Finding the poles of our differential

  We first want to find all the poles of this differential.

  At $𝑥=\pm 1$, we have a pole. However, we need to look at the entire differential, not just its coefficient. So, near $𝑥=1$, choose $𝑥=1-{𝑡}^2$, and hence $d𝑥=-2𝑡d𝑡$. Then:

  $$1-{𝑥}^2={𝑡}^2\left(2-{𝑡}^2\right)$$

  and

  $$𝑦=\pm 𝑡\sqrt{\left(2-{𝑡}^2\right)\left(1-{𝑒}^2{\left(1-{𝑡}^2\right)}^2\right)}$$

  Substituting everything into our differential gives

  $$𝜔=\pm 2\sqrt{\frac{1-{𝑒}^2{\left(1-{𝑡}^2\right)}^2}{2-{𝑡}^2}}d𝑡$$

  Notice that the factors of $𝑡$ have cancelled. Setting $𝑡=0$ is totally fine so at $𝑥=1$ we have no pole. Same logic applies for the other value, $𝑥=-1$.

  At $𝑥=\pm \frac{1}{𝑒}$, the differential is also regular. To see why, let ${𝑥}_0$ be $\frac{1}{𝑒}$ for example. Since ${𝑦}^2$ has a simple zero at ${𝑥}_0$, we can write:

  $${𝑦}^2=𝐶\left(𝑥-{𝑥}_0\right)+\text{a bunch of higher order terms}$$

  where $𝐶$ is a nonzero constant. We can do this legitimately because we’re near a zero, and a polynomial zoomed in looks like a line; in our case, $𝐶\left(𝑥-{𝑥}_0\right)$.

  Rearranging and differentiating gives

  $$d𝑥=\left(2\frac{𝑦}{𝐶}+\dots \right)d𝑦$$

  Note that $1-{𝑥}_0^2$ is nonzero; this is because if ${𝑥}_0=\pm \frac{1}{𝑒}$, and $0<𝑒<1$, then $\frac{1}{{𝑒}^2}>1$. Our differential contains a factor of $𝑦$ from its numerator, and another from $𝜔=\frac{𝑦}{1-{𝑥}^2}d𝑥$. So when we reconstruct $𝜔$ we’ll get something like:

  $$𝜔=\left(2\frac{{𝑦}^2}{𝐶\left(1-{𝑥}^2\right)}+\dots \right)d𝑦$$

  And as $𝑦⟶0$, this of course tends to 0 instead of becoming a pole. So I hope you’re convinced there are no poles at any finite point.

  But we’re not quite done! We still have to consider infinity. Let’s switch to $𝑧=\frac{1}{𝑥}$, so that instead of considering $𝑥⟶\infty$, we can consider 0 instead. To find the behaviour of our differential, first calculate its coefficient:

  $$\begin{array}{ccc}{\left(\frac{𝑦}{1-{𝑥}^2}\right)}^2& =\frac{{𝑦}^2}{{\left(1-{𝑥}^2\right)}^2}& =\frac{1-{𝑒}^2{𝑥}^2}{1-{𝑥}^2}\\ {\left(\frac{𝑦}{1-{𝑥}^2}\right)}^2& =\frac{{𝑧}^2-{𝑒}^2}{{𝑧}^2-1}& \\ \frac{𝑦}{1-{𝑥}^2}& =\pm \sqrt{\frac{{𝑧}^2-{𝑒}^2}{{𝑧}^2-1}}& \end{array}$$

  We can expand the expression inside the square root using the geometric series expansion:

  $$\frac{{𝑧}^2-{𝑒}^2}{{𝑧}^2-1}={𝑒}^2+\left({𝑒}^2-1\right){𝑧}^2+\left({𝑒}^2-1\right){𝑧}^4+\dots$$

  Taking the square root gives

  $$\frac{𝑦}{1-{𝑥}^2}=\pm \left(𝑒+\frac{{𝑒}^2-1}{2𝑒}{𝑧}^2+\dots \right)$$

  Only even powers of $𝑧$ appear! And since:

  $$d𝑥=\frac{-d𝑧}{{𝑧}^2}$$

  then:

  $$𝜔=\pm \left(\frac{𝑒}{{𝑧}^2}+\frac{{𝑒}^2-1}{2𝑒}+\dots \right)d𝑧$$

  There is a pole of order 2, but notice how there’s no term involving $\frac{1}{𝑧}$. Therefore, its residue is zero, so we have now shown that every residue of $𝜔$ is zero.

  So far so good. Now the proof gets maybe a little bit more complicated…

  What if an elementary antiderivative existed?

  FTSOC assume that our integral has an elementary primitive. By Liouville’s theorem, its differential must satisfy the decomposition described earlier:

  $$𝜔=d𝐴+\sum _{𝑖}{𝑐}_{𝑖}\frac{d{𝑢}_{𝑖}}{{𝑢}_{𝑖}}$$

  Recall that derivatives always have residue zero, whereas logarithmic derivatives have integer residues.

  Let

  $$𝜂=\sum _{𝑖}{𝑐}_{𝑖}\frac{d{𝑢}_{𝑖}}{{𝑢}_{𝑖}}$$

  be the sum of the logarithmic terms. Since both $𝜔$ and $d𝐴$ have residue zero everywhere, $𝜂$ must also have residue zero everywhere.

  Each function ${𝑢}_{𝑖}$ can only have finitely many zeros and poles. Let these points be ${𝑃}_1,\dots ,{𝑃}_{𝑚}$. At each point ${𝑃}_{𝑗}$, the residue of $𝜂$ is

  $$\sum _{𝑖}{𝑐}_{𝑖}{𝑛}_{𝑖𝑗}=0$$

  where ${𝑛}_{𝑖𝑗}$ is the order of the zero or pole of ${𝑢}_{𝑖}$ at ${𝑃}_{𝑗}$. These numbers are integers.

  Why do the logarithmic terms disappear?

  We now have to use a little bit of linear algebra for the sake of rigor. It’d be all too easy if we could just conclude all ${𝑐}_{𝑖}=0$, as it’s completely possible for some of them to cancel each other out.

  The equations above form a system of linear equations with integer coefficients. Consider something like $2{𝑐}_1-4{𝑐}_2=0$ which has $\left(2,1\right)$ as a solution. But note that any scalar multiple of this solution is also valid, and this generalises for bigger linear systems too. So its solutions can be written as linear combinations of vectors with rational values.

  More rigorously, there are rational vectors ${𝑚}_1,\dots ,{𝑚}_{𝑟}$ and constants ${𝜆}_1,\dots ,{𝜆}_{𝑟}$ such that the vector of coefficients $\left({𝑐}_1,\dots ,{𝑐}_{𝑛}\right)$ is a linear combination of these rational vectors. For each rational vector, we can multiply by a suitable integer to make all its entries integers, by virtue of the definition of a rational number. Call the resulting integer vectors ${𝑚}_{𝑘}=\left({𝑚}_{{𝑘}_1},\dots ,{𝑚}_{{𝑘}_{𝑛}}\right)$.

  The idea from here is to use the fact that if we multiply two functions together, their orders add, and if we raise a function to a power, its order is multiplied.

  Each of these integer vectors gives us a product

  $${𝑈}_{𝑘}=\prod _{𝑖}{𝑢}_{𝑖}^{{𝑚}_{{𝑘}_{𝑖}}}$$

  whose zeros and poles cancel at every point. Thus ${𝑈}_{𝑘}$ has no zeros or poles anywhere on our curve. But a nonconstant algebraic function on our curve must have zeros and poles.

  Why?

  On the compact version (it doesn’t have any missing end points or places where it stretches to infinity) of our curve, suppose we have a holomorphic function $𝑓$. Since the curve is compact, $\left|𝑓\right|$ must have a maximum somewhere. We want to show that a holomorphic function can only have a maximum of its absolute value if it is constant.

  Near any point ${𝑧}_0$, a holomorphic function can be written as a power series. If it isn’t locally constant, let ${𝑎}_{𝑛}{ℎ}^{𝑛}$ be the first nonconstant term. If $𝑓\left({𝑧}_0\right)$ isn’t zero, we can choose the direction of $ℎ$ so that this term points in the same direction as $𝑓\left({𝑧}_0\right)$. For sufficiently small $ℎ$, the higher order terms are smaller, so the modulus increases. This contradicts having a maximum. If $𝑓\left({𝑧}_0\right)=0$, a maximum of modulus zero means that $𝑓$ is zero in a neighbourhood anyway. Thus, a maximum of the modulus forces the function to be locally constant, and the Identity Theorem tells us it is constant everywhere on our connected curve.

  Therefore, if a meromorphic function has no poles, it is holomorphic everywhere and must be constant. The same argument applies to its reciprocal: if a function has no zeroes, its reciprocal has no poles, and must be constant. So a nonconstant meromorphic function must have both zeroes and poles. In fact, the number of poles and the number of zeroes is equal, counting multiplicity!

  Therefore, using the contrapositive, each ${𝑈}_{𝑘}$ is constant, meaning

  $$\frac{d{𝑈}_{𝑘}}{{𝑈}_{𝑘}}=\sum _{𝑖}{𝑚}_{{𝑘}_{𝑖}}\frac{d{𝑢}_{𝑖}}{{𝑢}_{𝑖}}=0$$

  Since the original logarithmic sum is a linear combination of these expressions, it must also be zero. Hence, all the logarithmic terms are gone, leaving

  $$𝜔=d𝐴$$

  Characterising A

  Since $d𝐴=𝜔$, and $𝜔$ has poles of order 2 at infinity and nowhere else, $𝐴$ can only have simple poles at infinity and no poles anywhere else. This is because of the property we established earlier, that differentiation increases pole order.

  Now let’s find all algebraic functions with this property. Since ${𝑦}^2$ is a polynomial in $𝑥$, we can express any function on our curve, which are functions of $𝑥,𝑦$ that involve:

  $$𝐴=𝑝\left(𝑥\right)+𝑞\left(𝑥\right)𝑦$$

  where $𝑝$ and $𝑞$ are functions of $𝑥$.

  Because $𝐴$ has no poles at finite points, neither $𝑝$ nor $𝑞$ can have a finite pole. To see why, consider the two points above any ordinary finite value of $𝑥$: one has $𝑦$ and the other has $-𝑦$ by symmetry. Adding and subtracting the two values of $𝐴$ shows that both $𝑝$ and $𝑞$ must be finite.

  At the branch points, the same conclusion follows by using $𝑥-{𝑥}_0={𝑡}^2$: the terms involving $𝑝$ contain even powers of $𝑡$, while those involving $𝑞𝑦$ contain odd powers, so their poles cannot cancel.

  Therefore, $𝑝$ and $𝑞$ must be polynomials.

  Why does this lead to a contradiction?

  Now consider what happens at infinity. Since

  $${𝑦}^2=\left(1-{𝑥}^2\right)\left(1-{𝑒}^2{𝑥}^2\right)$$

  the leading term of $𝑦$ is either $+𝑒{𝑥}^2$ or $-𝑒{𝑥}^2$, depending on which point at infinity we approach.

  If $𝑞$ is nonzero, then $𝑞\left(𝑥\right)𝑦$ has a pole of order at least 2 at infinity. The polynomial $𝑝\left(𝑥\right)$ cannot cancel this at both points at infinity, because the leading term of $𝑦$ has opposite signs at the two points. This contradicts the fact that $𝐴$ can have at most simple poles at infinity. Therefore, $𝑞=0$, leaving

  $$𝐴=𝑝\left(𝑥\right)$$

  Since $𝐴$ can have at most simple poles at infinity, $𝑝$ can have degree no more than 1. Hence:

  $$𝐴=𝛼𝑥+𝛽$$

  Differentiating gives

  $$d𝐴=𝛼d𝑥$$

  But this would mean

  $$\begin{array}{cc}\frac{𝑦}{1-{𝑥}^2}d𝑥& =𝛼d𝑥\\ 𝑦& =𝛼\left(1-{𝑥}^2\right)\end{array}$$

  Squaring and cancelling:

  $$\begin{array}{cc}\left(1-{𝑥}^2\right)\left(1-{𝑒}^2{𝑥}^2\right)& ={𝛼}^2{\left(1-{𝑥}^2\right)}^2\\ 1-{𝑒}^2{𝑥}^2& ={𝛼}^2\left(1-{𝑥}^2\right)\end{array}$$

  Comparing coefficients gives

  $${𝛼}^2=1\text{and}{𝛼}^2={𝑒}^2$$

  Hence ${𝑒}^2=1$, which is impossible for an ellipse with $0<𝑒<1$ $↯$. Therefore, the elliptic integral has no elementary antiderivative.

• Although it has no elementary form, we can approximate it numerically.

What are we doing?

We fix the semi-major axis at $𝑎=1$ and vary the semi-minor axis $𝑏$ between $0.0025$ and $1$, for a total of 400 shapes. Of course, f64 and single-threaded. For each of 400 shapes we compute the perimeter in several different ways and compare each result against a reference good to  32 digits.

The power mean

If I told you right now to guess the perimeter of an ellipse with a semimajor of 5 units and a semiminor of 1, I suspect a good number of you might think to average the two to get 3 and try considering the perimeter of a circle with radius 3. It’s not a terrible strategy, and I think it is worth considering. Though, there are lots of averages, so we’ll consider the most flexible which is the power mean, as we can optimise it a lot as you’ll soon see.

The arithmetic mean formula is at $𝑝=1$, the root mean square at $𝑝=2$, and the geometric mean as $𝑝\to 0$. They all don’t do too well though, and in fact the theoretical best $𝑝$ via a bit of algebra turns out to be $\frac{ln\left(2\right)}{ln\left(\frac{𝜋}{2}\right)}$ which is just over 1.53. The usual choice for this is $𝑝=3/2$ - it’s close enough, but I was curious how much $𝑝$ really matters.

The practical minimum is at $𝑝=1.527$, where the worst error is $2.3\times {10}^{\text{−3}}$, 4.5x better than at $3/2$. All things considered, given $𝑝=3/2$ differs by less than one order of magnitude, and dramatically simplifies our calculation, it is definitely worth using. With $𝑠=\left(𝑎\sqrt{𝑎}+𝑏\sqrt{𝑏}\right)/2$,

so it is two square roots and a cube root. The cube root is, as you might expect, the expensive part: with the library root the call costs 17.34 ns, more than the whole AGM. Some bit-hacking and Newton-Raphson a la Quake:

const MAGIC: u64 = 682 << 52;

fn cbrt(x: f64) -> f64 {
    let mut y = f64::from_bits(x.to_bits() / 3 + MAGIC);
    y = y * (2.0 / 3.0) + x / (3.0 * y * y);
    y
}

Dividing the floating-point bit representation by three gives a surprisingly good initial estimate of the cube root. This works because the exponent encodes the logarithmic scale of the number, so dividing the representation by three approximately divides its logarithm by three. An iteration of Newton-Raphson afterwards brings our error from 6% to just 0.3%, and altogether brings the call from 17.34 ns to just 3.75 ns.

In my research I learnt of another iterative method, called Halley’s method. Basically, the main differences are that Halley uses cubic convergence (NR is quadratic) and Halley also needs the second derivative (NR needs only the first).

As you can see, Halley is pretty much a free win, $1.21\times {10}^{\text{−4}}$ against the $3.15\times {10}^{\text{−3}}$ of Newton for 5% more time, since both are one division and a few multiplies. The perimeter column looks like it disagrees, but it is not measuring the root: at $𝑝=3/2$ the power mean itself is only good to about ${10}^{-2}$ at the flat end and it comes out low there, so any root better than that lands on the same number. Newton only looks better because one step from the seed always lands above the true root, by about 0.3% here, and that bias cancels part of the formula’s error, pulling the worst case to $7.28\times {10}^{\text{−3}}$ while Halley sits at the formula’s own $1.03\times {10}^{\text{−2}}$.

Ramanujan’s first

Ramanujan published this formula (and another one) in 1914, of course without derivation.

p = pi * (3*(a + b) - sqrt((3*a + b)*(a + 3*b)))

It needs one square root and a few muls. At $\frac{𝑏}{𝑎}=0.5$ the error is $2.8\times {10}^{-6}$, at $0.1$ it is $8.4\times {10}^{-4}$, and its worst case is $3.95\times {10}^{\text{−3}}$ at the flatter end of the graph near $𝑒=0$. That’s still 2.5 orders better than the geometric mean for the same square root, and 0.5x the error of the power mean at 0.6x the cost. A call takes 2.17 ns.

Ramanujan’s second

The second formula moves the square root inside the correction and adds a division.

t = (a - b) / (a + b)
h = t*t
p = pi * (a + b) * (1 + 3*h / (10 + sqrt(4 - 3*h)))

It’s not too bad! Up to 9 digits for a division and a square root, though it is 7-8 digits over most of the range and degrades at the flat end, where the correction saturates as $ℎ\to 1$: the worst error is $3.47\times {10}^{\text{−4}}$.

A call takes 4.48 ns, 2x the first formula, for 3 or 4 more digits in the middle of the range. It is a good place to stop if a few digits are enough; though Ramanujan left one more correction in his notes.

Ramanujan’s second with his correction

The last formula in the notebooks is the second one with a term added inside the bracket, $3{ℎ}^5/2^{17}$. This is a correction found by considering the difference between the exact perimeter series and his approximation.

t = (a - b) / (a + b)
h = t*t
h5 = h*h*h*h*h
p = pi * (a + b) * (1 + 3*h / (10 + sqrt(4 - 3*h)) + 3*h5 / 131072)

It costs five more multiplications but no new square root. At 2:1 it takes the error from $4.56\times {10}^{\text{−10}}$ to $7.89\times {10}^{\text{−11}}$, which is 6x better and 10 digits, but at the flat end it is more or less the same: $3.3\times {10}^{\text{−4}}$ vs $3.47\times {10}^{\text{−4}}$.

So regrettably the gain isn’t too helpful: 4.79 ns against 4.48 ns is  6% more time for  5% less error. But I would be remiss to not at least consider it.

The obvious series

So far everything we’ve looked at is a basically an elaborate guess. The alternative is to expand the integral itself. We can expand for eccentricity as follows:

m = 1 - (b/a)*(b/a)
g = 1
mp = 1
s = 0
for n in 1..=max_terms:
    g *= (2*n - 1) / (2*n)
    mp *= m
    term = g*g * mp / (2*n - 1)
    s += term
    if term < 1e-18: break
p = 2 * pi * a * (1 - s)

As expected, though the expansion is exact in the limit, it’s impractical at the flat end. Its terms go like ${𝑒}^{2𝑛}/{𝑛}^2$, so the sum converges quickly only while $𝑒$ is small, which fails as the ellipse flattens. At $\frac{𝑏}{𝑎}=0.5$ it needs quite a ridiculous 106 terms to reach the last bit; at $0.1$, 2393, and at $0.01$,  400,000; and at the flat end of the sweep,  6 million. The bench gives it 5000 and it is still wrong by only the 5th decimal there. So all in all, in its current state, it’s basically useless. But the principle itself can be improved!

Brute force

Thousands of terms is not particularly appealing. Another option is numerical integration.

end = pi / 2
h = end / 64
total = f(0) + f(end)
for i in 1..64:
    y = f(i * h)
    total += (i odd ? 4 : 2) * y
p = 4 * total * h / 3

With 64 panels, Simpson evaluates a sine, a cosine and a square root 65 times, so it is again quite expensive. It is the most accurate method here where the ellipse is nice and round, $1.83\times {10}^{\text{−16}}$ at 2:1. Better than any of the formulas above, but the error climbs unforunately quite rapidly at a sharp corner on the graph at about $𝑒=0.08$: by $𝑏/𝑎=0.01$ it is $1.48\times {10}^{\text{−5}}$.

Panels are somewhat crude. We can do better using Gauss-Legendre, which to put it simply, places them where they do the most good. A node costs about what a panel costs (11 ns of sine, cosine and square root). Error-wise, 4 nodes against 4 panels is $8.82\times {10}^{\text{−4}}$ against $4.26\times {10}^{\text{−3}}$, 16 against 16 is $4.78\times {10}^{\text{−6}}$ against $2.63\times {10}^{\text{−4}}$, and 32 against 32 is $2.35\times {10}^{\text{−7}}$ against $6.57\times {10}^{\text{−5}}$.

A final method is Carlson’s symmetric integrals. Instead of sampling the integrand, they write the perimeter as a symmetric integral of three arguments, $𝐿=8𝑎{𝑅}_{𝐺}\left(0,1-{𝑘}^2,1\right)$, and it uses the same duplication loop logic as AGM. It is the most accurate method here at $6.66\times {10}^{\text{−16}}$ for 18.13 ns, about the same as the AGM.

A better variable

Let’s go back to the converging series. The series in ${𝑒}^2$ converges slowly not because the integral is hard or something but because $𝑒$ is a poor variable. Expanding in $ℎ={\left(\left(𝑎-𝑏\right)/\left(𝑎+𝑏\right)\right)}^2$ instead gives a series with positive coefficients:

t = (a - b) / (a + b)
h = t*t
c = 1
hp = 1
s = 1
for n in 1..=max_terms:
    c *= (n - 1.5)*(n - 1.5) / (n*n)
    hp *= h
    term = c * hp
    s += term
    if term < 1e-18 * s: break
p = pi * (a + b) * s

$ℎ$ is bounded by 1 and is usually far smaller than $𝑒$: at $\frac{𝑏}{𝑎}=0.5$ it is only $\frac{1}{9}$. Table 2 compares the two. At $\frac{𝑏}{𝑎}=0.5$ the $ℎ$ series needs 15 terms where the series in ${𝑒}^2$ needs 106, and at $0.1$ it is 66 against 2393. It is the same integral and the same loop with one substitution, and at the flat end of the sweep the two differ by 10 orders of magnitude on the same term budget.

Though it is exact, it is still not cheap; 3-4x the AGM for the same digits.

Remez

Another method we can use is the Chebyshev minimax, in particular the Remez exchange algorithm; this basically tries to fit a polynomial to minimise our maximum error. For what it’s worth, I think it is a very beautiful way to approximate a whole host of functions. If you’re interested, check out the Sollya project for more.

We fit the ratio $𝑃/\left(𝜋\left(𝑎+𝑏\right)\right)$ that Ramanujan’s formulas use, not $𝑃$ itself, and the criterion is relative error in that ratio: the algorithm moves $\left({𝑃}_{\text{fit}}-𝑃\right)/𝑃$ to the same size at every extremum, alternating in sign. Since $𝜋\left(𝑎+𝑏\right)$ is exact for a given shape, that is relative error in $𝑃$ too, which is the number the tables report, so nothing needs weighting by $1/𝑃$. $\left[5/5\right]$ in $ℎ$:

num = p0 + h*(p1 + h*(p2 + h*(p3 + h*(p4 + h*p5))))
den = 1 + h*(q1 + h*(q2 + h*(q3 + h*(q4 + h*q5))))
p = pi * (a + b) * num / den

The fit is one-off, over 801 evenly spaced points in whichever variable is being fitted, $ℎ$, $𝑠$ or $𝑢$, across the whole sweep ($𝑏/𝑎$ from 0.0025 to 1).

With 10 multiply-adds and a division, from $\left[3/3\right]$ to $\left[5/5\right]$ we have $1\times {10}^{\text{−6}}$, $9.49\times {10}^{\text{−8}}$, $1.01\times {10}^{\text{−8}}$. We can improve this by making the fit be the square root. In fact, Ramanujan’s second formula is a rational function of $\sqrt{4-3ℎ}$ rather than of $ℎ$, which is why it holds up at the flat end. In

the singularity ($ℎ=1$, i.e. the degenerate ellipse $𝑏=0$) moves from 1% outside the interval to 10% outside it. We see way more bang for our buck this way.

You’ll notice I also tested cube root for this, so that $𝑢={\left(1-ℎ\right)}^{1/3}={𝑠}^{2/3}$, which reduces our error again. $\left[5/5\right]$ reaches $1.49\times {10}^{\text{−13}}$ against $4.28\times {10}^{\text{−12}}$. Unfortunately, it’s unsurprisingly slower than sqrt: 15.88 ns against 7.83 ns. The AGM at 16.95 ns is cheaper and 3 orders more accurate.

The arithmetic geometric mean

I mentioned this one a couple times in passing in the previous methods. I suspect that this is the most widely used approximation. It’s pretty simple, and very accurate, but not too fast. It essentially replaces a pair of numbers by their arithmetic and geometric means and repeats; it rapidly converges too.

where ${𝑐}_0=𝑘$ and ${𝑐}_{𝑛}=\frac{{𝑎}_{𝑛-1}-{𝑏}_{𝑛-1}}{2}$. Then $𝐿=4𝑎𝐸\left(𝑒\right)$, as before.

r = b / a
k = sqrt(1 - r*r)
x = 1
y = r
c = k
s = k*k / 2
w = 0.5
for turn in 1..=8:
    xn = (x + y) / 2
    yn = sqrt(x*y)
    c = (x - y) / 2
    x = xn
    y = yn
    w *= 2
    s += w * c*c
    if w * c*c < 1e-18 * s: break
p = 4 * a * (pi / (2*x)) * (1 - s)

Incredibly, each turn roughly doubles the number of correct digits. 5 turns suffice at 2:1 and 6 for anything flatter in the sweep, and the worst error anywhere is $1.11\times {10}^{\text{−15}}$, which we can attribute to f64 rounding.

16.95 ns, 4x Ramanujan’s second but still less than either series or Simpson.

Which one to use

If speed is the priority, Ramanujan’s first is the best of the cheapest ones; use the second if you realise you want a bit more accuracy. If accuracy is the priority, AGM and Carlson (which is what Boost C++ uses) will both serve you well. Otherwise, for a best of both worlds, I would suggest Remez (sqrt).

Sources

• B. Conrad, Impossibility theorems for elementary integration, 2005. Example 4.3 gives the general idea of the argument used above.
• S. Ramanujan, Modular equations and approximations to $𝜋$, Quarterly Journal of Mathematics 45 (1914), 350-372.
• J. M. Borwein and P. B. Borwein, Pi and the AGM, Wiley, 1987. The arithmetic geometric mean and the elliptic integrals.
• E. Ya. Remez, Sur la determination des polynomes d’approximation de degre donnee, Communications de la Societe Mathematique de Kharkov 10 (1934), 41-63.
• B. C. Carlson, Numerical computation of real or complex elliptic integrals, Numerical Algorithms 10 (1995), 13-26. The symmetric integrals and the duplication behind them.
• NIST Digital Library of Mathematical Functions, chapter 19, dlmf.nist.gov/19. The relations between the symmetric and Legendre forms.
• Boost.Math, special_functions, for the ${𝑅}_{𝐺}$ branch and the stored polynomials described above.