A not so useful HP-16C program - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html) +--- Forum: General Forum (/forum-4.html) +--- Thread: A not so useful HP-16C program (/thread-10552.html) |
A not so useful HP-16C program - Gerson W. Barbosa - 04-21-2018 11:37 PM Well, it will run on the HP-11C as well (only replace DSZ with DSE). 001- LBL A 002- STO I 003- 4 004- ENTER 005- ENTER 006- 2 007- ENTER 008- 4 009- LBL 0 010- ENTER 011- + 012- R↓ 013- × 014- √x 015- ENTER 016- R↓ 017- x↔y 018- + 019- LSTx 020- R⬆ 021- × 022- LSTx 023- R↓ 024- x↔y 025- ÷ 026- ENTER 027- + 028- ENTER 029- R↓ 030- R↓ 031- DSZ 032- GTO 0 033- R↓ 034- × 035- √x 036- R/S 037- ENTER 038- R⬆ 039- ÷ 040- ENTER 041- × 042- CHS 043- 1 044- + 045- √x 046- 3 047- × 048- CHS 049- 1 050- 6 051- + 052- × 053- R⬆ 054- + 055- + 056- 1 057- 5 058- ÷ 059- RTN No numbered registers, only the index register and the stack. No attempt has been made to make it shorter. Hopefully this might have some didactic value (to demonstrate an algorithm). Gerson. RE: A not so useful HP-16C program - Joe Horn - 04-22-2018 01:08 AM Ok, I give up. What is it for? What does it do? How is it used? Some hints, please. RE: A not so useful HP-16C program - Gerson W. Barbosa - 04-22-2018 01:14 AM (04-22-2018 01:08 AM)Joe Horn Wrote: Ok, I give up. What is it for? What does it do? How is it used? Some hints, please. Don’t give up, at least of yet. Enter an integer number. Start with 1 and run the program. Increase it to 2, then to 3 and see what you get in X and Y. Repeat this again pressing the R/S key afterwards. I’ll give an reference later. Painstakingly trying to fix errors on my smartphone while watching to Lost in Space on Netflix (Danger, Will Robinson!). Who said touchscreen was a good idea? RE: A not so useful HP-16C program - rprosperi - 04-22-2018 01:37 AM So it's easy to see the input is the number of iterations, but not so clear to see how the resulting X and Y converge. Look forward to the insights once folks have suffered long enough. RE: A not so useful HP-16C program - Gerson W. Barbosa - 04-22-2018 02:11 AM (04-22-2018 01:37 AM)rprosperi Wrote: So it's easy to see the input is the number of iterations, but not so clear to see how the resulting X and Y converge. Look forward to the insights once folks have suffered long enough. Back to my good old desktop computer! That's Archimedes' method to approximate Pi, except that instead of starting with the hexagon, I start with inscribed and circumscribed squares to a radius 1/2 circumference (thus saving a few steps as the initial constants, 2 and 4, don't involve surds). Without a positional number system and without a consistent math notation, that was such a feat, 0.6 digits per iteration! Not bad around 250 BC. Pressing R/S accelerates the convergence by a factor of 3. Basically, I use formula 2.6 in this paper in terms of a and b (perimeters of the inscribed and circumscribed n-gons to a radius 1/2 circumference ). I came up with a similar precision formula seven years ago, albeit an empirically-obtained one. Currently, I have a method that yields 3 times as much digits per iteration, compared to the basic Archimedes' algorithms, but I think it can go up to more than 8 times as much (5 digits per iteration). No intention to compete with modern methods, though. Square root extractions are a time-consuming task, be it done my hand or by machine. P.S.: BTW, the perimeter of the circumscribed 96-gon, Archimedes’ second bound, is about 21.9990021999/7 (well, actually 22.9990021975, but the former is nicer). This is also the place to discuss near integers like 2*(e - atan(e)) = 2.9999978 (that ‘s gonna be our 3 in that nerd’s clock!). RE: A not so useful HP-16C program - Dieter - 04-22-2018 01:08 PM (04-22-2018 02:11 AM)Gerson W. Barbosa Wrote: P.S.: BTW, the perimeter of the circumscribed 96-gon, Archimedes’ second bound, is about 21.9990021999/7 (well, actually 22.9990021975, but the former is nicer). This is also the place to discuss near integers like 2*(e - atan(e)) = 2.9999978 (that ‘s gonna be our 3 in that nerd’s clock!). This reminds me of e^pi – pi. And you should also read the caption. ;-) Dieter RE: A not so useful HP-16C program - Gerson W. Barbosa - 04-22-2018 05:59 PM (04-22-2018 01:08 PM)Dieter Wrote:(04-22-2018 02:11 AM)Gerson W. Barbosa Wrote: P.S.: BTW, the perimeter of the circumscribed 96-gon, Archimedes’ second bound, is about 21.9990021999/7 (well, actually 22.9990021975, but the former is nicer). This is also the place to discuss near integers like 2*(e - atan(e)) = 2.9999978 (that ‘s gonna be our 3 in that nerd’s clock!). I've tried to "fix" that at least a couple of times :-) \({e}^{\pi }-\pi +\frac{9^{2}}{89998-{10}^{5}\cdot \left ( {\frac{9^{2}}{89998}} \right )^{2}}=19.99999999999999295470\) \({e}^{\pi }-\pi+\left(\frac{3}{10^{2}}\right)^{2}+\frac{1}{\left ( \ln (2)\cdot 10^{4}+\frac{\sqrt{10}}{6} \right )^{2}}=20.00000000000000072951\) (04-22-2018 01:08 PM)Dieter Wrote: And you should also read the caption. ;-) "Also, I hear the 4th root of (9^2 + 19^2/22) is pi." Only slightly better, but many 2's and too many 9's, although not so much at 6's and 7's: \(\frac{2\left ( 16\sqrt{2}+1 \right )}{15+\frac{1}{24-\frac{9999}{2^{20+\frac{22552}{99999}}}}}=3.1415926535876\) Gerson. RE: A not so useful HP-16C program - Gerson W. Barbosa - 04-22-2018 09:31 PM (04-21-2018 11:37 PM)Gerson W. Barbosa Wrote: Hopefully this might have some didactic value (to demonstrate an algorithm). RPN, especially when making extensive use of the stack, is not adequate to demonstrate or share algorithms. So, here is a Decimal BASIC version: OPTION ARITHMETIC DECIMAL_HIGH ! 1000 digits precision LET r = TIME LET n = 6 ! number of sides of the first circumscribed polygon, a hexagon in this case LET m = 330 ! number of iterations ( number of sides of the last polygon: 6*2^m ) LET b = 2*SQR(3) ! perimeter of the circumscribed hexagon LET a = 3*b/4 ! perimeter of the inscribed triangle FOR i = 1 TO m LET a = SQR(a*b) ! current a = GM(previous a, current b); GM = geometric mean LET b = 2*a*b/(a + b) ! next b = HM(current a, previous b); HM = harmonic mean LET n = n + n ! double the number of sides at each iteration NEXT i LET a = SQR(a*b) ! final a LET t = a/n LET q = (2*b+a*(16-3*SQR(1-t*t)))/15 ! Chakrabarti-Hudson approximation to pi LET t = q/n LET c = t*t LET t = c*(329868000 + c*(42226800 + c*(4619230 + 481213*c))) LET u = 1164240000 + t LET p = 1/4656960000*(a*(1164240000 + t) + SQR(a*(-9313920000*b*(-1164240000 + t) + a*u*u))) ! my approximation to pi PRINT p PRINT TIME - r END 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679 8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196 4428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273 7245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094 3305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912 9833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132 0005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235 4201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859 5024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303 598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420199 1.55 This yields 999 digits of pi after 330 iterations. Notice that q is pi to about 600 digits, and a and b are pi to about 200 digits. The following is a more faithful version of the RPN program, which starts with a square instead of a hexagon: LET n = 4 ! start with a square LET m = 5 ! 5 iterations -> 128-gon LET b = 4 LET a = 2 FOR i = 1 TO m LET a = SQR(a*b) LET b = 2*a*b/(a+b) LET n = n + n NEXT i LET a = SQR(a*b) ! Archimedes method (128-gon lower bound) PRINT a LET t = a/n LET q = (2*b+a*(16-3*SQR(1-t*t)))/15 ! Chakrabarti-Hudson approximation (three-time acceleration) PRINT q END 3.14127725093278 3.14159265359634 ============================================ Update: (04-25-2018 04:47 PM) Here is a version without the Chakrabarti-Hudson approximation. It requires an intermediate step, but saves one square root extraction. Not sure whether this is more efficient, though. OPTION ARITHMETIC DECIMAL_HIGH ! 1000 digits precision LET s = TIME LET n = 6 ! number of sides of the first circumscribed polygon, a hexagon in this case LET m = 330 ! number of iterations ( number of sides of the last polygon: 6*2^m ) LET b = 2*SQR(3) ! perimeter of the circumscribed hexagon LET a = 3*b/4 ! perimeter of the inscribed triangle FOR i = 1 TO m LET a = SQR(a*b) ! current a = GM(previous a, current b); GM = geometric mean LET b = 2*a*b/(a + b) ! next b = HM(current a, previous b); HM = harmonic mean LET n = n + n ! double the number of sides at each iteration NEXT i LET a = SQR(a*b) ! a -> pi to 199 digits LET r = (2*a + b)/3 ! r -> pi to 399 digits LET t = r/n LET c = t*t LET q = a*(1 + (c*(60 + 7*c))/360) ! q -> pi to 598 digits LET t = q/n LET c = t*t LET t = c*(329868000 + c*(42226800 + c*(4619230 + 481213*c))) LET u = 1164240000 + t LET p = 1/4656960000*(a*(1164240000 + t) + SQR(a*(-9313920000*b*(-1164240000 + t) + a*u*u))) PRINT TIME - s;"seconds" PRINT p ! p -> pi to 999 digits END .83 seconds 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679 8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196 4428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273 7245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094 3305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912 9833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132 0005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235 4201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859 5024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303 598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198 ============================================ Update: (04-27-2018 05:43 PM) Another version which involves two cubic root extractions: OPTION ARITHMETIC DECIMAL_HIGH ! 1000 digits precision SUB CBR(x) ! Cubic root subroutine IF x<>0 THEN LET cb = EXP(LOG(x)/3) DO LET w = cb LET cb = (2*cb + x/(cb*cb))/3 LOOP UNTIL ABS(cb - w) < 1e-1000 LET x = cb ELSE LET x = 0 END IF END SUB LET s = TIME LET n = 6 ! number of sides of the first circumscribed polygon, a hexagon in this case LET m = 331 ! number of iterations ( number of sides of the last polygon: 6*2^m ) LET b = 2*SQR(3) ! perimeter of the circumscribed hexagon LET a = 3*b/4 ! perimeter of the inscribed triangle FOR i = 1 TO m LET a = SQR(a*b) ! current a = GM(previous a, current b); GM = geometric mean LET b = 2*a*b/(a + b) ! next b = HM(current a, previous b); HM = harmonic mean LET n = n + n ! double the number of sides at each iteration NEXT i LET a = SQR(a*b) ! a -> pi to 199 digits LET b2 = b*b LET n2 = n*n LET n4 = n2*n2 LET n6 = n2*n4 LET t = b*(30*n2 - 13*b2) + SQR(500*n6 + b2*(600*n4 + b2*(165*b2 - 720*n2))) CALL CBR(t) LET cr2 = 2 CALL CBR(cr2) ! cb2 = cubic root of 2 LET q = (b + (b2 - 5*n2)*cr2/t + t/cr2)/3 ! q -> pi to 599 digits LET t = q/n LET c = t*t LET t = c*(329868000 + c*(42226800 + c*(4619230 + 481213*c))) LET u = 1164240000 + t LET p = 1/4656960000*(a*(1164240000 + t) + SQR(a*(-9313920000*b*(-1164240000 + t) + a*u*u))) PRINT TIME - s;"seconds" PRINT p ! p -> pi to 996 digits END .91 seconds 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679 8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196 4428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273 7245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094 3305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912 9833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132 0005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235 4201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859 5024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303 5982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164202 |