The Museum of HP Calculators

HP Forum Archive 16

[ Return to Index | Top of Index ]

Yet another slice of pie
Message #1 Posted by Katie Wasserman on 30 Oct 2006, 1:20 a.m.

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.

      
Re: Yet another slice of pie
Message #2 Posted by Paul Dale on 30 Oct 2006, 2:34 a.m.,
in response to message #1 by Katie Wasserman

...and it doesn't work in real mode :-)

Still, a nice algorithm and program.

- Pauli

            
Re: Yet another slice of pie
Message #3 Posted by Katie Wasserman on 30 Oct 2006, 10:30 a.m.,
in response to message #2 by Paul Dale

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

      
Re: Yet another slice of pie
Message #4 Posted by Crawl on 30 Oct 2006, 6:47 a.m.,
in response to message #1 by Katie Wasserman

BTW, that formula is just the Taylor series of arcsine, evaluated at one-half.

            
Re: Yet another slice of pie
Message #5 Posted by Katie Wasserman on 30 Oct 2006, 10:45 a.m.,
in response to message #4 by Crawl

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

      
Re: Yet another slice of pie
Message #6 Posted by Namir on 30 Oct 2006, 9:29 a.m.,
in response to message #1 by Katie Wasserman

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.

            
Re: Yet another slice of pie
Message #7 Posted by Katie Wasserman on 30 Oct 2006, 10:51 a.m.,
in response to message #6 by Namir

Namir,

It sure looks a lot cleaner in BASIC, very pretty! Thanks for presenting it so that it's so understandable.

-Katie

      
The fastest algorithm for Pi????!!!!
Message #8 Posted by Namir on 30 Oct 2006, 9:13 p.m.,
in response to message #1 by Katie Wasserman

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.

            
Re: The fastest algorithm for Pi????!!!!
Message #9 Posted by Katie Wasserman on 31 Oct 2006, 12:52 a.m.,
in response to message #8 by Namir

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


[ Return to Index | Top of Index ]

Go back to the main exhibit hall