HP Forums
Yet another slice of pie - Printable Version

+- HP Forums (https://archived.hpcalc.org/museumforum)
+-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html)
+--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html)
+--- Thread: Yet another slice of pie (/thread-101774.html)



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