(41/42) AGM - Arithmetic-Geometric Mean
06-09-2020, 06:48 AM (This post was last modified: 06-09-2020 06:48 AM by Werner.)
Post: #1
 Werner Senior Member Posts: 648 Joined: Dec 2013
(41/42) AGM - Arithmetic-Geometric Mean
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
06-09-2020, 01:14 PM
Post: #2
 Albert Chan Senior Member Posts: 1,658 Joined: Jul 2018
RE: (41/42) AGM - Arithmetic-Geometric Mean
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.
06-09-2020, 03:51 PM (This post was last modified: 06-20-2020 04:27 PM by Albert Chan.)
Post: #3
 Albert Chan Senior Member Posts: 1,658 Joined: Jul 2018
RE: (41/42) AGM - Arithmetic-Geometric Mean
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
06-09-2020, 11:55 PM
Post: #4
 Albert Chan Senior Member Posts: 1,658 Joined: Jul 2018
RE: (41/42) AGM - Arithmetic-Geometric Mean
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
06-10-2020, 04:10 AM (This post was last modified: 06-10-2020 04:21 AM by Gerson W. Barbosa.)
Post: #5
 Gerson W. Barbosa Senior Member Posts: 1,419 Joined: Dec 2013
RE: (41/42) AGM - Arithmetic-Geometric Mean
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).
06-10-2020, 08:02 AM
Post: #6
 Werner Senior Member Posts: 648 Joined: Dec 2013
RE: (41/42) AGM - Arithmetic-Geometric Mean
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
06-10-2020, 11:23 AM
Post: #7
 Gerson W. Barbosa Senior Member Posts: 1,419 Joined: Dec 2013
RE: (41/42) AGM - Arithmetic-Geometric Mean
(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.
06-10-2020, 01:16 PM (This post was last modified: 06-10-2020 02:23 PM by Werner.)
Post: #8
 Werner Senior Member Posts: 648 Joined: Dec 2013
RE: (41/42) AGM - Arithmetic-Geometric Mean
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
06-10-2020, 02:14 PM (This post was last modified: 06-11-2020 01:00 AM by Gerson W. Barbosa.)
Post: #9
 Gerson W. Barbosa Senior Member Posts: 1,419 Joined: Dec 2013
RE: (41/42) AGM - Arithmetic-Geometric Mean
(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.
06-10-2020, 02:24 PM
Post: #10
 Werner Senior Member Posts: 648 Joined: Dec 2013
RE: (41/42) AGM - Arithmetic-Geometric Mean
(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
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)