Code: Select all
;Many years ago when I used Quickbasic, I came across a program to calc pi. Of course this has no practical
;purpose, but just to fight winter boredom I modified the program (see comments included below) to run
;under Purebasic, now can use 10 digits per word instead of 4. Also added a bit of assembly to avoid calculating
;the remainder, this only makes a slight difference since most of time is taken by dividing.
;Using threads sped up by a factor of 3. Note calc atan(1/5) takes about twice as long as atan(1/239).
;Programmed by Jeff Wyatt, Highlands Ranch, Colorado (1/2019).
;On my Athelon II X4 635 the program can calculate 100,000 digits in about 8 secs with debugging off.
; >>>>>>>>>>>>>>>>>>>>>>>COMPILE With THREAD SAFE BOX CHECKED<<<<<<<<<<<<<<<<<<<<<<<<<<<<
;Based upon Machin's (1706) formula: pi/4=arctan(1)=4*arctan(1/5)-arctan(1/239), he did 100 digits by hand
;which would take about 8 usec with this program. Time is prop. to digits^2 so need to use another method
;for very large number of digits. Good way to test your computer's integer arithmatic.
;See https://www.angio.net/pi/ for further info and searchable listing of pi.
;I used it to check that my program worked correctly for the first 2,000,000 digits.
;Previous comments from authors 20 years ago
;Program To calculate pi, version 4.8
;A major rewrite of version 4.2, this uses only two arrays instead of
;three, and includes a host of speedups based on a similar C program.
;A sampler: all the carries are reserved until the end, the divide and
;add routines are combined, two terms are added at a time, and the number
;of function calls is minimized. It's a big change for a small gain, since
;the compiled version requires 28.6 seconds for 5000 digits on my 486 66 MHz
;computer, a 10 gain over version 4.2; like before, it's capable of about
;150,000 digits of pi.
;
;This program has come a long way from version 1.0; thanks are due to
;Larry Shultis, Randall Williams, Bob Farrington and Adrian Umpleby.
;One final note for speed freaks: this program will run about 6 times faster
;if written in C using an optimizing compiler. Likewise, if you can figure
;out a way to do integer division and use both the quotient and remainder,
;this program can easily be sped up by 25. -jasonp@isr.umd.edu
Declare atan239(dummy): Declare atan51(dummy): Declare atan52(dummy)
Declare PrintOut(): Declare WriteOut(filename$)
Global words, big,t239.f,t52.f,btime.f,ctime.f,Digits,DigitsperWord
DigitsperWord=10 ;change if you wish, >10 print out will be incorrect calc may overflow,
big=1 ;<10 print will have lots of extraneous zeros
For i=1 To DigitsperWord ;1 followed by 10 zeros, note if you use 14 digits/word
big=big*10 ;arctan(1/5) will overflow about digit 64,760
Next ;when 1.4*digits*big~denom*big >9e18
OpenConsole()
Print("how many digits ")
Digits=Val(Input())
words = Digits / DigitsperWord+2
Global Dim sum(words + 2),Dim sum1(words+2), Dim sum2(words+2)
CPUs=CountCPUs()
start = ElapsedMilliseconds()
;--------------- -4*atan(1/239)
If CPUs>1 ;note procedure must have one parameter for threading
a239=CreateThread(@atan239(),1)
Else
atan239(dummy)
EndIf
;------------ 16*atan(1/5) Do in two parts so can be threaded
atime.f=(ElapsedMilliseconds()-start)/1000.0
If CPUs>2
a51=CreateThread(@atan51(),1)
Else
atan51(dummy)
EndIf
atan52(dummy)
If IsThread(a51):WaitThread(a51):EndIf
btime.f=ElapsedMilliseconds()/1000.0-atime
If IsThread(a239):WaitThread(a239):EndIf
sum1(2)=sum1(2)-big/5 ;remove since calculate this term twice, intially first term
;is 16/5=3.2, so sum1(2) and sum(2) will contain .2 which is
;big/5, now can sum the three sums
For x=2 To words
sum(x)=sum(x)+sum1(x)+sum2(x)
Next
For x = words To 2 Step -1 ;finish up
If sum(x) < 0 ;release borrows
quotient = sum(x) / big
sum(x) = sum(x) - (quotient - 1) * big
sum(x - 1) = sum(x - 1) + quotient - 1
EndIf
If sum(x) >= big ;and carries
quotient = sum(x) / big
sum(x) = sum(x) - quotient * big
sum(x - 1) = sum(x - 1) + quotient
EndIf
Next
ctime.f=(ElapsedMilliseconds()-start)/1000.0
PrintOut()
filename$=GetCurrentDirectory()+"PiData.txt"
PrintN("Will write a file to "+filename$)
Print(" Save to file (y/n)")
If UCase(Input())="Y"
WriteOut(filename$)
EndIf
End
;------------------------------------------------------------------
Procedure atan239(dummy) ;arctan(x) = x-x^3/3+x^5/5-x^7 ....
t239=ElapsedMilliseconds()
Dim term(words +2)
remainder=4
For x = 2 To words
dividend = remainder * big ;crunch out 1st term
term(x) = dividend / 239
remainder = dividend - term(x) * 239
sum2(x) = sum2(x) - term(x) ;subtract as we want -atan(1/239)
Next x
denom=3:firstword=2
Repeat ;do two more terms, add first, subtract second
remainder1=0
remainder2=0
For x = firstword To words
temp = term(x)
dividend = remainder1 * big + temp
temp = dividend / 57121
remainder1 = dividend - temp * 57121
term(x) = temp
dividend = remainder2 * big + temp
temp = dividend / denom
remainder2 = dividend - temp * denom
sum2(x) = sum2(x) + temp
Next
If term(firstword) = 0 : firstword = firstword + 1:EndIf
denom = denom + 2
remainder1 = 0: remainder2 = 0
For x = firstword To words
temp = term(x)
dividend = remainder1 * big + temp
temp = dividend / 57121
remainder1 = dividend - temp * 57121
term(x) = temp
dividend = remainder2 * big + temp
temp = dividend / denom
remainder2 = dividend - temp * denom
sum2(x) = sum2(x) - temp
Next x
If term(firstword) = 0 : firstword = firstword + 1:EndIf
denom = denom + 2
Until firstword >=words
t239=ElapsedMilliseconds()-t239
EndProcedure
;-------------------------------------------------------------------
Procedure atan51(dummy) ; atan(1/5) one half of calc terms 1/5, 1/9, 1/13 ......
Dim term(words +2)
denom = 5: firstword = 1: lastword =3
sum(1) = 3: term(1) = 3: sum(2) = big/5: term(2) = sum(2)
Repeat
remainder1=0
remainder2=0
For x = firstword To lastword + 1
temp = term(x)
dividend = remainder1 * big + temp
temp = dividend / 625 ; 625 = 5^4
;remainder1 = dividend - temp * 625
EnableASM ;Purebasic makes doing this so easy, couldn't resist as remainder
mov remainder1, rdx ;already is in the rdx register
DisableASM
term(x) = temp
dividend = remainder2 * big + temp
temp = dividend / denom
;remainder2 = dividend - temp * denom
EnableASM
mov remainder2, rdx
DisableASM
sum(x) = sum(x) +temp
Next
For x = lastword + 2 To words
dividend = remainder2 * big
temp = dividend / denom
;remainder2 = dividend - temp * denom
EnableASM
mov remainder2, rdx
DisableASM
sum(x) = sum(x) + temp
Next
If term(lastword + 1) > 0 And lastword < words : lastword = lastword + 1:EndIf
If term(firstword) = 0 : firstword = firstword + 1:EndIf
denom = denom + 4
Until firstword >= words
EndProcedure
Procedure atan52(dummy) ; atan(1/5) other half of calc terms 1/3, 1/7, 1/11 ....
t52=ElapsedMilliseconds()
Dim term(words +2)
denom = 3: firstword = 1: lastword =3
sum1(1) = 3: term(1) = 3: sum1(2) = big/5: term(2) = sum1(2)
remainder1=0
remainder2=0
For x = firstword To lastword + 1 ;term 1/3
temp = term(x)
dividend = remainder1 * big + temp
temp = dividend / 25
;remainder1 = dividend - temp * 625
EnableASM
mov remainder1, rdx
DisableASM
term(x) = temp
dividend = remainder2 * big + temp
temp = dividend / denom
;remainder2 = dividend - temp * denom
EnableASM
mov remainder2, rdx
DisableASM
sum1(x) = sum1(x) - temp
Next
For x = lastword + 2 To words
dividend = remainder2 * big
temp = dividend / denom
;remainder2 = dividend - temp * denom
EnableASM
mov remainder2, rdx
DisableASM
sum1(x) = sum1(x) - temp
Next
If term(lastword + 1) > 0 And lastword < words : lastword = lastword + 1:EndIf
If term(firstword) = 0 : firstword = firstword + 1:EndIf
denom = denom + 4
Repeat ;do the rest
remainder1 = 0: remainder2 = 0
For x = firstword To lastword + 1
temp = term(x)
dividend = remainder1 * big + temp
temp = dividend / 625
;remainder1 = dividend - temp * 625
EnableASM
mov remainder1, rdx
DisableASM
term(x) = temp
dividend = remainder2 * big + temp ;eventually this will overflow >9e18
temp = dividend / denom
;remainder2 = dividend - temp * denom
EnableASM
mov remainder2, rdx
DisableASM
sum1(x) = sum1(x) - temp
Next x
For x = lastword + 2 To words
dividend = remainder2 * big
temp = dividend / denom
;remainder2 = dividend - temp * denom
EnableASM
mov remainder2, rdx
DisableASM
sum1(x) = sum1(x) -temp
Next x
If term(lastword + 1) > 0 And lastword < words: lastword = lastword + 1:EndIf
If term(firstword) = 0 : firstword = firstword + 1:EndIf
denom = denom + 4
Until firstword >= words
t52=ElapsedMilliseconds()
EndProcedure
;------------------------------------------------------------------
Procedure PrintOut()
PrintN("")
p$="pi = 3."
i=2
Repeat
For j=i To i+4
If j>words
p$=p$+Space(11)
Else
p$=p$+RSet(Str(sum(j)),10,"0")+" "
EndIf
Next
p$=p$+RSet(Str(10*i+30),7)
PrintN(p$)
p$=" "
i=i+5
Until i>words
PrintN("")
PrintN(" arctan(1/5) computation time: "+StrF(btime,2)+" seconds")
PrintN(" arctan(1/239) computation time: "+StrF(t239/1000,2)+" seconds")
PrintN(" total computation time: "+StrF(ctime,2)+" seconds")
EndProcedure
Procedure WriteOut(filename$)
CreateFile(1,filename$)
WriteStringN(1,"pi data "+Str(digits)+" Digits "+Str(DigitsperWord)+" Digits/Word pi = 3.")
i=2
Repeat
p$=""
For j=i To i+4
If j>words
p$=p$+Space(11)
Else
p$=p$+RSet(Str(sum(j)),10,"0")+" "
EndIf
Next
p$=p$+RSet(Str(10*(i)+30),7)
WriteStringN(1,p$)
i=i+5
Until i>words
WriteStringN(1,"")
WriteStringN(1," arctan(1/5) computation time: "+StrF(btime,2)+" seconds")
WriteStringN(1," arctan(1/239) computation time: "+StrF(t239/1000.0,2)+" seconds")
WriteStringN(1," total computation time: "+StrF(ctime,2)+" seconds")
CloseFile(1)
EndProcedure