'===========================================================================
' Subject: GREGORY'S SERIES                  Date: 02-28-99 (23:41)       
' Author:  Jason Stratos Papadopoulos        Code: QB, QBasic, PDS        
' Origin:  jasonp@Glue.umd.edu             Packet: ALGOR.ABC
'===========================================================================
DECLARE SUB lastmult (digits&)
DECLARE SUB trickmult (j&)
DECLARE SUB addterm (j&)
DECLARE SUB printout (sum&())
CLS
DEFLNG A-Z

' You may be familiar with Gregory's series:
'
'    pi       1   1   1   1
'    -- = 1 - - + - - - + - - ...
'     4       3   5   7   9
'
' It's a prime example of reducing a mysterious number like pi to boring
' arithmetic. Unfortunately it's also a prime example of a *slow* infinite
' series! Since the terms alternate and decrease in size, the error involved
' in stopping at a certain point is bounded by the first term ignored; to get
' 5 digits of pi you'd have to add up 2*10^5 terms of the series, since each
' extra decimal means 10 times the accuracy means 10 times the work. For this
' reason nobody uses Gregory's series to compute anything, until now.
'
' This program lets you peer into the sum of Gregory's series far deeper than
' is ordinarily possible, by allowing one to add up enormous numbers of terms
' quickly. The starting point is a strange series for pi/4:
' 
'  pi   5     11  1*2      17  1*2   3*4      23  1*2   3*4    5*6
'  -- = -  -  -- -----  +  -- ----- -----  -  -- ----- ----- -------  + ...
'  4    6     14 3*5*2     22 3*5*2 7*9*2     30 3*5*2 7*9*2 11*13*2
'
' which can be found by a transformation of Gregory's series. Each term in
' this series is about 1/8 the size of the previous one, so it converges at
' about .9 digits per term...not bad, but not spectacular. What's very sur-
' prising is that the above can be slightly modified to converge to a
' *finite* sum of the terms in Gregory's series. For n even,
'
'        1   1   1          1         2n   / 5     11 1*2*(n-2) 
'    1 - - + - - - + ... - ----  =   ---- |  -  -  -- ---------- 
'        3   5   7         2n-1      2n+1  \ 6     14 3*5*(2n+3)
'
'    17 1*2*(n-2)  3*4*(n-4)      23 1*2*(n-2)  3*4*(n-4)   5*6*(n-6)      \
' +  -- ---------- ----------  -  -- ---------- ---------- ----------- +... |
'    22 3*5*(2n+3) 7*9*(2n+5)     30 3*5*(2n+3) 7*9*(2n+5) 11*13*(2n+7)    /
'
' Of course, as n becomes infinite the second series reduces to the first.
' It looks like a huge mess, but it's pretty neat: the above says that you
' can dial in a value of n, and you have a series that converges to the sum
' of n terms of Gregory's series at .9 digits per term, no matter how big
' n is. Hence 100 terms can get 90 digits of the sum, whether n is a thousand
' or a trillion!
'
' Combined with some multiple precision arithmetic and a few tricks for
' dividing by huge numbers, the result is pretty impressive. This program
' can calculate (to tens of thousands of digits) the sum of an incredible
' number of terms in Gregory's series, 10^5 or 10^50 or 10^1000 or whatever;
' and it doesn't need the age of the universe to do so. So the next time
' some C programmer gives you lip about how slow Basic is, ask him/her to
' time adding up 10^20 terms of Gregory's series, to 1000 digits of accuracy;
' my 486 using the compiled version of this program needs 59 seconds.
'
' A well-documented phenomenon you'll find is that when adding up numbers
' of terms that are multiples of 2 and 5, what you'll see is one wrong digit,
' then lots of correct ones, then one wrong digit, then lots of correct ones,
' and so on for a long way past the decimal point.

PRINT "This program adds up 4 * Gregory's series (1 - 1/3 + 1/5 - 1/7 +...)"
INPUT "How many digits do you wish to carry the sum to"; digits: PRINT
PRINT "How many terms do you wish added up? Express your answer as"
PRINT "two numbers, m and e, where the number of terms is m * 10^e;"
PRINT "so for 600,000 terms type ' 6 , 5 '. Please make m<10 and e>=5"
start:
INPUT mantissa%, exponent%: PRINT

IF exponent% < 5 or mantissa% > 9 THEN
   PRINT "Read the instructions, and try again..."
   GOTO start
END IF

SELECT CASE exponent% MOD 4                'divide exponent by 4; if
   CASE 0                                  'it isn't evenly divisible by 4,
      exponent% = exponent% \ 4            'absorb the difference into the
      mantissa% = mantissa% * 2            'mantissa
   CASE 1
      exponent% = (exponent% - 1) \ 4
      mantissa% = mantissa% * 20
   CASE 2
      exponent% = (exponent% - 2) \ 4
      mantissa% = mantissa% * 200
   CASE 3
      exponent% = (exponent% - 3) \ 4
      mantissa% = mantissa% * 2000
