Posts: 1,000
Threads: 50
Joined: May 2022
Reputation:
27
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
Posts: 20
Threads: 2
Joined: Apr 2022
Reputation:
3
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"
Posts: 1,000
Threads: 50
Joined: May 2022
Reputation:
27
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! 
Just a small correction:
Code: (Select All) $Console:Only
DefInt C-Q
DefLng Z
fmax = 5020
. . .
Posts: 3,984
Threads: 178
Joined: Apr 2022
Reputation:
222
Only 10 more days...
b = b + ...
Posts: 80
Threads: 7
Joined: Jan 2023
Reputation:
10
Here's my dartboard solution 
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 
RokCoder - dabbling in QB64pe for fun
Posts: 46
Threads: 7
Joined: Apr 2022
Reputation:
2
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
Posts: 3,984
Threads: 178
Joined: Apr 2022
Reputation:
222
03-06-2023, 06:18 PM
(This post was last modified: 03-06-2023, 07:56 PM by bplus.)
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
b = b + ...
Posts: 80
Threads: 7
Joined: Jan 2023
Reputation:
10
(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!
RokCoder - dabbling in QB64pe for fun
Posts: 200
Threads: 5
Joined: Apr 2022
Reputation:
22
03-07-2023, 12:34 AM
(This post was last modified: 03-07-2023, 12:55 AM by JRace.)
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.
|
