[futurebasic] Re: [FB] Working Prime Numbers Program

Message: < previous - next > : Reply : Subscribe : Cleanse
Home   : March 2000 : Group Archive : Group : All Groups

From: Jay Reeve <jktr@...>
Date: Thu, 23 Mar 00 22:53:33 -0600
Hi all,

Here is my next (and final) entry in the Prime Numbers sweepstakes. This 
version is twice as fast and requires 30% less memory than the previous 
one--and better documented. It won't compete with Mathematica, but I 
think it's a neat bit of code anyway.

 =J= a  y

// Prime finder            //
// Adapted from Phil Yates //
// by Jay Reeve            //
//Requires FBII runtime

_MaxNumber = 450'20000000

//Set _MaxNumber to whatever you like. 
//Max Number allowable is 2,147,483,547
//Memory required is approx 41 k per million integers.

_FlagBytes = (_MaxNumber) / 24 + 1'bytes needed for flags
_max3rd    = _MaxNumber/3

dim gPrimeFlags;_FlagBytes'Allocate _FlagBytes bytes for flags
dim gRunTime!'for timing
dim column,row'for printing
end globals

WINDOW 1,"Primes",(10,30)-(800,640),3

local fn printNumber(theNum&,theColor)
def cycle (1,10,column)
if column = 1 then inc(row)
color theColor
call moveto((Column * 80) + 20,Row * 13)
print TheNum&;" ";
end fn

local fn PrintArray
DIM StartPoint&
TEXT _geneva,12,_BoldBit%,0
Column = 0
Row = 0
fn printNumber(2, _zRed)'A prime (special)
fn printNumber(3, _zRed)'A prime (special)
FOR StartPoint& = 4 TO _MaxNumber step 3
LONG IF fn bittst(gPrimeFlags,StartPoint&/3)
long if StartPoint& and 1
fn printNumber(StartPoint&, _zRed)'A prime
fn printNumber(StartPoint&+1, _zBlack)'Not a prime(even)
fn printNumber(StartPoint&, _zBlack)'Not a prime(even)
fn printNumber(StartPoint&+1, _zRed)'A prime
end if
fn printNumber(StartPoint&, _zBlack)'Not a prime
fn printNumber(StartPoint&+1, _zBlack)'Not a prime
end if
fn printNumber(StartPoint&+2, _zBlack)'Not a prime(Mult of 3)
next StartPoint&
end fn

This routine uses the Sieve of Eratosthenes, which reveals higher 
primes by eliminating all multiples of lower primes. I have added a 
twist that increases speed and reduces memory required.

For any 3 consecutive integers greater than 3, we know that
one will be a multiple of 3. One of the other 2 will be even (i.e., a  
multiple of 2), so it will suffice to record only whether the remaining 
odd integer is or is not a prime number. One bit is adequate to record 
whether a group of 3 includes a prime or not, so a byte will hold enough  
data for 24 consecutive integers, and 40.7kb is all that's required for  
each million integers.

We start by marking every odd integer which is not a multiple of 3 
as prime, then go about clearing that mark for those which are 
multiples of other primes.

When a prime is found, we have already cleared multiples of all the 
lower primes, so the first multiple we need to clear is the square of
the new prime.

From there, we need to clear higher multiples of the new prime which 
are not even and are not multiples of 3. We can do this by alternately
incrementing the square by the prime * 2 and the prime * 4. It works 
out that if the group number is odd, it follows one pattern, and if 
it's even it follows another pattern.

One interesting feature of the code below is that it handles the sieve
without ever calculating the actual values represented by the group 
numbers, thus avoiding time-consuming divisions. The actual value can
be calculated with (num*3+1) OR 1, which you can see in FN PrintPrimes
at the end of the program.

On my 240mHz G3, this will find and record the primes from 1 to 
20,000,000 in 2.5 seconds.

clear local
local fn FindPrimes
dim mult&,num&,carry,increment&
dim StartPoint&,sqr3rd&,StartTick&,EndTick&
StartTick& = FN TickCount'Start timing
sqr3rd& = sqr(_MaxNumber)/3'Number of last group we need to check
def blockfill (@gPrimeFlags,_FlagBytes,&FF)'Initialize all flags as primes

for num = 1 to sqr3rd&'check groups of 3 up to sqrt of _MaxNumber
long if fn bittst(gPrimeFlags,num&)'is there a prime in this group?
long if num and 1'will there be a carry?
mult& = 3*num*num+(num<<2)+1'square it
increment = num+num+1
carry = 0

do'Clear all odd mults of number that aren't mults of 3
call bitclr(gPrimeFlags,mult&)
mult = mult + Increment'add prime * 2
if carry then mult = mult + increment + 1'add prime * 4
carry = 1 - carry
until mult > _max3rd

mult& = 3*num*num+(num<<1)'square it
increment = num+num
carry = 1
call bitclr(gPrimeFlags,mult&)
mult = mult + Increment + 1
if carry then mult = mult + increment
carry = 1 - carry
until mult > _max3rd
end if
end if
IF _MaxNumber <= 450 then fn PrintArray'disable if required

EndTick& = FN TickCount
gRunTime! = (EndTick& - StartTick&) / 60
end fn'FindPrimes
clear local
DIM StartPoint&,Count&
local fn PrintPrimes 'from last 10,000 integers
cursor _watchcursor
print "Counting primes found..."
count& = 2'for 2 and 3
for StartPoint& = 1 to _max3rd
Count& = Count& + fn bittst(gPrimeFlags,StartPoint&)
next StartPoint&

TEXT _monaco,9,0,0
WIDTH = _textWrap
COLOR = _zBlack
print "Found" count& " primes between 2 and" _Max3rd*3+2
delay 1500'time to read

StartPoint& = (_MaxNumber - 10000)/3
long if StartPoint& < 1
StartPoint& = 1
print 2;"  ";3"  ";'special cases
end if
FOR StartPoint& = StartPoint& TO _max3rd
LONG IF fn bittst(gPrimeFlags,StartPoint&)'Is it prime?
print StartPoint&*3 + 1 or 1 ;"  ";'calc & show value
end if
next StartPoint&
cursor _arrowcursor
end fn

FN FindPrimes
if _MaxNumber <= 450 then do : until LEN(INKEY$)'time to look
PRINT "Time Taken for primes up to ";_Max3rd*3+2
PRINT "Is ";gRunTime!;" Seconds"
do : until LEN(INKEY$)
FN PrintPrimes'disable if required
do : until LEN(INKEY$)