QB64 Phoenix Edition
Program to calculate pi - Printable Version

+- QB64 Phoenix Edition (https://staging.qb64phoenix.com)
+-- Forum: QB64 Rising (https://staging.qb64phoenix.com/forumdisplay.php?fid=1)
+--- Forum: Code and Stuff (https://staging.qb64phoenix.com/forumdisplay.php?fid=3)
+---- Forum: Programs (https://staging.qb64phoenix.com/forumdisplay.php?fid=7)
+---- Thread: Program to calculate pi (/showthread.php?tid=1510)

Program to calculate pi - Kernelpanic - 02-26-2023

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


RE: Program to calculate pi - david_uwi - 03-04-2023

I wrote one in QB45 it will work until long integers overflow about 5000 digits. It uses Machin's formula. It seems marginally slower that the one you posted.
Code: (Select All)
DefInt C-Q
DefLng Z
fmax = 5020
Dim i(fmax) As Integer, k(fmax) As Integer
Dim ii(fmax) As Integer, kk(fmax) As Integer
tt = Timer
fl = -1: ii(1) = 10: l = 1: m = -1
'calculating arctan(1/5) to fmax d.p.s
While l < fmax - 20
    m = m + 2
    fl = fl * -1
    If m = 1 Then q = 5 Else q = 25
    ir = 0
    For j = l To fmax
        ii(j) = ii(j) + ir * 10
        ir = ii(j) Mod q
        ii(j) = ii(j) \ q
    2 Next j
    zir = 0
    For l = l To fmax
        If ii(l) <> 0 GoTo 5
    Next l
    l = l - 20
    5 If l < 40 Then l = 1
    If m > 3200 Then
        For j = l To fmax
            i(j) = ii(j)
            z = i(j) + zir * 10
            zir = z Mod m
            i(j) = z \ m
        Next j
    End If
    If m < 3201 Then
        ir = 0
        For j = l To fmax
            i(j) = ii(j)
            i(j) = i(j) + ir * 10
            ir = i(j) Mod m
            i(j) = i(j) \ m
        Next j
    End If
    cf = 0
    'add routine
    If fl = 1 Then
        For j = fmax To l Step -1
            If i(j) > 0 Or cf > 0 Then
                k(j) = k(j) + i(j) + cf
                cf = 0
                If k(j) > 9 Then
                    cf = k(j) \ 10
                    k(j) = k(j) Mod 10
                End If
            End If
        Next j
    End If
    If fl = -1 Then
        'subtract routine
        For j = fmax To l Step -1
            k(j) = k(j) - i(j) - cf
            cf = 0
            If k(j) < 0 Then
                cf = 1
                k(j) = 10 + k(j)
            End If
        Next j
    End If
Wend
For j = 1 To fmax
    kk(j) = k(j)
    k(j) = 0: i(j) = 0: ii(j) = 0
Next j
'now let's calculate arctan(1/239)
fl = -1: ii(1) = 10: l = 1: m = -1
While l < fmax - 20
    m = m + 2
    fl = fl * -1
    If m = 1 Then q = 1 Else q = 2
    For n = 1 To q
        ir = 0
        For j = l To fmax
            ii(j) = ii(j) + ir * 10
            ir = ii(j) Mod 239
            ii(j) = ii(j) \ 239
        Next j
    Next n
    For l = l To fmax
        If ii(l) <> 0 GoTo 7
    Next l
    7 l = l - 20
    If l < 40 Then l = 1
    ir = 0
    For j = l To fmax
        i(j) = ii(j)
        i(j) = i(j) + ir * 10
        ir = i(j) Mod m
        i(j) = i(j) \ m
    Next j
    cf = 0
    'add routine
    If fl = 1 Then
        For j = fmax To l Step -1
            k(j) = k(j) + i(j) + cf
            cf = 0
            If k(j) > 9 Then
                cf = k(j) \ 10
                k(j) = k(j) Mod 10
            End If
        Next j
    End If
    If fl = -1 Then
        'subtract routine
        For j = fmax To l Step -1
            k(j) = k(j) - i(j) - cf
            cf = 0
            If k(j) < 0 Then
                cf = 1
                k(j) = 10 + k(j)
            End If
        Next j
    End If
Wend
'multipling atn(1/5) by 16
For m = 1 To 4
    cf = 0
    For j = fmax To 1 Step -1
        kk(j) = kk(j) + kk(j) + cf
        cf = kk(j) \ 10
        kk(j) = kk(j) Mod 10
    Next j
Next m
'multipling atn(1/239) by 4
For m = 1 To 2
    cf = 0
    For j = fmax To 1 Step -1
        k(j) = k(j) + k(j) + cf
        cf = k(j) \ 10
        k(j) = k(j) Mod 10
    Next j
Next m
'16*atn(1/5)-4*atn(1/239) = pi
For j = fmax To 1 Step -1
    kk(j) = kk(j) - k(j) - cf
    cf = 0
    If kk(j) < 0 Then
        cf = 1
        kk(j) = 10 + kk(j)
    End If
Next j
Print ""
Print Using "&"; "3.";
l = -1
For j = 1 To fmax - 20
    l = l + 1
    If l = 50 Then Print "": Print "  ";
    If l = 50 Then l = 0
    Print Using "#"; kk(j);
Next j
Print: Print "time taken= "; Timer - tt; " sec"


RE: Program to calculate pi - Kernelpanic - 03-04-2023

Quote:@david_uwi - I wrote one in QB45 it will work until long integers overflow about 5000 digits. It uses Machin's formula. It seems marginally slower that the one you posted.

Good program! Pi won't let go of people!  Big Grin

Just a small correction:

Code: (Select All)
$Console:Only

DefInt C-Q
DefLng Z
fmax = 5020
. . .


RE: Program to calculate pi - bplus - 03-04-2023

Only 10 more days...

RE: Program to calculate pi - RokCoder - 03-06-2023

Here's my dartboard solution  Cool

Code: (Select All)
RANDOMIZE TIMER
CONST N& = 10000000
FOR i& = 1 TO N&
    x! = RND * 2 - 1
    y! = RND * 2 - 1
    h& = h& - (x! * x! + y! * y! < 1)
NEXT i&
PRINT h& / N& * 4

I was going to add graphics of darts being thrown at the dartboard but I can procrastinate well enough without going to those extremes!

And, yes, I know it's a rubbish solution but it's kinda fun Smile

RE: Program to calculate pi - BSpinoza - 03-06-2023

Long time ago... 
"Esimating of pi by using the Monte Carlo Method" with graphics...
Code: (Select All)
_TITLE " Estimating of PI using the Monte Carlo Method "
'written May 08, 2018 by Bruno Schaefer, Losheim am See, Germany
'
' A simple Monte Carlo method to compute the value of Pi uses some very basic geometric relationships.
' Consider a circle of radius R inscribed in a square with a side length of 2R.
' If R=1, them the area of the circle is Pi*R^2 = Pi. Now consider the upper right quadrant.
' The area of the quarter circle is (Pi·R2)/4 = Pi/4. The area of the quadrant is R^2 = 1.
' Consider a random point (x,y) in the quadrant. Note the probability it is below the curve is simply
' the area of the region under the curve compared to the total area of the quadrant.
' This is simply Pi/4.
' To compute Pi using this method, one needs to generate N random points (x,y) in the quadrant.
' Count the number of these points that are below the curve (x^2 + y^2 < R^2), a "hit",
' here this number is called H.
' Pi/4 can be approximated as the number of hits divided by the total number of points generated by
' noting that H/G approaches Pi/4 as N goes to infinity.
DIM mcpi AS _FLOAT, x AS _FLOAT, y AS _FLOAT
N = 10000000 '! number of points on the whole area
H = 0 ' ! number of points in the quarter of a circle
RANDOMIZE TIMER
SCREEN _NEWIMAGE(630, 690, 256)
CLS , 15
COLOR 0, 15
WINDOW (-40, -30)-(590, 640)
FOR y = 40 TO 590 STEP 55
    LINE (15, y)-(565, y), 7 '  helplines parallel to x
NEXT y
FOR x = 15 TO 565 STEP 55
    LINE (x, 40)-(x, 590), 7 '   helplines parallel to y
NEXT x
FOR y = 40 TO 590 STEP 55
    LINE (12, y)-(18, y), '  ticks y-axis
NEXT y
FOR x = 15 TO 565 STEP 55 '  Ticks x-axis
    LINE (x, 37)-(x, 43) ', 7
NEXT x
LINE (15, 40)-(565, 40) 'x-axis
LINE (15, 40)-(15, 590) 'y-axis
LINE (15, 590)-(565, 590) 'upper line
LINE (565, 40)-(565, 590) 'right line
LOCATE 4, 3: PRINT " 1.0"
LOCATE 11, 3: PRINT " 0.8"
LOCATE 18, 3: PRINT " 0.6"
LOCATE 25, 3: PRINT " 0.4"
LOCATE 32, 3: PRINT " 0.2"
LOCATE 39, 3: PRINT " 0.0"
LOCATE 40, 7: PRINT "0.0"
LOCATE 40, 20: PRINT "0.2"
LOCATE 40, 34: PRINT "0.4"
LOCATE 40, 48: PRINT "0.6"
LOCATE 40, 61: PRINT "0.8"
LOCATE 40, 75: PRINT "1.0"
LOCATE 1, 30: PRINT " " + CHR$(227) + " = " + STR$(_PI(1)) + " "
FOR i = 1 TO N
    vx = RND
    vy = RND
    xWert = 15 + vx * 550
    yWert = 40 + vy * 550
    IF (hypot(vx, vy) <= 1.0) THEN
        LET H = H + 1
        PSET (xWert, yWert), 12
    ELSE
        PSET (xWert, yWert), 9
    END IF
    mcpi = 4 * H / i
    LOCATE 2, 30: PRINT " N = " + STR$(i) + " "
    LOCATE 3, 30: PRINT " " + CHR$(227) + "(est.) = " + STR$(mcpi) + "             "
NEXT i
END

FUNCTION hypot (x, y)
    hypot = SQR(x ^ 2 + y ^ 2)
END FUNCTION


RE: Program to calculate pi - bplus - 03-06-2023

FYI QB64 has _Hypot(x, y) function!
https://qb64phoenix.com/qb64wiki/index.php/HYPOT

For 2 points (x1, y1) and (x2, y2)
dist = _Hypot(x1 - x2, y1 - y2)

Edit: Correction there should be a minus sign between y1 and y2 not a comma! apologies

RE: Program to calculate pi - RokCoder - 03-06-2023

(03-06-2023, 06:18 PM)bplus Wrote: FYI QB64 has _Hypot(x, y) function!
https://qb64phoenix.com/qb64wiki/index.php/HYPOT

For 2 points (x1, y1) and (x2, y2)
dist = _Hypot(x1 - x2, y1, y2)

I didn't know that either. I imagine there are loads of utility functions added since the original BASIC that I should be using. In this instance it's surely quicker to stick with x.x + y.y and no square root or hypot function at all. Nice to see a less lazy version than my quick one though! I didn't realise this technique actually had a name!

RE: Program to calculate pi - JRace - 03-07-2023

I KNEW I had an old Pi around here somewhere.  I could smell it.
This was originally written for the TRS-80.  I recently stumbled over it on one of my old disk images.

No extended-precision tricks here.  This program just closes in on Pi to the limits of double-precision accuracy, which was slow enough on 8-bit machines.


This first version uses Basic's built-in SQR function:
Code: (Select All)
1 'Originally written for TRS-80 Model III Basic; spaces added for later Basics.
10 DefDbl A-Z
20 PI = 2
30 S = Sqr(S + 2): P = PI: PI = 2 * P / S
35 Print P; Tab(32); PI
40 If P <> PI Then 30
50 Print "PI ="; PI
Okay, THAT version of the program is good enough for all of the PC Basics that I've tried it on (QB64PE & FB under Windoze, and QBasic, GWBasic, & Turbo Basic under DOS).


Buuut: the TRS-80's SQR only calculated to single-precision accuracy, so here's a workaround for that:
Code: (Select All)
1 'Originally written for TRS-80 Model III Basic; spaces added for later Basics.
10 DefDbl A-Z
20 PI = 2: E = .000000000000001
30 N = R + 2: GoSub 100: PP = PI: PI = PI * 2 / R
35 Print PP; Tab(32); PI
40 If Abs(PP - PI) > E Then 30
50 Print "PI ="; PI
60 End

100 A = 0: B = N
120 S = (B - A) / 10
130 If Abs(A * B - N) < E Then R = (A + B) / 2: Return
140 If Abs(A * A - N) < E Then R = A: Return
150 If Abs(B * B - N) < E Then R = B: Return
160 I = A
170 If N < I * I Then B = I: A = I - S: GoTo 120
180 I = I + S: If I > B Then B = I + S: GoTo 120
190 GoTo 170
300 'Lines 100-190 form a more accurate square root subroutine.

301 'The custom square-root subroutine was for the TRS-80, to improve accuracy.

302 'It is not needed for QB64PE, FB, QBasic, GWBasic, or Turbo Basic,
303 'where it actually ruins the accuracy of the last digit of Pi.
304 'The Basics listed above can use the built-in SQR function and nail
305 'every single digit under double-precision.
If you actually run this version on any of the Basics listed above, you will see that this home-brewed SQR function falls a bit short and the last calculated digit of Pi is wrong.  ¯\_(ツ)_/¯


I haven't tried compiling these programs in BasCom, nor shall I.  After spending a couple hours one day tracking down info on the goofy long-forgotten command-line switches for that compiler & linker, I'm sick of early MS dev tools now.