Yet another slice of pie « Next Oldest | Next Newest »

 ▼ Katie Wasserman Posting Freak Posts: 1,477 Threads: 71 Joined: Jan 2005 10-30-2006, 01:20 AM In looking for a 'natural' way to compute pi on a 16C I came across the following formula attributed to Newton: ```pi/6 = 1/2 + (1)/(2*3)/(2^3) + (1*3)/(2*4*5/(2^5) + (1*3*5)/(2*4*6*7)/(2^7) + ... ``` Which can be written in a more algorithmic-like form: ```pi * 1000....0 = 3 * 1000...0 * (x(0) + x(1) + x(2) + ...) where X(0) = 1; i=1 X(n) = X(n-1) * (i*i)/(i+1)/(i+2)/2/2; i=i+2 ``` Here it is implemented on a 16C using no registers other than 'I' and no assumptions about number base, word size or complement mode. The first 20 lines of the program are just to get the largest power of 10 (times 3) that will fit into the current word size without knowing the current number base. ``` LBL A 1 ; STO I ; starting value of i SL ; LST x ; 3 (base 10) + ; LBL 9 1 ; SL ; ENTER ; 10 (base 10) SL ; SL ; + ; X<>Y ; SWAP x ; times F? 5 ; Check for overflow - overflow is used to determine max word length GTO 8 GTO 9 ; keep multiplying by 10 LBL 8 LST x ; Restore maximal value 300....0 LBL 7 ; Top of main loop LST x ; X(n-1) RCL I ; x ; RCL I ; x ; *(i*i) ISZ ; RCL I ; / ; /(i+1) ISZ ; RCL I ; / ; /(i+2) SR SR ; /2/2 x=0 GTO 6 ; check for end, nothing more to add in + ; accumulate X(n) GTO 7 ; keep going... LBL 6 RDN ; roll down to result RTN ``` Convergence is so-so at about .6 digits per iteration. With a 64-bit word size, unsigned mode the calculator gives 19 digits of pi as 3141592653589793226 (the last 2 are wrong due to truncation) in about 90 seconds. It's not fast, but it is portable and fairly simple. Edited: 30 Oct 2006, 1:23 a.m. ▼ Paul Dale Posting Freak Posts: 3,229 Threads: 42 Joined: Jul 2006 10-30-2006, 02:34 AM ...and it doesn't work in real mode :-) Still, a nice algorithm and program. - Pauli ▼ Katie Wasserman Posting Freak Posts: 1,477 Threads: 71 Joined: Jan 2005 10-30-2006, 10:30 AM It doesn't take much to modify it for floating point mode, but still using integers: ```LBL A 1 STO I 3 EEX ; 3,000,000,000 is maximal -- no obvious way to check for overflow 9 x LBL 7 LST x RCL I x RCL I x ISZ RCL I / ISZ RCL I / 4 / 1 x>y ; check for less than 1; is same as 0 check in binary mode GTO 6 RND ; drop the '1' + GTO 7 LBL 6 RDN ; drop the '1' RDN RTN ``` Crawl Senior Member Posts: 306 Threads: 3 Joined: Sep 2009 10-30-2006, 06:47 AM BTW, that formula is just the Taylor series of arcsine, evaluated at one-half. ▼ Katie Wasserman Posting Freak Posts: 1,477 Threads: 71 Joined: Jan 2005 10-30-2006, 10:45 AM Yes indeed. Although Newton did find this formula for arcsin he actually (according to A History of PI) used a more complicated, faster converging formula to calculate pi = 3.1415926535897928... (last 2 digits wrong due to rounding). Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-30-2006, 09:29 AM Katie, Another cool algorithm!!! I implemented in QBASIC (for a real quick test) and it does real well. Took 21 iterations in QBASIC to match all of the digits of pi (calculate using 4 * ATN(1)) except the last one. Namir PS: Here is the QBASIC listing: ```' Estimate PI, as suggested by Kat Wasserman. ' HP Museum post on 10/30/2006 ' ' In looking for a 'natural' way to compute pi on a 16C I came across the following formula attributed to Newton: ' ' pi/6 = 1/2 + (1)/(2*3)/(2^3) + (1*3)/(2*4*5/(2^5) + (1*3*5)/(2*4*6*7)/(2^7) + ... ' ' Which can be written in a more algorithmic-like form: ' ' pi * 1000....0 = 3 * 1000...0 * (x(0) + x(1) + x(2) + ...) ' where X(0) = 1; i=1 ' X(n) = X(n-1) * (i*i)/(i+1)/(i+2)/2/2; i=i+2 ' ' DIM X AS DOUBLE DIM I AS INTEGER, J AS DOUBLE DIM PI AS DOUBLE, PI2 AS DOUBLE PI2 = 4 * ATN(1) X = 1 PI = X I = 1 J = 0 DO J = J + 1 X = X * I * I / (I + 1) / (I + 2) / 4 PI = PI + X ' actually update the estimate for Pi/3 I = I + 2 LOOP UNTIL ABS(3 * PI - PI2) < 1E-15 CLS PI = 3 * PI PRINT "Iterations = "; J PRINT "Estimated Pi = "; PI PRINT "Actual Pi = "; PI2 PRINT "Error = "; (PI - PI2) PRINT "% Error = "; 100 * (PI - PI2) / PI2 ``` Edited: 30 Oct 2006, 9:32 a.m. ▼ Katie Wasserman Posting Freak Posts: 1,477 Threads: 71 Joined: Jan 2005 10-30-2006, 10:51 AM Namir, It sure looks a lot cleaner in BASIC, very pretty! Thanks for presenting it so that it's so understandable. -Katie Namir Posting Freak Posts: 2,247 Threads: 200 Joined: Jun 2005 10-30-2006, 09:13 PM Katie, I found some fast algorithms to calculate Pi here Namir PS: Here is the QBASIC listing for the fast algorithm to calculate Pi. It akes 2 iterations to achieve very very good accuracy. ```DECLARE FUNCTION LnFact# (N AS INTEGER) ' Srinivasa Ramanujan Algorithm for fast calculations of Pi DIM Sum AS DOUBLE, K AS INTEGER DIM Fact1 AS DOUBLE, fact2 AS DOUBLE DIM Power AS DOUBLE, Term AS DOUBLE DIM Pi AS DOUBLE, Pi2 AS DOUBLE DIM T1 AS DOUBLE, T2 AS DOUBLE, T3 AS DOUBLE, T4 AS DOUBLE CLS Sum = 1103 K = 0 DO K = K + 1 T1 = LnFact(4 * K) T2 = LnFact(K) T3 = LOG(1103# + 26390# * K) T4 = K * LOG(396#) Term = T1 + T3 - 4 * (T2 + T4) Term = EXP(Term) Sum = Sum + Term PRINT Sum, Term LOOP UNTIL ABS(Term) < 1E-11 Pi = 1 / (2# * SQR(2#) / 9801# * Sum) Pi2 = 4 * ATN(1) PRINT "Iterations = "; K PRINT "Estimated Pi = "; Pi PRINT "Actual Pi = "; Pi2 PRINT "Error = "; (Pi - Pi2) PRINT "%Error = "; 100 * (Pi - Pi2) / Pi2 END FUNCTION LnFact# (N AS INTEGER) DIM Res AS DOUBLE DIM I AS INTEGER Res = 0 FOR I = 2 TO N Res = Res + LOG(I) NEXT I LnFact = Res END FUNCTION ``` Edited: 30 Oct 2006, 11:04 p.m. ▼ Katie Wasserman Posting Freak Posts: 1,477 Threads: 71 Joined: Jan 2005 10-31-2006, 12:52 AM That's fast, but I sure wouldn't want to have to implement that on a 16C in integer mode! There's a lot more about this algorithm and related ones on mathworld

 Possibly Related Threads... Thread Author Replies Views Last Post Happy National Pie Day Mark Hardman 2 310 01-23-2011, 07:11 PM Last Post: Martin Pinckney Pie Recipe Chuck 1 193 10-23-2010, 01:59 PM Last Post: Gerson W. Barbosa One more slice of PI GWB 11 509 11-27-2004, 05:40 PM Last Post: Larry Corrado Another Slice of PI Bill (Smithville, NJ) 3 253 05-13-2004, 11:13 AM Last Post: Gene

Forum Jump: 