here's a program to calculate pi, version 3.0
Author Message
here's a program to calculate pi, version 3.0

Hello again. Here's the next version of my QuickBasic program to
crunch out pi. Give it to anyone you want, I don't care.

jasonp

-------------CUT HERE INCLUDING THIS LINE--------------------------

DECLARE SUB PrintOut (words%)
DECLARE SUB Multiply (term%(), mult&)
DECLARE SUB Divide (term%(), denom&)
DECLARE SUB FastDivide (denom&)

'Program to calculate pi, version 3.0
'I've changed the algorithm to Machin's formula using Euler's atan series
'and cleaned up the subroutine calls. All in all, this version is about 25%
'faster than version 2.0; when compiled it requires 65 seconds for 5000
'digits on a 486 66MHz computer. Also, since Machin's formula converges much
'faster than Euler's, the program can find many more digits of pi before
'hitting overflow. In fact, I think the number of digits is limited only by
'the size of the 64K segment QuickBasic assigns for arrays (tops out around
'120,000 digits)
'
'{*filter*}s who want even more digits from this program can run QuickBasic
'with the /AH option for huge arrays, and change the program to work with
'3 digits at a time instead of 4. Just change any instance of 10000 to 1000,
'change \4 to \3 in the variable "words" below, and fiddle with the PrintOut
'subroutine.
'
'This program has come a long way from version 1.0; thanks are due to
'Christian Goldbach, Randall Williams, and Bob Farrington for good ideas.
'One final note for speed freaks: this program will run about 4 times faster
'if written in C using an optimizing comiler. I also plan to code the Divide
'and Multiply SUBs in 386 assembly sometime to get monster speed.

DEFINT A-Z
CLS
INPUT "how many digits"; digits&

words = digits& \ 4 + 4
DIM SHARED sum(words), term(words)

'--------------------16*atan(1/5)
PRINT TIME\$:  denom& = 3: firstword = 1
CALL FastDivide(26)
CALL Multiply(term(), 80)

DO UNTIL firstword = words
CALL Divide(term(), denom&)
CALL Multiply(term(), denom& - 1)
CALL Divide(term(), 26)
denom& = denom& + 2

IF term(firstword) = 0 THEN firstword = firstword + 1
LOOP

'-------------4*atan(1/239)
denom& = 3: firstword = 2
CALL FastDivide(57122)
CALL Multiply(term(), 956)

DO UNTIL firstword = words
CALL Divide(term(), denom&)
CALL Multiply(term(), denom& - 1)
CALL Divide(term(), 57122)
denom& = denom& + 2

IF term(firstword) = 0 THEN firstword = firstword + 1
LOOP

CALL PrintOut(words)
END

'--------------------------------------------------------------------
SHARED words, firstword

IF sign = 1 THEN

FOR x = words TO firstword STEP -1
sum(x) = sum(x) + term(x)
IF sum(x) >= 10000 THEN
sum(x - 1) = sum(x - 1) + 1
sum(x) = sum(x) - 10000
END IF
NEXT x

ELSE

'subtract it off
FOR x = words TO firstword STEP -1
sum(x) = sum(x) - term(x)
IF sum(x) < 0 THEN
sum(x - 1) = sum(x - 1) - 1
sum(x) = sum(x) + 10000
END IF
NEXT x

END IF
END SUB

'-------------------------------------------------------------------
SUB Divide (term(), denom&)
SHARED words, firstword

FOR x = firstword TO words
dividend& = remainder& * 10000 + term(x)
quotient = dividend& \ denom&
term(x) = quotient
remainder& = dividend& - quotient * denom&
NEXT x

END SUB

'---------------------------------------------------------------------
SUB FastDivide (denom&)
'not really a fast divide, but there are fewer operations
'since dividend& below doesn't have term(x) added on (always 0)

SHARED words
remainder& = 1
FOR x = 2 TO words
dividend& = remainder& * 10000
quotient = dividend& \ denom&
term(x) = quotient
remainder& = dividend& - quotient * denom&
NEXT x

END SUB

'---------------------------------------------------------------------
SUB Multiply (term(), mult&)
SHARED words, firstword

FOR x = words TO firstword STEP -1
product& = mult& * term(x) + carry
carry = product& \ 10000
term(x) = product& - carry * 10000&
NEXT x

END SUB

'------------------------------------------------------------------
SUB PrintOut (words)

PRINT : PRINT "pi=3."
i = 2
DO UNTIL i = words - 1

PRINT " " + RIGHT\$("000" + LTRIM\$(STR\$(sum(i))), 4);
IF (i - 1) MOD 15 = 0 THEN PRINT
i = i + 1

