(41/42) AGM - Arithmetic-Geometric Mean

[Werner]

Upon request ;-)

Code:

00 { 28-Byte Prgm }
01▸LBL "AGM"
02 ENTER
03▸LBL 02 @     .       x       y       .
04 ENTER
05 - @          .       x       y       y
06 R↓
07 STO ST T @   x       y       y       x
08 STO- ST Z @  x       y       y-x     x
09 × @          x*y     y-x     x       x
10 SQRT
11 X<>Y
12 2
13 ÷ @         (y-x)/2 sqrt(x*y) x      x
14 R↑
15 STO+ ST Y @  x       x'      y'      x
16 X≠Y?
17 GTO 02
18 END

Cheers, Werner

[Albert Chan]

Thanks, Werner

Great idea on testing convergence of arithmetic mean.
I think this always work (no infinite loops).

My argument is like this (AM, GM for arithmetic and geometric means)

We know AM ≥ GM, but with rounding error, it is possible the rule fail.
At the final stage of AGM, with rounding error, we may have AM ≥ GM - 1ULP
(gap cannot be more than 1 ULP, otherwise next sequence will reduce it)

Example, with IEEE double, final stage of AGM:

AGM(a,b) = AGM((a+b)/2, √(ab)) = ...

AGM(1,4) ⇒ AGM(2.2430285802876027, 2.2430285802876031) ; AM = GM - 1 ULP
AGM(1,5) ⇒ AGM(2.6040081905309407, 2.6040081905309402) ; AM = GM + 1 ULP
AGM(1,6) ⇒ AGM(2.9513287423905727, 2.9513287423905727) ; AM = GM

This shows that testing convergence AM = GM will not work.

Note that, even with rounding error, AM sequence is decreasing.
(We required GM ≥ AM + 2 ULP to reverse this trend)

If we assume non-zero AGM arguments, AM sequence always converge.
This applied to intermediates too. Example, AGM(-1,1) should not converge (*)

Above (max of ± 1 ULP) may apply to complex numbers too.
Again, using IEEE double for real and imaginary parts.

AGM(1+1j, 1+5j):
AM ⇒ (1.1429851941908882 + 2.680170008949454j)
GM ⇒ (1.142985194190888 + 2.6801700089494545j) ; AM = GM + (1,-1) ULP

This suggest your setup also work for complex numbers

(*) Actually, AGM(-1,1) converged for IEEE double, when AM underflow to 0.0
I wonder testing for GM convergence may be better, quiting much earlier.

[Albert Chan]

Using your idea, this is my updated agm2(a,b)
Timings suggested it run just as fast as my old code, even with b=0 test included.
(For b=0, s value is useless)

Code:

from cmath import sqrt, pi

def agm2(a,b):
    t = -1.
    s = a*a
    while b:
        a1 = (a+b) * 0.5
        if a1 == a: break
        k = a1-a
        b = sqrt(a*b)
        a = a1
        t += t
        s += t*k*k
    return b, s

Now, redo EK(2), that AGM2 crashes, and had to rigged it to get it running

>>> b = 1j
>>> x, y = agm2(1+b, 2*sqrt(b))
>>> x,y
((1.1981402347355923+1.1981402347355923j), 1.8277863241778547j)

>>> k = pi/x
>>> e = k*y/4
>>> e             # E(m=2)
(0.59907011736779614+0.59907011736779614j)
>>> k             # K(m=2)
(1.3110287771460598-1.3110287771460598j)

FYI, this was my rig, which might not work in all cases ...

>>> x0, y0 = agm2(1, b)
>>> 2*x0, 2*(y0 + b*b)              # to get above x, y
((1.1981402347355923+1.1981402347355923j), 1.8277863241778547j)

P.S. how to convert above agm2 code for Free-42 ?

Update: I added updated agm2 code, that checked convergence of GM's.
It should produce same result as above code, since it also returns GM's.

Update2: redefined agm2, starting with s = 0 instead of a*a
this produces much more accurate E(m) - K(m), see here

[Albert Chan]

Hi, Werner

