The AGM on PCalc
February 26, 2022 at 1:16 PM by Dr. Drang
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 and :
- Start by taking and as our initial pair of estimates, and .
Calculate their arithmetic mean and use it as one number in the next iteration:
Calculate their geometric mean and use it as the other number in the next iteration:
- 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 , and so they prefer to reverse the formulas for and , 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 and on the stack and returns and , ready for the next iterative step. If and are equal, then it returns just the one converged number.
A few things to note:
- PCalc doesn’t have any looping commands, so there’s no way to have a single function call go through all the iterations necessary to converge. You might think—as I do whenever I start writing a PCalc function that would benefit from looping—that you could simulate looping by giving a negative number of steps to the Skip command, but negative numbers aren’t allowed.
- The check on equality in the first step is the sort of thing you’re always warned against doing with floating point numbers. And I’ve run into cases where PCalc never considers the two numbers on the stack to be equal. But since PCalc doesn’t loop, I’m in control of the iteration. If the two numbers look the same to me, it doesn’t matter what PCalc thinks—I’ll just stop running the function and drop the redundant value from the stack.
- The function uses Registers 0 and 1 for intermediate calculations and doesn’t clear them. This may look dangerous, but because I never write functions that assume starting values for the registers, it isn’t.
To calculate the AGM, I put the two numbers on the stack and call the function repeatedly until convergence. Fast and simple.
If you want the function and are too lazy to key in the steps yourself, you can download it.
