Program to calculate pi

[Kernelpanic]

Pete would love the program. For others, one or the other might be of interest.

Code: (Select All)

'===========================================================================
' 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

' Mit Locate und mit Color Ausgabe angepasst - 26. Feb. 2023
'===========================================================================
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

'----------- Intro Screen by (c) Marc Antoni, Oct. 2, 2000 -----------------
Color 7, 0
Cls
Locate 10: Print "  Pi-Berechnung nach Euler (1707 - 1783)"
Locate 12: Print "    (Pi^2)/8 = 1/1^2 + 1/3^2 + 1/5^2 + 1/7^2 + ..."
Locate 20: Print "              Programming by Jason Stratos Papadopoulos"
Locate 24: Print "              ... weiter mit beliebiger Taste"
Do: Loop While InKey$ = ""

'----------- End of Intro Screen -------------------------------------------
Cls
Locate 2, 2
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-------
Locate CsrLin + 1, 2
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)
Do: Loop While InKey$ = ""
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:
  Locate CsrLin + 1, 2
  Color 4, 0
  Print "pi=3."

  'Wieder zuruecksetzen
  Color 7, 0
  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:
  Locate CsrLin + 2, 2
  Print Time$

End Sub