I modified your code to test for convergence of GM instead
This meant AGM(x,0) will quit right the way.

Line 8 and line 14 changed, line 15 deleted.

Code:

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

[Gerson W. Barbosa]

15 steps, 31 bytes:

Code:


00 { 31-Byte Prgm }
01▸LBL "AGM"
02 STO ST T
03 X<>Y
04 +
05 2
06 STO÷ ST Y
07 X<> ST L
08 X=Y?
09 GTO 01
10 RCL× ST T
11 SQRT
12 GTO "AGM"
13▸LBL 01
14 R↓
15 END

Example:

agm(1 + 2i, 3 + 4i) + agm(5 + 6i, 7 + 8i)

1 ENTER 2 ⬏ COMPLEX 3 ENTER 4 ⬏ COMPLEX XEQ AGM
5 ENTER 6 ⬏ COMPLEX 7 ENTER 8 ⬏ COMPLEX XEQ AGM +

→

7.835947925677309435700566636717429 +
9.886442230519331091330438898700769i

P.S.: 42-only (the 41 lacks recall arithmetic).

[Werner]

If I can do AM=(a+b)/2 instead of a+(b-a)/2, and I test for GM convergence, then:

Code:

00 { 22-Byte Prgm }
01▸LBL 02
02 R↑
03▸LBL "AGM" @    a       b
04 ENTER
05 RCL+ ST Z
06 2
07 ÷ @            AM      a       b       b
08 R↓
09 ×
10 SQRT
11 X≠Y? @         GM      b       AM      AM
12 GTO 02
13 END

to make it 41-compatible, replace lines 4 and 5 with
RCL ST Y
RCL ST Y
+
at the expense of 1 byte.

If you don't like that the global label gets excuted in the loop, put it before LBL 02 and place a Rv in between, again at the expense of 1 byte

Cheers, Werner

[Gerson W. Barbosa]

(06-10-2020 08:02 AM)Werner Wrote:  If you don't like that the global label gets excuted in the loop, put it before LBL 02 and place a Rv in between, again at the expense of 1 byte

No, that’s not an issue. But I would like the program to preserve the previous content of stack register X so that chained calculations involving AGM are possible without the use of numbered registers to store intermediate results.

Example:

agm(1, 2) + agm(3, 4)

1 ENTER 2 XEQ AGM 3 ENTER 4 XEQ AGM +

→

4.938818707406477275808395332421768

Anyway, the low byte count and least number of steps are very nice!

Gerson.

[Werner]

Seeing that yours really is only 29 bytes (if you replace the GTO "AGM" by a short numeric label), it was not easy to do better. But not impossible, and 41-compatible:

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

Cheers, Werner

[Gerson W. Barbosa]

(06-10-2020 01:16 PM)Werner Wrote:  Seeing that yours really is only 29 bytes (if you replace the GTO "AGM" by a short numeric label), it was not easy to do better. But not impossible, and 41-compatible:

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 00
15 R↓
16 END

Well done!

With a 1-character label, like “M”, mine is actually 27 bytes long, 15 steps (including LBL & END). I’ve labeled yours “W” here (for Werner, or Winner :-), 25 bytes, 16 steps and 41-compatible. That’s what I’ll keep on my 41C and 41CV.

Just a small typo in your listing: line 14 should be GTO 02 (or line 01 changed to LBL 00).

Regards,

Gerson.

——-

P.S.: Just for the record,

Code:


00 { 30-Byte Prgm }
01▸LBL "AGM"
02 STO ST T
03 X<>Y
04 +
05 2
06 STO÷ ST Y
07 X<> ST L
08 X=Y?
09 GTO 01
10 R↑
11 ×
12 SQRT
13 GTO "AGM"
14▸LBL 01
15 R↓
16 END

30 bytes (26 bytes for 1-character global labels), 16 steps, 41-compatible.

[Werner]

(06-10-2020 02:14 PM)Gerson W. Barbosa Wrote:  Just a small typo in your listing: line 14 should be GTO 02 (or line 01 changed to LBL 00).

Thanks, I corrected it!
Werner