[net.math] Calculating Pi to Excessive Accuracy

tstorm@vu44.UUCP (Theo van der Storm) (08/03/85)

In article <423@ist.UUCP> bco@ist.UUCP (Brian Collins) writes:
>Following Ken Yap's posting, which soaks up 2 hours of CPU to
>calculate Pi to 1000 decimal places, may I offer the
>following sh/bc script that takes just over 6 mins on my
>Pyramid.
>
>bc <<!
>	/* Computes Pi to $1 places */
>	t = 3 * 10^${1-100}
>	for (i = 1; t != 0; i++) {
>		p = p + t
>		t = t * (2 * i - 1) * (2 * i - 1) / (8 * i * (2 * i + 1))
>	}
>	p
>!
>
>	...differential of asin( x) ;	asin( 0.5) = 	pi / 6 ...
>
>Any advance on 6 minutes?
>
>Brian Collins		...mcvax!ist!bco (bco@ist)
>Imperial Software Technology, London

Using Ken Yap's formula (!) here is a bc program
(written by me), which is almost twice as fast as YOUR program:
[Replace =+, =/, etc. by +=, /=, etc. if necessary]

bc << EOF
n = $1
t = 16 * (10^n) / 5
e = 1
s = 0
while (t != 0){
    s =+ t / e
    t =/ -25
    e =+ 2
}
t = 4 * (10^n) / 239
e = 1
while (t != 0){
    s =- t / e
    t =/ -57121
    e =+ 2
}
s

Apparently, Ken Yap's implementation was not very efficient :-)
-- 
	Theo van der Storm,    052 19'08"N / 004 51'16"E
UUCP:	       {seismo|decvax|philabs}!mcvax!vu44!tstorm
UUCP:	{ark|botter|klipper|tjalk|vu45|vu60}!vu44!tstorm

bco@ist.UUCP (Brian Collins) (08/07/85)

Following Ken Yap's posting, which soaks up 2 hours of CPU to
calculate Pi to 1000 decimal places, may I offer the
following sh/bc script that takes just over 6 mins on my
Pyramid.
[ I can't be bothered to shar it just for 9 lines! ]

bc <<!
	/* Computes Pi to $1 places */
	t = 3 * 10^${1-100}
	for (i = 1; t != 0; i++) {
		p = p + t
		t = t * (2 * i - 1) * (2 * i - 1) / (8 * i * (2 * i + 1))
	}
	p
!

The cost is guaranteed better than quadratic in the number of
decimal places, since each term is less than a quarter of the
previous one.

The recurrence relation is easily deduced from first principles
by:
	integrating the
		binomial expansion of the
			differential of
				asin( x)
	one term at a time, and using this to compute
		asin( 0.5)
	which is (of course)
		pi / 6

Any advance on 6 minutes?

-- 
Brian Collins		...mcvax!ist!bco (bco@ist)
Imperial Software Technology, London