A not so useful HP-16C program
04-22-2018, 09:31 PM (This post was last modified: 04-27-2018 05:53 PM by Gerson W. Barbosa.)
Post: #8
 Gerson W. Barbosa Senior Member Posts: 1,537 Joined: Dec 2013
RE: A not so useful HP-16C program
(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
 « Next Oldest | Next Newest »

 Messages In This Thread A not so useful HP-16C program - Gerson W. Barbosa - 04-21-2018, 11:37 PM RE: A not so useful HP-16C program - Joe Horn - 04-22-2018, 01:08 AM RE: A not so useful HP-16C program - Gerson W. Barbosa - 04-22-2018, 01:14 AM RE: A not so useful HP-16C program - rprosperi - 04-22-2018, 01:37 AM RE: A not so useful HP-16C program - Gerson W. Barbosa - 04-22-2018, 02:11 AM RE: A not so useful HP-16C program - Dieter - 04-22-2018, 01:08 PM RE: A not so useful HP-16C program - Gerson W. Barbosa - 04-22-2018, 05:59 PM RE: A not so useful HP-16C program - Gerson W. Barbosa - 04-22-2018 09:31 PM

User(s) browsing this thread: 1 Guest(s)