LOOP
PRINT : PRINT : PRINT TIME\$

END SUB

Tue, 06 Jul 1999 03:00:00 GMT
here's a program to calculate pi, version 3.0

Quote:

> >Hello again. Here's the next version of my QuickBasic program to
> >crunch out pi. Give it to anyone you want, I don't care.

> >jasonp

> Jason,  getting better with each new version.
> Thanks,
> D.C.

Now how do we write a program that will turn out
digits as fast as the one built into clisp.   I got about
51000 digits in less than 5 minutes.  Anyone know
where to get the source code for the German version
of clisp. It might be fun to try a basic program of that
pi generator.
Larry

Thu, 08 Jul 1999 03:00:00 GMT
here's a program to calculate pi, version 3.0

Quote:
> Hello again. Here's the next version of my QuickBasic program to
> crunch out pi. Give it to anyone you want, I don't care.

> jasonp

[snip]

Quote:

> 'Program to calculate pi, version 3.0
> 'I've changed the algorithm to Machin's formula using Euler's atan series
> 'and cleaned up the subroutine calls. All in all, this version is about
25%
> 'faster than version 2.0; when compiled it requires 65 seconds for 5000
> 'digits on a 486 66MHz computer. Also, since Machin's formula converges
much
> 'faster than Euler's, the program can find many more digits of pi before
> 'hitting overflow. In fact, I think the number of digits is limited only
by
> 'the size of the 64K segment QuickBasic assigns for arrays (tops out
around
> '120,000 digits)

Since any array will need to fit into lower memory and since variables
are limited to two bytes, there are some problems.  Because variables
are limited, even with huge arrays, you will have to use several dimensions
on for your array. I have worked with arrays of several million bits as
huge arrays. They had to be multiply dimensioned.
So you should be able to get somewhat more than 120000 digits.

Great program, it runs about 4 times faster than mine did.

Quote:
> '
> '{*filter*}s who want even more digits from this program can run QuickBasic
> 'with the /AH option for huge arrays, and change the program to work with
> '3 digits at a time instead of 4. Just change any instance of 10000 to
1000,
> 'change \4 to \3 in the variable "words" below, and fiddle with the
PrintOut
> 'subroutine.
> '
> 'This program has come a long way from version 1.0; thanks are due to
> 'Christian Goldbach, Randall Williams, and Bob Farrington for good ideas.

name is Larry Shultis. Goldbach is a handle I got for wasting so much
time  on the Goldbach conjecture.

Quote:
> 'One final note for speed freaks: this program will run about 4 times
faster
> 'if written in C using an optimizing comiler. I also plan to code the
Divide
> 'and Multiply SUBs in 386 assembly sometime to get monster speed.

[snip]

Thu, 08 Jul 1999 03:00:00 GMT
here's a program to calculate pi, version 3.0

: Hello again. Here's the next version of my QuickBasic program to
: crunch out pi. Give it to anyone you want, I don't care.

If you want, I'll convert it to C and run it on the nice fast computer
at the Univ here. On the other hand, the algorithms you're using here
are not the ones that would be chosen for a computation to many millions
of decimal places.

The problem is that computing n terms of the arctan series only provides
k.n correct digits of pi, for some value of k around 3 or 4. There are
much cleverer algorithms (Gauss-Salamin is the name for one of them,
Borwein's Quartic Iteration another), in which computing n terms gives
k.2^n correct digits (so a billion correct digits requires only about 30
terms).

These algorithms, however, require multiplication and division of long-numbers
by long-numbers; they are constrained by the speed of your multiplication
routines. It turns out that, for very long numbers, it is much faster to do

A' = FFT(A)
B' = FFT(B)
C = inverse_FFT(A' * B')

where A', B' are arrays, A'*B' is the array with c[i]=A'[i]*B'[i], and the
FFT is a standard operation taking time about n*log(n) which happens to
have the correct properties to make the above routine do multiplication.

Using FFTs, you can do some very clever things ... I have an FFT-based
program which will check that 2^216091-1 is prime within 4 hours, and
compute pi to *four million* decimal places in about the same time.

bits.

--
Tom

The Eternal Union of Soviet Republics lasted seven times longer than
the Thousand Year Reich

Fri, 09 Jul 1999 03:00:00 GMT
here's a program to calculate pi, version 3.0

I'll never again lose sleep over the value of Pi, cheers Jas.

Sun, 11 Jul 1999 03:00:00 GMT

 Page 1 of 1 [ 5 post ]

Relevant Pages