END SELECT

digits = digits \ 4 + 3                                      'initialize
DIM SHARED sum(digits + 1)
DIM SHARED term1(digits + 1)
DIM SHARED term2(digits + 1)
DIM SHARED term3(digits + 1)
j = 1: firstword = 1: term1(1) = 1
sum(1) = 3: start! = TIMER

FOR i = 2 TO digits: sum(i) = 3333: NEXT i                   'first term

WHILE firstword < digits                               'main computation
    CALL addterm(j)
    j = j + 2
WEND

CALL lastmult(digits)

FOR j = digits TO 2 STEP -1                'release carries and borrows
   IF sum(j) < 0 THEN
      quotient = sum(j) \ 10000
      sum(j - 1) = sum(j - 1) + quotient - 1
      sum(j) = sum(j) - 10000 * quotient + 10000
   END IF
   IF sum(j) >= 10000 THEN
      quotient = sum(j) \ 10000
      sum(j - 1) = sum(j - 1) + quotient
      sum(j) = sum(j) - 10000 * quotient
   END IF
NEXT j

finish! = TIMER
CALL printout(sum())
PRINT "Computation finished in "; finish! - start!; "seconds"
END

'------------------------------------------------------------------------
SUB addterm (j)
SHARED firstword, digits

CALL trickmult(j)

denom1 = 4 * j - 1: denom2 = denom1 + 2: denom3 = denom1 + 4
remainder1 = 0: remainder2 = 0
remainder3 = 0: remainder4 = 0

FOR i = firstword TO digits
      tempterm = term1(i)                         'multiply by (2j-1)/(4j-1)
      dividend = remainder1 * 10000 + tempterm
      temp = dividend \ denom1
      remainder1 = dividend - denom1 * temp
      tempterm = tempterm - temp

      dividend = remainder2 * 10000 + tempterm    'multiply by (2j+1)/(4j+3)
      temp = dividend \ denom2
      remainder2 = dividend - denom2 * temp
      tempterm = tempterm - temp

      dividend = remainder3 * 10000 + tempterm                  'divide by 8
      tempterm = dividend \ 8
      remainder3 = dividend - 8 * tempterm
     
      dividend = remainder4 * 10000 + tempterm        'scale the term and 
      temp = dividend \ denom3                        'subtract from the sum
      remainder4 = dividend - denom3 * temp
      sum(i) = sum(i) - temp - tempterm - tempterm - tempterm
      term1(i) = tempterm
NEXT i
 
denom1 = denom1 + 4: denom2 = denom1 + 2: denom3 = denom1 + 4
remainder1 = 0: remainder2 = 0
remainder3 = 0: remainder4 = 0

CALL trickmult(j + 1)

FOR i = firstword TO digits
      tempterm = term1(i)                         'multiply by (2j+3)/(4j+7)
      dividend = remainder1 * 10000 + tempterm
      temp = dividend \ denom1
      remainder1 = dividend - denom1 * temp
      tempterm = tempterm - temp

      dividend = remainder2 * 10000 + tempterm    'multiply by (2j+5)/(4j+9)
      temp = dividend \ denom2
      remainder2 = dividend - denom2 * temp
      tempterm = tempterm - temp

      dividend = remainder3 * 10000 + tempterm                 'divide by 8
      tempterm = dividend \ 8
      remainder3 = dividend - 8 * tempterm
  
      dividend = remainder4 * 10000 + tempterm          'scale the term and 
      temp = dividend \ denom3                          'add to the sum
      remainder4 = dividend - denom3 * temp
      sum(i) = sum(i) + temp + tempterm + tempterm + tempterm
      term1(i) = tempterm
NEXT i

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

END SUB

'-------------------------------------------------------------------------
SUB lastmult (digits)
SHARED mantissa%, exponent%

' as the name implies, this subroutine performs the final scaling
' operation, by multiplying sum() by 1 - 1/(m*10^e + 1). The code looks
' quite a bit like the code in trickmult, but is much simpler since no
' checking is needed for overflow.

sourceword = 1: remainder = 0
destword = 1 + exponent%

FOR i = destword TO digits
   dividend = remainder * 10000 + sum(sourceword)
   temp = dividend \ mantissa%
   remainder = dividend - temp * mantissa%
   term2(sourceword) = temp
   sourceword = sourceword + 1
NEXT i

sourceword = 1
FOR i = destword TO digits
   sum(i) = sum(i) - term2(sourceword)
   sourceword = sourceword + 1
NEXT i

sourceword = 1: remainder = 0
destword = destword + exponent%

DO UNTIL destword > digits
  
   FOR i = destword TO digits
      temp = term2(sourceword)
      dividend = remainder * 10000 + temp
      temp = dividend \ mantissa%
      remainder = dividend - temp * mantissa%
      sum(i) = sum(i) + temp
      term2(sourceword) = temp
      sourceword = sourceword + 1
   NEXT i

   sourceword = 1: remainder = 0
   destword = destword + exponent%
 
   FOR i = destword TO digits
      temp = term2(sourceword)
      dividend = remainder * 10000 + temp
      temp = dividend \ mantissa%
      remainder = dividend - temp * mantissa%
      sum(i) = sum(i) - temp
      term2(sourceword) = temp
      sourceword = sourceword + 1
   NEXT i

   sourceword = 1: remainder = 0
   destword = destword + exponent%
                 
LOOP

END SUB

'------------------------------------------------------------------------
SUB printout (sum())
SHARED digits

PRINT "the sum is 3+."
FOR i = 1 TO digits \ 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 * (digits \ 3) + 2 TO digits
   PRINT RIGHT$("0000" + LTRIM$(STR$(sum(i))), 4);
NEXT i
PRINT : PRINT

END SUB

'--------------------------------------------------------------------------
SUB trickmult (j)
SHARED mantissa%, exponent%, digits, firstword

' this subroutine multiplies term1() by 1 - (6j+1)/(m * 10^e + 2j+1). It
' isn't easy, since m*10^e can be a huge number like 10^50; rather than using
' multiple precision long division (which is quite slow) I use a dirty trick
' that relies on the huge number being a power of 10; in that case, dividing
' by 10^50 is easy, just shift the decimal point 50 places to the left (sym-
' bolically)! Of course, this program works in base 10000 instead of 10, but
' the principle's the same. One last dirty trick is to expand the multiplier
' as a geometric series:
'
'      6j+1        6j+1   /     2j+1       2j+1   2     2j+1   3     \
' ------------- = ------ | 1 - ------ + ( ------ ) - ( ------ ) + ... |
' m*10^e + 2j+1   m*10^e  \    m*10^e     m*10^e       m*10^e        /
'
' thus reducing the problem to multiplies and divides by small numbers
' (m and 2j+1), and shifts by a (possibly huge) exponent. It still isn't
' blazing fast, but it works nicely.

num = 6 * j + 1: remainder = 0
sourceword = firstword
destword = firstword + exponent%
term2(destword - 1) = 0

FOR i = destword TO digits
   dividend = remainder * 10000 + num * term1(sourceword)
   temp = dividend \ mantissa%
   remainder = dividend - temp * mantissa%           'create temporary term
   sourceword = sourceword + 1
   term2(i) = temp                           'term2 <- term1 * num / m*10^e
NEXT i

FOR i = destword TO digits
   temp = term2(i)
   IF temp >= 10000 OR temp <= -10000 THEN
      quotient = temp \ 10000                             'release carries 
      term2(i - 1) = term2(i - 1) + quotient        
      term1(i - 1) = term1(i - 1) - quotient
      temp = temp - 10000 * quotient
   END IF
   term1(i) = term1(i) - temp                       'term1 <- term1 - term2
   term2(i) = temp
NEXT i
IF term2(destword - 1) <> 0 THEN destword = destword - 1

num = 2 * j + 1
sourceword = destword
destword = destword + exponent%

DO UNTIL destword > digits
   term3(destword - 1) = 0: remainder = 0
  
   FOR i = destword TO digits                'term3 <- term2 * num / m*10^e 
      temp = term2(sourceword)                   
      dividend = remainder * 10000 + num * temp
      temp = dividend \ mantissa%
      remainder = dividend - temp * mantissa%
      IF temp >= 10000 OR temp <= -10000 THEN         'release carries
         quotient = temp \ 10000
         term3(i - 1) = term3(i - 1) + quotient
         term1(i - 1) = term1(i - 1) + quotient
         temp = temp - 10000 * quotient
      END IF                                    'term1 <- term1 + term3
      term1(i) = term1(i) + temp
      term3(i) = temp
      sourceword = sourceword + 1
   NEXT i
   IF term3(destword - 1) <> 0 THEN destword = destword - 1

   sourceword = destword: remainder = 0
   destword = destword + exponent%
   IF destword <= digits THEN term2(destword - 1) = 0
  
   FOR i = destword TO digits                 'term2 <- term3 * num / m*10^e  
      temp = term3(sourceword)
      dividend = remainder * 10000 + num * temp
      temp = dividend \ mantissa%
      remainder = dividend - temp * mantissa%
      IF temp >= 10000 OR temp <= -10000 THEN              'release carries
         quotient = temp \ 10000
         term2(i - 1) = term2(i - 1) + quotient
         term1(i - 1) = term1(i - 1) - quotient
         temp = temp - 10000 * quotient
      END IF
      term1(i) = term1(i) - temp                    'term1 <- term1 + term2
      term2(i) = temp
      sourceword = sourceword + 1
   NEXT i
   IF term2(destword - 1) <> 0 THEN destword = destword - 1

   sourceword = destword
   destword = destword + exponent%

LOOP

END SUB