AGM (HP-33S & HP-15C) « Next Oldest | Next Newest »

 ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-12-2010, 05:09 PM ``` 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. ▼ Lyuka Member Posts: 169 Threads: 12 Joined: Aug 2007 10-12-2010, 10:06 PM 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 ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-13-2010, 11:52 AM 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. ▼ Lyuka Member Posts: 169 Threads: 12 Joined: Aug 2007 10-13-2010, 08:05 PM 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 ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-14-2010, 06:38 AM 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. ▼ Lyuka Member Posts: 169 Threads: 12 Joined: Aug 2007 10-14-2010, 07:56 AM 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 David Hayden Senior Member Posts: 528 Threads: 40 Joined: Dec 2008 10-14-2010, 06:51 AM 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 Paul Dale Posting Freak Posts: 3,229 Threads: 42 Joined: Jul 2006 10-13-2010, 02:50 AM The 34s has an optional AGM function that does just this for real and complex arguments :-) - Pauli ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-13-2010, 11:54 AM 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. Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-13-2010, 04:43 AM 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 XY  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 Xy 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 ``` ▼ Jean-Michel Member Posts: 191 Threads: 41 Joined: Jun 2007 10-13-2010, 05:39 PM ```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. ``` ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-13-2010, 09:28 PM 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. Paul Dale Posting Freak Posts: 3,229 Threads: 42 Joined: Jul 2006 10-14-2010, 03:56 AM 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 ▼ Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-14-2010, 06:21 PM Very good of you Jean-Marc and Pauli. You did away with the tolerance issue and gave us a fast 41C implementation. Namir ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-14-2010, 07:55 PM 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 ▼ Francis Pierot Junior Member Posts: 7 Threads: 0 Joined: Oct 2010 10-15-2010, 04:19 AM More like A0014 GTO A on HP-33s ;-) ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-15-2010, 06:52 AM Fixed. Merci! :-) Werner Member Posts: 163 Threads: 7 Joined: Jul 2007 10-16-2010, 03:41 PM 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 ``` ▼ Werner Member Posts: 163 Threads: 7 Joined: Jul 2007 10-18-2010, 05:44 AM 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 ``` ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-18-2010, 06:43 PM 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 ``` Werner Member Posts: 163 Threads: 7 Joined: Jul 2007 10-14-2010, 04:56 AM 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 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 a g ‘IFTE( eps= END NIP »``` Regards, Gerson. Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-16-2010, 12:56 AM 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 ``` ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-16-2010, 04:41 PM 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). ``` ▼ Jean-Michel Member Posts: 191 Threads: 41 Joined: Jun 2007 10-16-2010, 05:13 PM 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. ▼ Lyuka Member Posts: 169 Threads: 12 Joined: Aug 2007 10-16-2010, 06:18 PM 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 Thomas Klemm Senior Member Posts: 735 Threads: 34 Joined: May 2007 10-16-2010, 06:23 PM 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 Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-16-2010, 07:31 PM Thank you, Jean-Marc! I should have imagined it :-) Best regards, Gerson. ▼ Jean-Michel Member Posts: 191 Threads: 41 Joined: Jun 2007 10-17-2010, 07:23 PM 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. Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 10-18-2010, 05:30 PM 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 ``` Marcus von Cube, Germany Posting Freak Posts: 3,283 Threads: 104 Joined: Jul 2005 10-14-2010, 05:31 PM 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 731 10-29-2013, 10:23 AM Last Post: Thomas Radtke Elliptic integrals of 1st and 2nd kind calculated by complex agm Gjermund Skailand 3 319 06-29-2013, 03:39 PM Last Post: Gjermund Skailand Maximum number of program steps in HP-42S, 33S, and 35S? Walter B 3 334 12-18-2012, 03:44 PM Last Post: Eric Smith Another technique of HP-33S label saving x34 6 419 10-11-2012, 11:09 AM Last Post: Luiz C. Vieira (Brazil) HP-33S Display Quirk Matt Agajanian 4 302 05-15-2012, 08:31 AM Last Post: Michael de Estrada Sudoku Solver for HP-33s Marcel Samek 3 287 05-11-2012, 03:01 PM Last Post: Lode WP34s: Complex AGM (Bug...) Eduardo Duenez 33 1,512 04-23-2012, 05:59 PM Last Post: Paul Dale HP 33s out of production? bill platt 36 1,793 03-06-2012, 05:04 AM Last Post: LMMT HP-33s bugs--Please clarify Matt Agajanian 3 240 03-06-2012, 01:45 AM Last Post: Thomas Radtke HP-33s Battery Switching--What's the gag? Matt Agajanian 0 144 03-03-2012, 07:09 AM Last Post: Matt Agajanian

Forum Jump: