Yet another slice of pie - Katie Wasserman - 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
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 - Paul Dale - 10-30-2006
...and it doesn't work in real mode :-)
Still, a nice algorithm and program.
- Pauli
Re: Yet another slice of pie - Crawl - 10-30-2006
BTW, that formula is just the Taylor series of arcsine, evaluated at one-half.
Re: Yet another slice of pie - Namir - 10-30-2006
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 - Katie Wasserman - 10-30-2006
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 - Katie Wasserman - 10-30-2006
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 - Katie Wasserman - 10-30-2006
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????!!!! - Namir - 10-30-2006
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????!!!! - Katie Wasserman - 10-31-2006
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
|