Yet another slice of pie



#10

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.


#11

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

Still, a nice algorithm and program.


- Pauli


#12

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

#13

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


#14

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

#15

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.


#16

Namir,

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

-Katie

#17

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.


#18

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: