DECLARE SUB atan239 (denom&)
DECLARE SUB atan5 (denom&)
DECLARE SUB PrintOut (words%)

'Program to calculate pi, version 4.8
'A major rewrite of version 4.2, this uses only two arrays instead of
'three, and includes a host of speedups based on a similar C program.
'A sampler: all the carries are reserved until the end, the divide and
'add routines are combined, two terms are added at a time, and the number
'of function calls is minimized. It's a big change for a small gain, since
'the compiled version requires 28.6 seconds for 5000 digits on my 486 66MHz
'computer, a 10% gain over version 4.2; like before, it's capable of about
'150,000 digits of pi.
'
'This program has come a long way from version 1.0; thanks are due to
'Larry Shultis, Randall Williams, Bob Farrington and Adrian Umpleby.
'One final note for speed freaks: this program will run about 6 times faster
'if written in C using an optimizing compiler. Likewise, if you can figure
'out a way to do integer division and use both the quotient and remainder,
'this program can easily be sped up by 25%.      -jasonp@isr.umd.edu

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

words = digits& \ 4 + 3
DIM SHARED sum&(words + 1), term(words + 1)
start! = TIMER
                                        
                                         '---------------16*atan(1/5)
denom& = 3: firstword = 1: lastword = 2
sum&(1) = 3: term(1) = 3: sum&(2) = 2000: term(2) = 2000

DO UNTIL firstword >= words
   CALL atan5(denom&)
   denom& = denom& + 2
LOOP

                                         '------------   -4*atan(1/239)
denom& = 3: firstword = 2: remainder& = 4

FOR x = 2 TO words
   dividend& = remainder& * 10000              'crunch out 1st term
   term(x) = dividend& \ 239&
   remainder& = dividend& - term(x) * 239&
   sum&(x) = sum&(x) - term(x)
NEXT x

DO UNTIL firstword >= words
   CALL atan239(denom&)
   denom& = denom& + 4
LOOP
                                      

FOR x = words TO 2 STEP -1                        '-------finish up
   IF sum&(x) < 0 THEN                                  'release carries
      quotient& = sum&(x) \ 10000                       'and borrows
      sum&(x) = sum&(x) - (quotient& - 1) * 10000
      sum&(x - 1) = sum&(x - 1) + quotient& - 1
   END IF
   IF sum&(x) >= 10000 THEN
      quotient& = sum&(x) \ 10000
      sum&(x) = sum&(x) - quotient& * 10000
      sum&(x - 1) = sum&(x - 1) + quotient&
   END IF
NEXT x

CALL PrintOut(words)
PRINT "computation time: "; TIMER - start!; " seconds"
END

'------------------------------------------------------------------
SUB atan239 (denom&)
SHARED words, firstword

remainder1& = term(firstword)                        'first divide implicitly
remainder2& = 0: remainder3& = 0: remainder4& = 0
denom2& = denom& + 2: firstword = firstword + 1

FOR x = firstword TO words
   temp& = term(x)
   dividend& = remainder1& * 10000 + temp&
   temp& = dividend& \ 57121
   remainder1& = dividend& - temp& * 57121

   dividend& = remainder2& * 10000 + temp&
   temp2& = dividend& \ denom&
   remainder2& = dividend& - temp2& * denom&
   sum&(x) = sum&(x) + temp2&

   dividend& = remainder3& * 10000 + temp&
   temp& = dividend& \ 57121
   remainder3& = dividend& - temp& * 57121

   dividend& = remainder4& * 10000 + temp&
   temp2& = dividend& \ denom2&
   remainder4& = dividend& - temp2& * denom2&
   sum&(x) = sum&(x) - temp2&
   term(x) = temp&
NEXT x

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

END SUB

'-------------------------------------------------------------------
SUB atan5 (denom&)
SHARED words, firstword, lastword

FOR x = firstword TO lastword + 1
   temp& = term(x)
   dividend& = remainder1& * 10000 + temp&
   temp& = dividend& \ 25
   remainder1& = dividend& - temp& * 25&
   term(x) = temp&

   dividend& = remainder2& * 10000 + temp&
   temp& = dividend& \ denom&
   remainder2& = dividend& - temp& * denom&
   sum&(x) = sum&(x) - temp&
NEXT x

FOR x = lastword + 2 TO words
   dividend& = remainder2& * 10000
   temp& = dividend& \ denom&
   remainder2& = dividend& - temp& * denom&
   sum&(x) = sum&(x) - temp&
NEXT x

IF term(lastword + 1) > 0 AND lastword < words THEN lastword = lastword + 1
IF term(firstword) = 0 THEN firstword = firstword + 1

denom& = denom& + 2
remainder1& = 0: remainder2& = 0

FOR x = firstword TO lastword + 1
   temp& = term(x)
   dividend& = remainder1& * 10000 + temp&
   temp& = dividend& \ 25
   remainder1& = dividend& - temp& * 25&
   term(x) = temp&

   dividend& = remainder2& * 10000 + temp&
   temp& = dividend& \ denom&
   remainder2& = dividend& - temp& * denom&
   sum&(x) = sum&(x) + temp&
NEXT x

FOR x = lastword + 2 TO words
   dividend& = remainder2& * 10000
   temp& = dividend& \ denom&
   remainder2& = dividend& - temp& * denom&
   sum&(x) = sum&(x) + temp&
NEXT x

IF term(lastword + 1) > 0 AND lastword < words THEN lastword = lastword + 1
IF term(firstword) = 0 THEN firstword = firstword + 1

END SUB

'------------------------------------------------------------------
SUB PrintOut (words)
PRINT "pi = 3+."
FOR i = 1 TO words \ 3
   PRINT " ";
   PRINT RIGHT$("0000" + LTRIM$(STR$(sum&(3 * (i - 1) + 2))), 4);
   PRINT RIGHT$("0000" + LTRIM$(STR$(sum&(3 * (i - 1) + 3))), 4);
   PRINT RIGHT$("0000" + LTRIM$(STR$(sum&(3 * (i - 1) + 4))), 4);
   IF i MOD 5 = 0 THEN PRINT "  :"; 12 * i
NEXT i

PRINT " ";
FOR i = 3 * (words \ 3) + 2 TO digits
   PRINT RIGHT$("0000" + LTRIM$(STR$(sum&(i))), 4);
NEXT i
PRINT : PRINT

END SUB