A not so useful HP-16C program

[Gerson W. Barbosa]

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.

[Joe Horn]

Ok, I give up. What is it for? What does it do? How is it used? Some hints, please.

[Gerson W. Barbosa]

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

[rprosperi]

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.

[Gerson W. Barbosa]

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

[Dieter]

(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

[Gerson W. Barbosa]

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

This reminds me of e^pi – pi.

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.

[Gerson W. Barbosa]

(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.141592653589793238462643383279502884197169399375105820974944592307816406286208​9986280348253421170679
   82148086513282306647093844609550582231725359408128481117450284102701938521105559​64462294895493038196
   44288109756659334461284756482337867831652712019091456485669234603486104543266482​13393607260249141273
   72458700660631558817488152092096282925409171536436789259036001133053054882046652​13841469519415116094
   33057270365759591953092186117381932611793105118548074462379962749567351885752724​89122793818301194912
   98336733624406566430860213949463952247371907021798609437027705392171762931767523​84674818467669405132
   00056812714526356082778577134275778960917363717872146844090122495343014654958537​10507922796892589235
   42019956112129021960864034418159813629774771309960518707211349999998372978049951​05973173281609631859
   50244594553469083026425223082533446850352619311881710100031378387528865875332083​81420617177669147303
   59825349042875546873115956286388235378759375195778185778053217122680661300192787​6611195909216420199 

   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.141592653589793238462643383279502884197169399375105820974944592307816406286208​9986280348253421170679
   82148086513282306647093844609550582231725359408128481117450284102701938521105559​64462294895493038196
   44288109756659334461284756482337867831652712019091456485669234603486104543266482​13393607260249141273
   72458700660631558817488152092096282925409171536436789259036001133053054882046652​13841469519415116094
   33057270365759591953092186117381932611793105118548074462379962749567351885752724​89122793818301194912
   98336733624406566430860213949463952247371907021798609437027705392171762931767523​84674818467669405132
   00056812714526356082778577134275778960917363717872146844090122495343014654958537​10507922796892589235
   42019956112129021960864034418159813629774771309960518707211349999998372978049951​05973173281609631859
   50244594553469083026425223082533446850352619311881710100031378387528865875332083​81420617177669147303
   59825349042875546873115956286388235378759375195778185778053217122680661300192787​6611195909216420198 



============================================

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.141592653589793238462643383279502884197169399375105820974944592307816406286208​9986280348253421170679
   82148086513282306647093844609550582231725359408128481117450284102701938521105559​64462294895493038196
   44288109756659334461284756482337867831652712019091456485669234603486104543266482​13393607260249141273
   72458700660631558817488152092096282925409171536436789259036001133053054882046652​13841469519415116094
   33057270365759591953092186117381932611793105118548074462379962749567351885752724​89122793818301194912
   98336733624406566430860213949463952247371907021798609437027705392171762931767523​84674818467669405132
   00056812714526356082778577134275778960917363717872146844090122495343014654958537​10507922796892589235
   42019956112129021960864034418159813629774771309960518707211349999998372978049951​05973173281609631859
   50244594553469083026425223082533446850352619311881710100031378387528865875332083​81420617177669147303
   59825349042875546873115956286388235378759375195778185778053217122680661300192787​66111959092164202