AGM (HP-33S & HP-15C)



#39

    Arithmetic-Geometric Mean (AGM)

HP-33S HP-15C

A0001 LBL A 001-42,21,11 f LBL A
A0002 R^ 002- 43 33 g R^
A0003 R^ 003- 43 33 g R^
A0004 STO A 004- 44 0 STO 0
B0001 LBL B 005-42,21, 0 f LBL 0
B0002 Rv 006- 33 Rv
B0003 Rv 007- 33 Rv
B0004 ENTER 008- 36 ENTER
B0005 Rv 009- 33 Rv
B0006 x<>y 010- 34 x<>y
B0007 + 011- 40 +
B0008 LASTx 012- 43 36 g LST x
B0009 R^ 013- 43 33 g R^
B0010 * 014- 20 *
BOO11 SQRTx 015- 11 SQRTx
B0012 x<>y 016- 34 x<>y
B0013 2 017- 2 2
B0014 / 018- 10 /
B0015 LOG 019- 43 13 g LOG
B0016 x<>y 020- 34 x<>y
B0017 LASTx 021- 43 36 g LST x
B0018 x<>y 022- 34 x<>y
B0019 LOG 023- 43 13 g LOG
B0020 x<>y 024- 34 x<>y
B0021 LASTx 025- 43 36 g LST x
B0022 Rv 026- 33 Rv
B0023 Rv 027- 33 Rv
B0024 - 028- 30 -
B0025 ABS 029- 43 16 g ABS
B0026 1E-11 030- 26 EEX
B0027 x<y? 031- 16 CHS
B0028 GTO B 032- 9 9
B0029 Rv 033-43,30, 8 g TEST 8
B0030 Rv 034- 22 0 GTO 0
B0031 RCL A 035- 33 Rv
B0032 x<>y 036- 33 Rv
B0033 RTN 037- 45 0 RCL 0
038- 34 x<>y
039- 43 32 g RTN
LN CK
A 12 7D30
B 123 7E2C

The program computes the arithmetic-geometric mean (AGM) of two positive real numbers. Due to rounding errors, there is some uncertainty in the last significant digit of the results. No optmization has been tried, so there is still room for improvement.

Examples:

1) AGM(3,4)

3 ENTER 4 XEQ A -> 3.48202767633 (The last digit should be 6)

2) AGM(12345,6789)

12345 ENTER 6789 XEQ A -> 9359.76103251 (The last digit should be 0)

3) sqrt(2*pi*sqrt(2*pi)/AGM(1,sqrt(2))) = Gamma(1/4)

2 pi * sqrt LASTx * 2 sqrt 1 XEQ A / sqrt -> 3.62560990822

4 1/x 1 - x! -> 3.62560990822

Gerson.


#40

Hi.

One of the most famous applications of AGM might be
Gauss-Legendre's AGM algorithm for Pi (which shows 2nd order convergence).

LBL "PI-GL"
1
STO 03
STO 02
STO 00
0.5
SQRT
STO 01
XEQ 00
XEQ 00
XEQ 00
XEQ 00
+
X^2
RCL/ 02
RTN
LBL 00
-
X^2
RCL* 03
STO- 02
RCL 03
RCL 00
RCL* 01
SQRT
RCL 00
RCL+ 01
0.5
*
STO 00
X<>Y
STO 01
RTN
END

REFERENCE:

J.Borwein and P.Borwein, Ramanujan and Pi, IEEE, Science and Applications, Supercomputing 88, Volume II, 117-128, 1988

Regards,

Lyuka


#41

Hello Lyuka,

I think there is a variation here:

http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv016.cgi?read=101118

(Message #11)

I'd like to do it in BASIC, on the HP-71B perhaps. I have implemented a kind of fixed-point arbitrary precision multiplication but the square root algorithm is still missing.

Regards,

Gerson.


#42

Hello Gerson,

How about these algorithm for the square root?

Quote:
Some sqrt(A) algorithm for 1 <= A <= 10, written in C-style.

/* Copyright (c) Takayuki HOSODA. All rights reserved.     */
/* sqrt2(a) , 2rd order convergence, Newton-Raphson method */
define sqrt2(a, epsilon) {
local x;

x = 0.8244 + 0.26072 * a; /* 1st order Chebyshev approximation, |err| <7.9% (1<=a<=10) */
while(abs(1 - a / (x * x)) > epsilon)
x = (x + a / x) * (1/2);
return(x);
}

/* sqrt3(a) , 3rd order convergence, Pade approximation */
define sqrt3(a, epsilon) {
local x, s;

x = 0.8244 + 0.26072 * a; /* 1st order Chebyshev approximation, |err| <7.9% (1<=a<=10) */
s = x * x;
while(abs(1 - a / s) > epsilon) {
x = x * (s + 3 * a) / (3 * s + a);
s = x * x;
}
return(x);
}

/* 1/sqrt(a) , 4th order convergence, taylor series, no division */
define sqrtinv4(a, epsilon) {
local x, h;

x = 1.0711 + a*(-0.17802 + a*0.0105597); /* 2nd order Chebyshev approximation, |err| <=9.7% (1<a<10) */
while(abs(h = 1 - a * x * x) > epsilon)
x *= 1 + h * (8 + h * (6 + h * 5)) * (1/16);
return(x);
}


Example: intermediate value and error

sqrt2(2, 1e-100)
1.415950375081733341 2.457732347050667709e-3
1.414214627565222532 1.506409720536721280e-6
1.414213562373496202 5.673167069204787484e-13
1.414213562373095049 8.046206148772844748e-26
1.414213562373095049 1.618535834713748356e-51
1.414213562373095049 6.549145620631325342e-103
1.414213562373095049

sqrt3(2, 1e-100)
1.414170564152398827 -6.080774244301628755e-5
1.414213562373085111 -1.405388097344160464e-14
1.414213562373095049 -1.734877563437081099e-43
1.414213562373095049 -3.263521730135614119e-130
1.414213562373095049

2*sqrtinv4(2, 1e-100)
7.070213571272571613e-1 -2.416011318629784659e-4
7.071067811865468656e-1 -1.863485008716057353e-15
7.071067811865475244e-1 -6.594648976029806730e-60
7.071067811865475244e-1 -1.034319719806940812e-237
1.414213562373095049

Regards,

Lyuka


#43

The functions need the following if statements to properly deal with a being 0:

if (a==0) return 0;

Otherwise, the functions look very cool!

Namir

Edited: 14 Oct 2010, 6:39 a.m.


#44

Hello Namir,

Thank you for your comment.

These functions are being optimized for use of decimal operation.

Therefore it assumes that the input 'A' is normalized within the range of 1 to 10, otherwise it will not function properly.

Regards,

Lyuka

#45

If anyone converts these to working C code, you'll need to change "(1/2)" to "(1./2)" and also "(1/16)" to "(1./16)" (or do something similar to ensure floating point division).

Without the change, the compiler will use integer division and create an integer result for both of these expressions. As a result, they will both be zero.

By adding the decimal point, you convert one argument to a floating point double, which causes the other argument, the divide operation, and the result to be floating point.

Dave

#46

The 34s has an optional AGM function that does just this for real and complex arguments :-)

- Pauli


#47

Quote:
The 34s has an optional AGM function that does just this for real and complex arguments :-)

Please let me know when it's available.

Regards,

Gerson.

#48

Here is the AGM version for the HP-41C:

LBL "AGM"
LBL A
"Y^X?"
PROMPT
LBL B
RCL Y
X<>Y
ST+ Z
*
0.5
ST* Z
Y^X
STO T
X<>Y
ST- T
RCL T
ABS
1E-10
CF 00
X<Y?
SF 00
RDN
RDN
FS? 00
GTO B
RTN

The program does not store any values in memory registers. Press A (in User mode) to get a prompt for the arguments. If you want to skip this prompt, enter the two values in the stack and then press B (in User mode). The program iterates and then displays the AGM in the X stack register. Notice that the tolerance for the solution is hard coded in the listing as 1E-10. I initially used 1E-11, but discovered that for some numbers, the iteration does not end, so I switched the tolerance value to 1E-10.

Edited: 13 Oct 2010, 7:51 a.m.


#49

Hello Namir,


Quote:
Notice that the tolerance for the solution is hard coded in the listing as 1E-10. I initially used 1E-11, but discovered that for some numbers, the iteration does not end, so I switched the tolerance value to 1E-10.

That's exactly what I was just testing. If you try AGM(4,5), you'll find the tolerance value has to be changed to 1E-08. Likewise for my second example, it will have to be changed to 1E-05. Finally, AGM(1E40,7E40) will loop forever unless the tolerance is switched to 1E32. But then the accuracy will be lowered in many cases. That's why I check the difference of the decimal logarithms of two consecutive means as a stop criterion. Perhaps there is a faster way, but that's the only one that has occurred to me. If you decide to change your program, I would suggest saving at least the stack register X to allow for easy one-level chained calculations.

Regards,

Gerson.

#50

How about using %ch to determine when to bail?

 01 LBL "AGM"
 02 CF 00
 03 LBL A
 04 "Y^X?"
 05 PROMPT
 06 LBL B
 07 RCL Y
 08 X<>Y
 09 ST+ Z
 10 *
 11 .5
 12 ST* Z
 13 Y^X
 14 %CH
 15 LASTX
 16 X<>Y
 17 ABS
 18 1 E-7
 19 X<Y?
 20 SF 00
 21 RDN
 22 RDN
 23 FS?C 00
 24 GTO B
 25 END

#51

Quote:
How about using %ch to determine when to bail?

Looks like a very good idea! :-)

Gerson.

#52

Agree ... %CH is a good way to go!

:-)

Namir

#53

Quote:
How about using %ch to determine when to bail?

Here is a new version using your suggestion. Faster and shorter! Incidentally all significant digits are correct now in examples 1 and 2.

Regards,

Gerson.

Arithmetic-Geometric Mean (AGM) - 2nd version

HP-33S HP-15C

A0001 LBL A 001-42,21,11 f LBL A
A0002 R^ 002- 43 33 g R^
A0003 R^ 003- 43 33 g R^
A0004 STO A 004- 44 0 STO 0
B0001 LBL B 005-42,21, 0 f LBL 0
B0002 Rv 006- 33 Rv
B0003 Rv 007- 33 Rv
B0004 ENTER 008- 36 ENTER
B0005 Rv 009- 33 Rv
B0006 x<>y 010- 34 x<>y
B0007 + 011- 40 +
B0008 LASTx 012- 43 36 g LST x
B0009 R^ 013- 43 33 g R^
B0010 * 014- 20 *
BOO11 SQRTx 015- 11 SQRTx
B0012 x<>y 016- 34 x<>y
B0013 2 017- 2 2
B0014 / 018- 10 /
B0015 %CHG 019- 43 15 g Delta%
B0016 LASTx 020- 43 36 g LST x
B0017 x<>y 021- 34 x<>y
B0018 ABS 022- 43 16 g ABS
B0019 1E-9 023- 26 EEX
B0020 x<y? 024- 16 CHS
B0021 GTO B 025- 7 7
B0022 Rv 026-43,30, 8 g TEST 8
B0023 Rv 027- 22 0 GTO 0
B0024 RCL A 028- 33 Rv
B0025 x<>y 029- 33 Rv
B0026 RTN 030- 45 0 RCL 0
031- 34 x<>y
032- 43 32 g RTN
LN CK
A 12 7D30
B 102 E161


#54

Hi,

here is another version for the HP-41

01 LBL "AGM"
02 ENTER^
03 LBL 01
04 CLX
05 RCL Z
06 RCL Y
07 +
08 2
09 /
10 RCL Y
11 SQRT
12 R^
13 SQRT
14 *
15 R^
16 X#Y?
17 GTO 01
18 END ( 30 bytes )

Though the tolerance is zero,
it seems there is no infinite loop.

Regards,
Jean-Marc.


#55

Hi Jean-Marc,

Works great! The flying goose never makes to the end of its short path and most of the times it will stop only a bit past halfway the display.

Regards,

Gerson.

#56

A one step shortening and one less square root per cycle:

01 LBL "AGM"
02 ENTER^
03 LBL 00
04 CLX
05 RCL Z
06 RCL Y
07 +
08 2
09 /
10 RCL Y
11 R^
12 *
13 SQRT
14 R^
15 X#Y?
16 GTO 00
17 END


- Pauli


#57

Very good of you Jean-Marc and Pauli. You did away with the tolerance issue and gave us a fast 41C implementation.

Namir


#58

Agreed!

    HP-33S                   HP-15C

