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 Add (sign%)

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)

CALL Add(1)

DO UNTIL firstword = words

CALL Divide(term(), denom&)

CALL Multiply(term(), denom& - 1)

CALL Divide(term(), 26)

CALL Add(1)

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)

CALL Add(-1)

DO UNTIL firstword = words

CALL Divide(term(), denom&)

CALL Multiply(term(), denom& - 1)

CALL Divide(term(), 57122)

CALL Add(-1)

denom& = denom& + 2

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

LOOP

CALL PrintOut(words)

END

'--------------------------------------------------------------------

SUB Add (sign)

SHARED words, firstword

IF sign = 1 THEN

'add it on

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