The Museum of HP Calculators

HP Forum Archive 19

 AGM (HP-33S & HP-15C)Message #1 Posted by Gerson W. Barbosa on 12 Oct 2010, 5:09 p.m. ``` 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 xy 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.

 Re: AGM for Pi (Free42 Decimal)Message #2 Posted by Lyuka on 12 Oct 2010, 10:06 p.m.,in response to message #1 by Gerson W. Barbosa 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

 Re: AGM for Pi (Free42 Decimal)Message #3 Posted by Gerson W. Barbosa on 13 Oct 2010, 11:52 a.m.,in response to message #2 by Lyuka Hello Lyuka, I think there is a variation here: (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.

 Re: algorithm for the square rootMessage #4 Posted by Lyuka on 13 Oct 2010, 8:05 p.m.,in response to message #3 by Gerson W. Barbosa 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 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

 Re: algorithm for the square rootMessage #5 Posted by Namir on 14 Oct 2010, 6:38 a.m.,in response to message #4 by Lyuka 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.

 Re: algorithm for the square rootMessage #6 Posted by Lyuka on 14 Oct 2010, 7:56 a.m.,in response to message #5 by Namir 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

 Re: algorithm for the square rootMessage #7 Posted by David Hayden on 14 Oct 2010, 6:51 a.m.,in response to message #4 by Lyuka 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

 Re: AGM (HP-33S & HP-15C)Message #8 Posted by Paul Dale on 13 Oct 2010, 2:50 a.m.,in response to message #1 by Gerson W. Barbosa The 34s has an optional AGM function that does just this for real and complex arguments :-) - Pauli

 Re: AGM (HP-33S & HP-15C)Message #9 Posted by Gerson W. Barbosa on 13 Oct 2010, 11:54 a.m.,in response to message #8 by Paul Dale 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.

 Re: AGM (HP-33S & HP-15C)Message #10 Posted by Namir on 13 Oct 2010, 4:43 a.m.,in response to message #1 by Gerson W. Barbosa 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

 Re: AGM (HP-33S & HP-15C)Message #11 Posted by Gerson W. Barbosa on 13 Oct 2010, 8:32 a.m.,in response to message #10 by Namir 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.

 Re: AGM (HP-33S & HP-15C)Message #12 Posted by Mark Storkamp on 13 Oct 2010, 10:01 a.m.,in response to message #10 by Namir 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

 Re: AGM (HP-33S & HP-15C)Message #13 Posted by Gerson W. Barbosa on 13 Oct 2010, 11:18 a.m.,in response to message #12 by Mark Storkamp Quote: How about using %ch to determine when to bail? Looks like a very good idea! :-) Gerson.

 Re: AGM (HP-33S & HP-15C)Message #14 Posted by Namir on 13 Oct 2010, 11:25 a.m.,in response to message #12 by Mark Storkamp Agree ... %CH is a good way to go! :-) Namir

 Re: AGM (HP-33S & HP-15C)Message #15 Posted by Gerson W. Barbosa on 13 Oct 2010, 3:37 p.m.,in response to message #12 by Mark Storkamp 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 xy 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 ```

 Re: AGM (HP-33S & HP-15C)Message #16 Posted by JMBaillard on 13 Oct 2010, 5:39 p.m.,in response to message #15 by Gerson W. Barbosa ```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. ```

 Re: AGM (HP-33S & HP-15C)Message #17 Posted by Gerson W. Barbosa on 13 Oct 2010, 9:28 p.m.,in response to message #16 by JMBaillard 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.

 Re: AGM (HP-33S & HP-15C)Message #18 Posted by Paul Dale on 14 Oct 2010, 3:56 a.m.,in response to message #16 by JMBaillard 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

 Re: AGM (HP-33S & HP-15C)Message #19 Posted by Namir on 14 Oct 2010, 6:21 p.m.,in response to message #18 by Paul Dale Very good of you Jean-Marc and Pauli. You did away with the tolerance issue and gave us a fast 41C implementation. Namir

 Re: AGM (HP-33S & HP-15C)Message #20 Posted by Gerson W. Barbosa on 14 Oct 2010, 7:55 p.m.,in response to message #19 by Namir 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

 Re: AGM (HP-33S & HP-15C)Message #21 Posted by Francis Pierot on 15 Oct 2010, 4:19 a.m.,in response to message #20 by Gerson W. Barbosa More like A0014 GTO A on HP-33s ;-)

 Re: AGM (HP-33S & HP-15C)Message #22 Posted by Gerson W. Barbosa on 15 Oct 2010, 6:52 a.m.,in response to message #21 by Francis Pierot Fixed. Merci! :-)

 Re: AGM (HP-33S & HP-15C)Message #23 Posted by Werner on 16 Oct 2010, 3:41 p.m.,in response to message #18 by Paul Dale 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 ```

 Re: AGM (HP-33S & HP-15C)Message #24 Posted by Werner on 18 Oct 2010, 5:44 a.m.,in response to message #23 by Werner 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 ```

 Re: AGM (HP-33S & HP-15C)Message #25 Posted by Gerson W. Barbosa on 18 Oct 2010, 6:43 p.m.,in response to message #24 by Werner 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 ```

 Re: AGM (HP-33S & HP-15C)Message #26 Posted by Werner on 14 Oct 2010, 4:56 a.m.,in response to message #1 by Gerson W. Barbosa 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

 Re: AGM (HP-33S & HP-15C)Message #27 Posted by Gerson W. Barbosa on 14 Oct 2010, 7:41 a.m.,in response to message #26 by Werner 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 XY 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

 Re: AGM (HP-33S & HP-15C)Message #28 Posted by C.Ret on 14 Oct 2010, 11:49 a.m.,in response to message #27 by Gerson W. Barbosa 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

 Re: AGM (HP-33S & HP-15C)Message #29 Posted by Gerson W. Barbosa on 15 Oct 2010, 11:14 a.m.,in response to message #28 by C.Ret 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.

 Re: AGM (HP-33S & HP-15C)Message #30 Posted by Gerson W. Barbosa on 16 Oct 2010, 12:56 a.m.,in response to message #27 by Gerson W. Barbosa 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 ```

 Re: AGM (HP-33S & HP-15C)Message #31 Posted by Gerson W. Barbosa on 16 Oct 2010, 4:41 p.m.,in response to message #30 by Gerson W. Barbosa 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). ```

 Re: AGM (HP-33S & HP-15C)Message #32 Posted by JMBaillard on 16 Oct 2010, 5:13 p.m.,in response to message #31 by Gerson W. Barbosa 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.

 Re: AGM (HP-33S & HP-15C)Message #33 Posted by Lyuka on 16 Oct 2010, 6:18 p.m.,in response to message #32 by JMBaillard 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

 Re: AGM (HP-33S & HP-15C)Message #34 Posted by Thomas Klemm on 16 Oct 2010, 6:23 p.m.,in response to message #32 by JMBaillard 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

 Re: AGM (HP-33S & HP-15C)Message #35 Posted by Gerson W. Barbosa on 16 Oct 2010, 7:31 p.m.,in response to message #32 by JMBaillard Thank you, Jean-Marc! I should have imagined it :-) Best regards, Gerson.

 Re: AGM (HP-33S & HP-15C)Message #36 Posted by JMBaillard on 17 Oct 2010, 7:23 p.m.,in response to message #35 by Gerson W. Barbosa 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.

 Re: AGM (HP-33S & HP-15C)Message #37 Posted by Gerson W. Barbosa on 18 Oct 2010, 5:30 p.m.,in response to message #30 by Gerson W. Barbosa 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 ```

 Re: AGM (HP-33S & HP-15C)Message #38 Posted by Marcus von Cube, Germany on 14 Oct 2010, 5:31 p.m.,in response to message #1 by Gerson W. Barbosa 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.

Go back to the main exhibit hall