A0001 LBL A 001-42,21,11 f LBL A
A0002 ENTER 002- 36 ENTER
A0003 Rv 003- 33 Rv
A0004 x<>y 004- 34 x<>y
A0004 + 005- 40 +
A0006 LASTx 006- 43 36 g LST x
A0007 R^ 007- 43 33 g R^
A0008 * 008- 20 *
AOO09 SQRTx 009- 11 SQRTx
A0010 x<>y 010- 34 x<>y
A0011 2 011- 2 2
A0012 / 012- 10 /
A0013 x>y? 013-43,30, 7 g TEST 7
A0014 GTO A 014- 22 11 GTO A
A0015 RTN 015- 43 32 g RTN

These are my short 33s and 15C versions. The difference to my first draft is only the test, x>y instead of x#y, which sometimes would cause endless loops. I wrongly assumed a tolerance factor might be a better solution. Anyway, I think this has been a good exercise ...

Still under test: AGM(7,8) loops forever on the HP-33s (but not on the 15C). Will test on the 42S.

Regards,

Gerson.

--------

Doesn't work on the HP-42S either, for (7,8); does work on Free42 Decimal, however. Perhaps an accuracy issue. The intermediate means increase towards the theoretical result so that after a few iterations two successive means should be the same. In most cases the one in register X gets slightly greater than the one in register Y (when they not converge to the same result), but in the (7,8) example this never happens.


Edited: 15 Oct 2010, 6:52 a.m. after one or more responses were posted


#59

More like A0014 GTO A on HP-33s ;-)


#60

Fixed.

Merci! :-)

#61

One byte shorter still:

01*LBL "AGM"
02 ENTER
03*LBL 02
04 ENTER
05 +
06 RDN
07 STO T
08 STO+ Z
09 *
10 SQRT
11 X<>Y
12 2
13 /
14 RUP
15 X#Y?
16 GTO 02
17 END

#62

More accurate, at the expense of only 2 bytes:

00 { 28-Byte Prgm }
01*LBL "AGM"
02 ENTER
03*LBL 02
04 ENTER
05 +
06 RDN
07 STO T
08 STO- Z
09 *
10 SQRT
11 X<>Y
12 2
13 /
14 RUP
15 STO+ Y
16 X#Y?
17 GTO 02
18 END

#63

Hello Werner,

This will run on the HP-42S only. The test in step 21 is made between successive arithmetic means, which are computed as a + (b - a)/2, as in your program. The original stack register X is preserved, but surely not a good idea trying to use only the stack. It would be much better to use a memory register instead, that's what they are meant for :-)

Regards,

Gerson.

HP-42S

00 { 44-Byte Prgm }
01*LBL "AGM"
02 R^
03*LBL 00
04 Rv
05 X=Y?
06 GTO 01
07 STO ST T
08 X<>Y
09 -
10 LASTX
11 RCL* ST T
12 SQRT
13 X<>Y
14 +/-
15 STO ST L
16 STO+ ST L
17 STO* ST X
18 RCL/ ST L
19 R^
20 STO+ ST Y
21 X#Y?
22 GTO 00
23 R^
24 X<>Y
25 R^
26*LBL 01
27 Rv
28 END


Examples:

1) pi x^2 16 ENTER 3 1/x y^x * -> 24.8698446781

3 sqrt 1 - 2 sqrt 2 * / 1 XEQ AGM -> 5.6747127659E-1

/ 3 1/x y^x 3 sqrt sqrt / -> 2.67893853471

3 1/x GAMMA -> 2.67893853471


2) 1 +/- sqrt pi XEQ AGM -> 1.43599262486 i8.82710446224E-1

#64

Gerson, you will not have been able to calulate example #3 as stated, as your programs use all four stack levels?

Here's my 42S implementation:

00 { 32-Byte Prgm }
01*LBL "AGM"
02 ENTER
03 ENTER
04*LBL 02
05 +
06 RDN
07 STO* ST Z
08 +
09 2
10 /
11 X<>Y
12 SQRT
13 RCL ST Y
14 %CH
15 1e-9
16 X<Y?
17 GTO 02
18 RUP
19 END

Cheers, Werner


#65

Hello Werner,

You can calculate the example #3 as stated if you store the stack register X in a memory register and recall it when AGM is done:

00 { 38-Byte Prgm }
01*LBL "AGM"
02 X<> ST Z (or RCL ST Z)
03 STO 00
04 RDN
05 ENTER
06 ENTER
07*LBL 02
08 +
09 RDN
10 STO* ST Z
11 +
12 2
13 /
14 X<>Y
15 SQRT
16 RCL ST Y
17 %CH
18 1E-9
19 X<Y?
20 GTO 02
21 RUP
22 RCL 00
23 X<>Y
24 END

Regards,

Gerson.

------------

P.S.: One step and one byte shorter:

00 { 37-Byte Prgm }
01*LBL "AGM"
02 X<> ST Z
03 STO 00
04 RDN
05 ENTER
06 ENTER
07*LBL 02
08 +
09 RDN
10 STO* ST Z
11 +
12 2
13 /
14 X<>Y
15 SQRT
16 RCL ST Y
17 %CH
18 1E-9
19 X<Y?
20 GTO 02
21 RCL 00
22 LASTX
23 END


Edited: 14 Oct 2010, 9:08 a.m.


#66

Dear all RPN’s enthusiasts !

I have to tell you that I am into studying your RPN programs with great interest. With great pleasure, I am playing to decrypt all this RPN codes and I am trying to understand how they differ or match from one HP model to another.

Behind all this line numbers, labels, tests, stack movements, I come to the surprising conclusion that they all scrub to converge to the expected result by recurrence, such as Imay have done it on an RPL calculator:

« -> a g 
‘IFTE( eps<a-g,
AGM((a+g)/2,SQRT(a*g)),
a)’
»
‘AGM’ STO

Where eps=1E-9 or whatever the available/expected precision.

OR, with RPL not so far from RPN stack manipulation :

« DO DUP2 + 2 / ROT ROT * SQRT UNTIL DUP2 - eps <= END » 


Edited: 16 Oct 2010, 12:52 p.m. after one or more responses were posted


#67

Hello,

I haven't been able to test your recursive version. Your second program will not leave the loop for a=1e7 and b=1e13, unless eps is greater than 10. By using %CH instead of subtraction, as per Mark Storkamp's suggestion, a fixed tolerance value will be enough:

« 
DO DUP2 * SQRT ROT ROT
+ 2 /
UNTIL DUP2 %CH
.000000001 <
END NIP
»

The following, which is a direct translation of my short versions, will loop forever in some cases, like 4 and 5. I haven't found non-working examples on the 15C, however.

« 
DO DUP2 * SQRT ROT ROT
+ 2 / DUP2
UNTIL >=
END NIP
»

Regards,

Gerson.

#68

Here is another attempt on the HP-42S, using ontly the stack and preserving register X, no tolerance factor. However occasionally the consecutive means will not converge to the same value as they should, perhaps due to the extra arithmetic in lines 15 and 16. I will test it on the HP-41 later.

00 { 34-Byte Prgm }
01*LBL "AGM"
02 R^
03*LBL 00
04 Rv
05 STO ST T
06 X<>Y
07 +
08 LASTX
09 R^
10 *
11 LASTX
12 Rv
13 *
14 LASTX
15 STO+ ST X
16 STO/ ST Y
17 CLX
18 LASTX
19 SQRT
20 R^
21 X#Y?
22 GTO 00
23 R^
24 24<>Y
25 END

Example:

sqrt(2*pi*sqrt(2*pi)/AGM(1,sqrt(2))) = Gamma(1/4)

2 pi * SQRT LASTx * 2 sqrt 1 XEQ AGM / SQRT -> 3.62560990822

4 1/x GAM -> 3.62560990822


#69

So far no problem on the HP-41, but I would have to check more examples. In these HP-42S and HP-41 versions the comparison is made between two consecutive AGMs, as is Jean-Marc Baillard's method, not between the arithmetic and geometric means. However, his program computes geometric means as sqrt(a)*sqrt(b), instead of sqrt(a*b), any special reason?


01*LBL'AGM
02 R^
03*LBL 00
04 RDN
05 STO T
06 X<>Y
07 +
08 LASTX
09 R^
10 *
11 LASTX
12 RDN
13 *
14 LASTX
15 STO+ X
16 STO/ Y
17 CLX
18 LASTX
19 SQRT
20 R^
21 X#Y?
22 GTO 00
23 R^
24 X<>Y
25 END

Examples:

1) AGM(3,4)
3 ENTER 4 XEQ AGM -> 3.4820277

2) AGM(12345,6789)
12345 ENTER 6789 XEQ AGM -> 9359.761031 (The last digit should be 3)

3) 2 pi * SQRT LASTx * 2 sqrt 1 XEQ AGM / SQRT -> 3.625609909 (The last digit should be 8)

4) AGM(1E07,1E13)

EEX 7 ENTER XEQ AGM -> 1.033295938E12 (Ok; on my shorter HP-15C version, which uses a different method, this example doesn't exit the loop).


#70

Hi Gerson,



In fact, I first used sqrt(a.b)

but in order to avoid an overflow if a.b > 9.999999999 E99

I finaly replaced sqrt(a.b) by sqrt(a).sqrt(b)



Best Regards,

Jean-Marc.


#71

Hi Jean-Marc,

And also in order to avoid an overflow if a + b > 9.999999999 E99
(a + b) / 2 would better be replaced by

a < b ? a + (b - a) / 2 : b + (a - b) / 2;

Regards,

Lyuka

#72

As the AGM is linear in its parameter (i.e. AGM(ka, kb) = k AGM(a, b)) I don't consider that a big issue. We could be happy to just calculate AGM(1, x) with 0 <= x <= 1 (assuming only real values).

Cheers

Thomas

#73

Thank you, Jean-Marc! I should have imagined it :-)

Best regards,

Gerson.


#74

Hi Gerson,



As Lyuka suggests it, (a+b)/2 may also overflow.

But practically, taking sqrt(a.b) and (a+b)/2

will work almost always.

So, Pauli's listing is better.



Best Regards,

Jean-Marc.

#75

By inserting a X<>Y instruction between lines 19 and 20 the program gets more stable on the real 42S:

.
.
.
20 X<>Y
21 R^
22 X#Y?
23 GTO 00
24 R^
25 X<>Y
26 END
#76

It's easy to implement AGM on the HP-38/39/40 series with the sequence applet.

In symb view:

U1(1)=A
U1(2)=(A+B)/2
U1(N)=(U1(N-1)+U2(N-1))/2
U2(1)=B
U2(2)=SQR(A*B)
U2(N)=SQR(U1(N-1)*U2(N-1))
On the home screen store the start values in A and B and then you can see the values in num view or plot the sequences in plot view. The graphics resolution isn't very good so you have to fiddle with the zoom options to see something. In num view, the number of digits is limited to 7. I don't know of a way to change that(*) but you can enter U1(6) or the like on the home screen to see a certain value.

At least it works.

(*) Just found out that you can move the cursor around and see the values on the bottom of the screen.


Edited: 15 Oct 2010, 4:55 a.m.


Possibly Related Threads...
Thread Author Replies Views Last Post
  Is the HP 33s the next hot HP collectible ? Michael de Estrada 12 701 10-29-2013, 10:23 AM
Last Post: Thomas Radtke
  Elliptic integrals of 1st and 2nd kind calculated by complex agm Gjermund Skailand 3 307 06-29-2013, 03:39 PM
Last Post: Gjermund Skailand
  Maximum number of program steps in HP-42S, 33S, and 35S? Walter B 3 324 12-18-2012, 03:44 PM
Last Post: Eric Smith
  Another technique of HP-33S label saving x34 6 411 10-11-2012, 11:09 AM
Last Post: Luiz C. Vieira (Brazil)
  HP-33S Display Quirk Matt Agajanian 4 294 05-15-2012, 08:31 AM
Last Post: Michael de Estrada
  Sudoku Solver for HP-33s Marcel Samek 3 276 05-11-2012, 03:01 PM
Last Post: Lode
  WP34s: Complex AGM (Bug...) Eduardo Duenez 33 1,460 04-23-2012, 05:59 PM
Last Post: Paul Dale
  HP 33s out of production? bill platt 36 1,735 03-06-2012, 05:04 AM
Last Post: LMMT
  HP-33s bugs--Please clarify Matt Agajanian 3 230 03-06-2012, 01:45 AM
Last Post: Thomas Radtke
  HP-33s Battery Switching--What's the gag? Matt Agajanian 0 139 03-03-2012, 07:09 AM
Last Post: Matt Agajanian

Forum Jump: