A not so useful HP-16C program

[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