HP Forums

Upon request ;-)


Cheers, Werner

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.

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)


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

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.

15 steps, 31 bytes:


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).

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


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 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.

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:


Cheers, Werner

(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,


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

(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