'===========================================================================
' Subject: CALCULATE PI                       Date: 01-05-97 (07:59)       
'  Author: Jason Stratos Papadopoulos         Code: QB, QBasic, PDS        
'  Origin: comp.lang.basic.misc             Packet: ALGOR.ABC
'===========================================================================
DECLARE SUB PrintOut (sum%(), words%)
DECLARE SUB Multiply (term%(), words%, mult&, firstword%)
DECLARE SUB Divide (term%(), words%, denom&, firstword%)
DECLARE SUB Add (sum%(), term%(), words%, sign%, firstword%)
DECLARE SUB FastDivide (term%(), words%, denom&)

'Program to calculate pi, version 2.0
'The algorithm used is Gregory's series with Euler acceleration.
'This program uses the optimal Euler 2/3 rule: rather than use Euler's
'series for all the terms, compute instead 1/3 of the terms using
'Gregory's series and the rest using Euler's. It can be shown that
'each term in this compound series cuts the error by a factor of 3,
'while using only Euler's series has each term cut the error by a
'factor of 2. This is a major timesaver: it reduces the number of terms
'to be added up by over 35%, and of the terms that remain 1/3 can
'be crunched out faster than normal! The code also includes some tricks
'to speed things up (like reducing the size of the arrays Euler's series
'works on).
'
'Converging faster also means more digits can be computed. Some tests
'show the program is capable of computing about 51,000 digits of pi,
'and is quite fast if compiled (5000 digits in about 90 seconds on
'a 486 66MHz computer). I'd be grateful if someone can help me code
'the Divide and FastDivide SUBs in assembly, which can probably make
'the program twice as fast. Comments or questions to jasonp@wam.umd.edu

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

words = digits& \ 4 + 4
terms& = CLNG(digits& / .477) \ 3 + 1
IF terms& MOD 2 > 0 THEN terms& = terms& + 1
DIM sum(words), term(words)

                                          'Gregory's Series-------
PRINT TIME$: sum(1) = 1: denom& = 3: sign = -1

FOR x& = 1 TO terms& - 1
  
   CALL FastDivide(term(), words, denom&)
   CALL Add(sum(), term(), words, sign, 2)
   denom& = denom& + 2: sign = -sign

NEXT x&
                                          'Euler's Acceleration---
firstword = 2: x& = 1
CALL FastDivide(term(), words, 2 * denom&)

DO UNTIL firstword = words

   denom& = denom& + 2
   CALL Add(sum(), term(), words, sign, firstword)
   CALL Divide(term(), words, denom&, firstword)
   CALL Multiply(term(), words, x&, firstword)
  
   IF term(firstword) = 0 THEN firstword = firstword + 1
   x& = x& + 1

LOOP
                                          'Finish up--------------
CALL Add(sum(), term(), words, sign, firstword)
CALL Multiply(sum(), words, 4, 1)
CALL PrintOut(sum(), words)
END

'--------------------------------------------------------------------
SUB Add (sum(), term(), words, sign, 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(), words, denom&, 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 (term(), words, denom&)
'not really a fast divide, but there are fewer operations
'since dividend& below doesn't have term(x) added on (always 0)

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(), words, mult&, firstword)

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

END SUB

'------------------------------------------------------------------
SUB PrintOut (sum(), words)

PRINT : PRINT "pi=3."
i = 2
DO UNTIL i = words - 1
   j = sum(i)
   IF j > 999 THEN
       PRINT " " + RIGHT$(STR$(j), 4);
   ELSEIF j > 99 THEN
       PRINT " 0" + RIGHT$(STR$(j), 3);
   ELSEIF j > 9 THEN
       PRINT " 00" + RIGHT$(STR$(j), 2);
   ELSE
       PRINT " 000" + RIGHT$(STR$(j), 1);
   END IF
  
   IF (i - 1) MOD 15 = 0 THEN PRINT
   i = i + 1
LOOP
PRINT : PRINT : PRINT TIME$

END SUB