The AGM on PCalc

I’ve been messing around recently with the arithmetic-geometric mean (AGM). Not that I care much about the AGM itself, but it’s a convenient way to calculate the elliptic integrals that come up in certain mechanics problems. If you follow Matt Parker on YouTube, you may have seen the AGM make a guest appearance in his video about calculating the area of a squircle from last fall (which also mentioned elliptic integrals).

If I were writing a script that needed elliptic integrals, I’d just use the ellipk function from the scipy.special library. But sometimes it’s nice to have a way to do things on your calculator. Hence the PCalc function below.

The procedure for calculating the AGM is iterative and usually converges in just a handful of steps. Here’s how we get the AGM of two numbers, which we’ll call aa and bb:

  1. Start by taking aa and bb as our initial pair of estimates, a0a_0 and b0b_0.
  2. Calculate their arithmetic mean and use it as one number in the next iteration:

    ai=ai1+bi12a_i = \frac{a_{i-1} + b_{i-1}}{2}
  3. Calculate their geometric mean and use it as the other number in the next iteration:

    bi=ai1bi1b_i = \sqrt{a_{i-1} b_{i-1}}
  4. Repeat Steps 2 and 3 until the two means converge.

As an example, here’s how we get the AGM of 1 and 2:

Step a b
0 1 2
1 1.50000 1.41421
2 1.45711 1.45648
3 1.45679 1.45679

As you can see, the arithmetic mean is always larger than the geometric mean. The MathWorld article linked above explains why that is and also why the two get closer together with each iteration. Some people like to define the problem with ai<bia_i < b_i, and so they prefer to reverse the formulas for aia_i and bib_i, but since both the arithmetic and geometric means use commutative operations, the order doesn’t matter.

This pairwise iteration works really well in PCalc, especially if, like me, you work with PCalc in RPN mode. Here’s the function I made that takes aia_i and bib_i on the stack and returns ai+1a_{i+1} and bi+1b_{i+1}, ready for the next iterative step. If aia_i and bib_i are equal, then it returns just the one converged number.

PCalc AGM Step

A few things to note:

To calculate the AGM, I put the two numbers on the stack and call the AGM Step function repeatedly until convergence. Fast and simple.

If you want the AGM Step function and are too lazy to key in the steps yourself, you can download it.