Download text file
#if 0
---------------------------- PowerBASIC/cc v2.0
---| DASoft |------------------------------------------
---------------------------- Code DATE: 2000-04-08
| FILE NAME PrimeGEN.bas | by
---------------------------- Don Schullian, Jr.
This code is released into the Public Domain
---------------------------------------------------------
No guarentee as to the viability, accuracy, or safty of
use of this code is implied, warrented, or guarenteed
---------------------------------------------------------
Use at your own risk!
---------------------------------------------------------
CONTACT AUTHOR AT d83@DASoftVSS.com
-------------------------------------------------------------------------
Hi,
First of all, I'm NOT a mathematician so bare with me here.
The code below contains ideas from, at least, three different sources plus one
or two of my own. The code is also optimized for speed over size as each clock
tick in the main loop could cost quadrillions of ticks throughout the run of
the program. And, lastly, when/if the user bails out or the program hits the
maximum number requested the next number in line to be checked will end in a
seven (7). This was done to facilitate restarting the program where it left off.
The overall theory of this prime number generator is to test only 40% of the
existing numbers. All even numbers (except 2) are NOT primes and all numbers
ending in 5 (except 5) are not primes. Hence only those numbers ending in 1, 3,
7 or 9 can qualify. Testing only odd numbers is quite common but I've added the
double loop that also skips those numbers divisible by 5.
Next, 1, 2, 3, and 5 are loaded and testing starts with 7. The generated file
then holds only 1/2 the difference between one prime and the following prime. This
assumes that no two prime numbers are more than 510 apart. (as both primes are odd,
the difference will always be even) So far, with all my testing this condition has
not even been neared and as I've gone past 40M primes I'd guess that it is safe.
There is an error test for this condition and the program will end if it is exceeded.
There are 2 constants (below) that manage the buffer sizes:
%BufSize determines the size of the 'testing primes'. 4k was sufficient to test
the first 40M primes, with 100 values not used.
%HoldSize determines how many new primes are found before the buffer is flushed
to the file. I'd suggest that you keep this value an increment of 4k to maximize
the speed which it will transfer the data.
If the highest prime you require will fit into a LONG INTEGER then change all
the QUADs to LONGs. (40Mth prime only in the range of 800M) This holds true also
for DWORDs (4B). Either of these simple changes will speed the program up by
quite a bit.
If you plan on running this program all night then make sure that you turn off
ALL screen savers, power savers, etc. or you'll end up with a program sitting around
doing nothing!
The sub PrimeView is incomplete so you can fill in the pertinent print command(s)
that you will require. Remember, if you go to disk then make sure you've got some
space available! 40M primes printed out in text will require over 760M of file space
and reams of paper!
#endif
'--------------------------------------------------'
' PRIME NUMBER GENERATOR '
'--------------------------------------------------'
%BufSize = 4096 ' 4k
%HoldSize = 32768 ' 32k
%DisplayReports = -1 '
'--------------------------------------------------'
' '
'--------------------------------------------------'
DECLARE SUB PrimeGen (BYVAL F AS STRING,BYVAL N AS QUAD,BYVAL S AS LONG)
DECLARE SUB PrimeView (BYVAL F AS STRING)
'--------------------------------------------------'
' '
'--------------------------------------------------'
FUNCTION PBmain () '
'
DIM MaxNbr AS LOCAL QUAD '
'
CONSOLE SCREEN 25, 80 '
CURSOR OFF '
'
MaxNbr = &h7FFFFFFFFFFFFFFF ' 4+ quintillian is max
'
PrimeGen "TEST.DAT", MaxNbr, 0 '
'
END FUNCTION '
'--------------------------------------------------'
' '
'--------------------------------------------------'
SUB PrimeGen ( BYVAL FileSpec AS STRING, _
BYVAL MaxNbr AS QUAD , _
BYVAL Restart AS LONG )
DIM B AS LOCAL LONG ' loop counter
DIM Bdata AS LOCAL STRING * %BufSize ' temporary data buffer
DIM Buf(%BufSize) AS LOCAL DWORD ' working buffer for primes
DIM Boff AS LOCAL LONG ' working buffer loop counter
DIM B_ptr AS LOCAL BYTE PTR ' pointer to temp buffer
DIM Count AS LOCAL QUAD ' found/saved prime counter
DIM D AS LOCAL STRING ' display string
DIM Etime AS LOCAL SINGLE ' elapsed time storage
DIM Hold AS LOCAL STRING * %HoldSize ' found prime # storage
DIM H_ptr AS LOCAL BYTE PTR ' pointer to above
DIM Hlast AS LOCAL LONG ' %HoldSize -1
DIM Hoff AS LOCAL LONG ' offset for H_ptr
DIM LastPrime AS LOCAL QUAD ' previously found prime number
DIM N AS LOCAL QUAD ' working number loop counter
DIM Nbr AS LOCAL QUAD ' sub loop working number
DIM S AS LOCAL DWORD ' sqr value of working number
DIM Start AS LOCAL QUAD ' 1st number to start search (7 is default)
DIM Ok AS LOCAL LONG ' prime found flag
DIM X AS LOCAL LONG ' junk variable for testing, etc.
'
B_ptr = VARPTR(Bdata) ' working buffer pointer
H_ptr = VARPTR(Hold) ' holding buffer for primes
Start = 7 '
LastPrime = 5 ' last 'found' prime number
'
IF Restart = 0 THEN '
IF LEN(DIR$(FileSpec)) > 0 THEN KILL FileSpec ' clear the existing file
OPEN FileSpec FOR BINARY AS #1 BASE=0 ' open new file
ELSE '
OPEN FileSpec FOR BINARY AS #1 BASE=0 ' open existing file
GET$ #1, 4, D ' read signature characters
N = LOF(1) - 4 ' get length of file
IF (N < 0 ) OR ( D <> CHR$(0,0,0,1) ) THEN ' test if our file or not
CLOSE #1 ' no... it isn't
PRINT FileSpec; "isn't one of my files" ' print warning message
WAITKEY$ ' wait for a keypress
EXIT SUB ' exit program
END IF '
WHILE N > 0 ' while data in file
Nbr = MIN(N,%HoldSize) ' compute next chunk size
N = N - Nbr ' decrease bytes left in file
GET$ #1, Nbr, Hold ' read chunk into buffer
FOR B = 0 TO Nbr-1 ' read each byte of buffer
LastPrime = LastPrime + ( 2 * @H_ptr[B] ) ' compute next prime number
NEXT '
WEND '
Start = LastPrime ' the next starting prime will end in a 7
DO ' so we search it out
Start = Start + 2 '
LOOP UNTIL (Start MOD 10) <> 7 '
END IF '
'
Buf(1) = 3 ' we'll start here
Boff = 1 ' working buffer offset
Ok = -1 ' prime found flag
'
FOR N = 7 TO &h7FFFFFFFFFFFFFFF STEP 10 ' create working primes
FOR Nbr = N TO (N + 6) STEP 2 ' 7, 9, 11, 13 skip numbers divisible by 5
S = SQR(Nbr) ' get sqr of N s
FOR B = 1 TO Boff ' check previous primes against N
IF Buf(B) > S THEN EXIT FOR ' this prime > SQR(N) so we're done
IF (Nbr MOD Buf(B)) = 0 THEN ' this N divides equally by Prime(B)
Ok = 0 ' turn off 'found' flag
EXIT FOR ' bail out of testing loop
END IF '
NEXT '
IF Ok = 0 THEN '
Ok = -1 ' This ain't a prime
ELSE '
INCR Boff ' set last found prime flag
Buf(Boff) = Nbr ' store found value in array
IF Boff = %BufSize THEN EXIT, EXIT ' the buffer is full
END IF '
NEXT '
NEXT '
'
Count = SEEK(1) '
Hold = CHR$(0,0,0,1) ' 1,2,3,5
Hoff = 3 ' hold buffer offset
Hlast = %HoldSize - 4 ' last hold buffer offset
Ok = -1 ' prime found flag
Boff = 3 ' working prime buffer offset for file input
'
Etime = TIMER ' start the timer
'
FOR N = Start TO MaxNbr STEP 10 ' This can be one HUMUMGUS! number
FOR Nbr = N TO (N + 6) STEP 2 ' 7, 9, 11, 13 skip numbers divisible by 5
S = SQR(Nbr) ' get sqr of N s
DO '
FOR B = 1 TO %BufSize ' check previous primes against N
IF Buf(B) > S THEN EXIT, EXIT ' this prime > SQR(N) so we're done
IF (Nbr MOD Buf(B)) = 0 THEN ' this N divides equally by Prime(B)
Ok = 0 ' turn off 'found' flag
EXIT, EXIT ' bail out of testing loop
END IF '
NEXT '
Boff = Boff + %BufSize ' we've run out of buffered primes
Buf(0) = Buf(%BufSize) ' set last known prime
GET #1, Boff, Bdata ' load in data
FOR B = 0 TO %BufSize -1 ' compute new prime set
Buf(B+1) = Buf(B) + (@B_ptr[B] * 2) '
NEXT '
LOOP '
'
IF Boff > 3 THEN ' if we've loaded a new set of primes
Buf(1) = 3 ' back up to 1st set of primes
Buf(2) = 7 ' notice that 5 is skipped
GET #1, 4, Bdata ' load 1st %BufSize primes
FOR B = 2 TO %BufSize -1 ' reset buffer values
Buf(B+1) = Buf(B) + (@B_ptr[B] * 2) '
NEXT '
SEEK #1, Count ' replace file pointer to end
Boff = 3 ' reset flag
END IF '
'
IF Ok = 0 THEN ' no, this number isn't a prime
Ok = -1 '
ITERATE '
END IF '
'
Ok = (Nbr - LastPrime) ' difference between last 2 primes
IF Ok > 510 THEN ' oops! Houston, we have a problem!
PRINT "Difference Overflow at"; Nbr '
EXIT, EXIT '
END IF '
INCR Hoff '
SHIFT RIGHT Ok, 1 ' \ difference by 2
LastPrime = Nbr ' set new last prime found
@H_ptr[Hoff] = Ok ' store difference
NEXT '
'
IF Hoff => Hlast THEN ' if hold buffer is full, write to disk
INCR Hoff ' this is outside the Nbr loop to make
PUT$ #1, LEFT$(Hold, Hoff) ' restarting easier to compute
Count = Count + Hoff '
Hoff = -1 '
IF INSTAT THEN ' if keyboard action
IF INKEY$ = CHR$(27) THEN EXIT FOR ' if such action was
END IF '
#if %DisplayReports '
LOCATE 1, 1 '
PRINT FORMAT$(Count, ",# Primes in"); '
PRINT FORMAT$(Nbr ," ,# Nmbrs at "); '
PRINT LEFT$(TIME$,5); '
#endif '
END IF '
NEXT '
'
IF Hoff > -1 THEN '
INCR Hoff ' clear out remaining primes in Hold buffer
PUT$ #1, LEFT$(Hold,Hoff) '
Count = SEEK(1) ' prime count
END IF '
SETEOF 1 '
CLOSE #1 '
'------------------------------------------------'
'------------------------------------------------'END OF FUNCTIONAL PORTION OF CODE
'------------------------------------------------'
Etime = TIMER - Etime '
IF Etime < 0 THEN Etime = Etime + 86400 '
LOCATE 1, 1 '
PRINT FORMAT$(Count, ",# Primes in"); '
PRINT FORMAT$(Nbr ," ,# Numbers at "); '
PRINT LEFT$(TIME$,5); '
PRINT '
PRINT FORMAT$(Etime,"ELAPSED TIME ,#.####sec") '
INPUT FLUSH '
WAITKEY$ '
'
END SUB '
'
'-------------------------------------------------------------------------------
'
SUB PrimeView ( BYVAL FileSpec AS STRING )
DIM Bytes AS LOCAL QUAD '
DIM D AS LOCAL STRING ' display string
DIM Hold AS LOCAL STRING * %HoldSize ' found prime # stroage
DIM H_ptr AS LOCAL BYTE PTR '
DIM Hlast AS LOCAL LONG '
DIM Hoff AS LOCAL LONG '
DIM LastPrime AS LOCAL QUAD '
'
FileSpec = "TEST2.DAT" '
H_ptr = VARPTR(Hold) '
'
IF LEN(DIR$(FileSpec)) = 0 THEN '
PRINT FileSpec; " has not been found." '
WAITKEY$ '
EXIT SUB '
END IF '
'
OPEN FileSpec$ FOR BINARY AS #1 '
GET$ #1, 4, D '
IF D <> CHR$(0,0,0,1) THEN '
CLOSE #1 '
PRINT FileSpec; " is not our file." '
WAITKEY$ '
EXIT SUB '
END IF '
SEEK #1, 4 '
Bytes = LOF(1) - 3 '
LastPrime = 3 '
WHILE Bytes > 0 '
Hlast = MIN(%HoldSize,Bytes) '
GET$ #1, Hlast, Hold '
FOR Hoff = 0 TO Hlast -1 '
LastPrime = LastPrime + ( @H_ptr[Hoff] * 2 ) '
' put print statement here '
NEXT '
Bytes = Bytes - Hlast '
WEND '
PRINT '
PRINT FORMAT$(LOF(1), "COUNT: ,#") '
PRINT FORMAT$(LastPrime," LAST: ,#") '
CLOSE #1 '
WAITKEY$ '
'
END SUB