(41C) and (42S) Arithmetic-Geometric Mean

[Eddie W. Shore]

Arithmetic-Geometric Mean

The program AGM calculates the arithmetic-geometric mean of two positive integers x and y. As the graphic above suggests, an iterative process is used to find the AGM, computing both the arithmetic mean and geometric mean until the two means converge.

a0 = x
g0 = y

Repeat:
Arithmetic Mean: a1 = (a0 + g0)/2
Geometric Mean: g1 = √(a0 * g0)
Transfer new to old: a0 = a1, g0 = g1
Until |a1 - g1| < tolerance

You can set the tolerance as low as you want. The programs presented on this blog set tolerance at 10^(-10) (1E-10), to fit the calculator's display.

HP 41C Program: AGM

01 LBL^T AGM
02 STO 01
03 X<>Y
04 STO 02
05 X<>Y
06 LBL 00
07 RCL 02
08 RCL 01
09 ENTER
10 R↑
11 R↑
12 X<>Y
13 R↓
14 ENTER
15 R↑
16 +
17 2
18 /
19 STO 01
20 R↓
21 *
22 SQRT
23 STO 02
24 R↑
25 -
26 ABS
27 1E-10
28 X≤Y?
29 GTO 00
30 CLA
31 ^T AGM =
32 ARCL 01
33 AVIEW
34 END

HP 42S/Swiss Micros DM42/Free42 Program AGM:

00 {53-Byte Prgm}
01 LBL "AGM"
02 STO 01
03 X<>Y
04 STO 02
05 X<>Y
06 LBL 00
07 RCL 02
08 RCL 01
09 ENTER
10 R↑
11 R↑
12 X<>Y
13 R↓
14 ENTER
15 R↑
16 +
17 2
18 /
19 STO 01
20 R↓
21 *
22 SQRT
23 STO 02
24 R↑
25 -
26 ABS
27 1E-10
28 X≤Y?
29 GTO 00
30 CLA
31 "AGM = "
32 ARCL 01
33 AVIEW
34 END

The instructions for both the HP 41C and 42S versions are same: enter X and Y on the respective stacks and XEQ AGM.

Example (ALL/STD mode is applied):

AGM(37, 78):
37, 78, XEQ AGM returns:
Alpha: AGM = 55.5947005279


Blog Post: http://edspi31415.blogspot.com/2020/07/h...metic.html

[Werner]

In a lengthy thread about bore hole volumes not so long ago, I came up with the following, thanks to Albert Chan and Gerson W. Barbosa (stack-only, 41/42 compatible and reg Z preserved):

Code:

00 { 27-Byte Prgm } @   X       Y       Z       T
01▸LBL 02
02 X<>Y
03 X<> ST T
04 2
05 ÷
06▸LBL "AGM" @          a       b       z
07 STO ST T @           a       b       z       a
08 X<>Y
09 STO+ ST T @          b       a       z       2AM
10 STO× ST Y
11 X<>Y
12 SQRT @               GM      b       z       2AM
13 X≠Y?
14 GTO 02
15 R↓
16 END

I you want to show "AGM=.." at the end, for Free42 and DM42 there's a way that does not use the Alpha register: just add

LSTO "AGM"
VIEW "AGM"

at the end.

Cheers, Werner

[Albert Chan]

Originally, Werner's code return converged AM sequence.

(06-09-2020 01:14 PM)Albert Chan Wrote:  If we assume non-zero AGM arguments, AM sequence always converge.

We can shown GM sequence also converge. (even with rounding errors)

Assuming it does not, but alternate between lo, hi.
In other words, GM sequence = GM₁, GM₂, ..., lo, hi, lo, hi, ...

But, √(lo * hi) = √(hi * lo)

Thus, we have only 2 possibilities:
1). GM₁, GM₂, ..., lo, hi, lo, lo, ...
2). GM₁, GM₂, ..., lo, hi, hi, hi, ...

→ assumption were wrong, GM *will* converge.

Bonus: GM convergence does not require non-zero AGM arguments.

→ AGM(x, 0) = AGM(0, x) = 0

[Werner]

But the next GM in the sequence is not calculated as sqrt(lo*hi)? but as SQRT(GMi*AMi)
Werner

[Albert Chan]

(07-29-2020 12:53 PM)Werner Wrote:  But the next GM in the sequence is not calculated as sqrt(lo*hi)? but as SQRT(GMi*AMi)

My mistake. I wrongly assumed AM converged to lo or hi too.

Proof (2nd attempt):

Again, assume GM sequence do not converge, but alternate between lo, hi
GM sequence = GM1, GM2, ..., lo, hi, lo, hi, ...

For big enough i, such that GM(i) = lo, and AM(i) converged to g:

GM(i+1) = hi = √(g*lo)
GM(i+2) = lo = √(g*hi)

If hi > lo, we have √(g*lo) > √(g*hi), which is impossible, even with rounding errors.

→ when AM converged, GM sequence is non-decreasing, will converge too.