+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (41/42) AGM - Arithmetic-Geometric Mean (/thread-15164.html)
(41/42) AGM - Arithmetic-Geometric Mean - Werner - 06-09-2020 06:48 AM
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 ENDCheers, Werner
RE: (41/42) AGM - Arithmetic-Geometric Mean - Albert Chan - 06-09-2020 01:14 PM
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.
RE: (41/42) AGM - Arithmetic-Geometric Mean - Albert Chan - 06-09-2020 03:51 PM
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, sNow, 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
RE: (41/42) AGM - Arithmetic-Geometric Mean - Albert Chan - 06-09-2020 11:55 PM
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 ENDRE: (41/42) AGM - Arithmetic-Geometric Mean - Gerson W. Barbosa - 06-10-2020 04:10 AM
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 ENDExample:
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).
RE: (41/42) AGM - Arithmetic-Geometric Mean - Werner - 06-10-2020 08:02 AM
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 ENDto 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
RE: (41/42) AGM - Arithmetic-Geometric Mean - Gerson W. Barbosa - 06-10-2020 11:23 AM
(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.
RE: (41/42) AGM - Arithmetic-Geometric Mean - Werner - 06-10-2020 01:16 PM
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 ENDCheers, Werner
RE: (41/42) AGM - Arithmetic-Geometric Mean - Gerson W. Barbosa - 06-10-2020 02:14 PM
(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 END30 bytes (26 bytes for 1-character global labels), 16 steps, 41-compatible.
RE: (41/42) AGM - Arithmetic-Geometric Mean - Werner - 06-10-2020 02:24 PM
(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