From: Jay Reeve <jktr@...>

Date: Thu, 23 Mar 00 22:53:33 -0600

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. 0"0 =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) xelse fn printNumber(StartPoint&, _zBlack)'Not a prime(even) fn printNumber(StartPoint&+1, _zRed)'A prime end if xelse 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 ' /* HOW IT WORKS 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 xelse mult& = 3*num*num+(num<<1)'square it increment = num+num carry = 1 do 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 next EndTick& = FN TickCount gRunTime! = (EndTick& - StartTick&) / 60 end fn'FindPrimes ' clear local DIM StartPoint&,Count& local fn PrintPrimes 'from last 10,000 integers cls 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& cls 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 BEEP : BEEP if _MaxNumber <= 450 then do : until LEN(INKEY$)'time to look cls 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$)