Prime number module

[Stephen Horne]

On 30 Sep 2003 13:00:39 GMT, bokr at oz.net (Bengt Richter) wrote:

>If the object is to be able to test primeness of an arbitrary candidate number
>in range(10**8), I think I might just make a bit map of the 10**8 on disk,
>which would only be
> >>> 10**8/8000000.
> 12.5
>megabytes.

If the object was simply to test primeness I would agree, but the
requirement included determining the number of primes below some given
value. It should be easy enough to fix, though - keep a cached count
for every n primes and only redo the count for the last few.


I think bitsets are the way to go for compression. The idea of only
storing every nth prime until you consider how you implement the seive
for the missing primes. You can't just check the interval - you need
to know the primes in earlier intervals to do the seive checks. You
may end up recalculating a lot of the smaller primes again.


It is not too difficult to count the bits in a single integer. One of
my favorite tricks exploits the fact that a single operation on an
integer can be seen as a parallel operation processing 32 or 64 bits
in one go.

Take an integer. Each bit can be considered a count of the number of
bits in that bit position. With a little bit of masking and shifting,
you can get a value where every pair of bits contains the number of
bits that were in that pair of bits in the original. In the 32 bit
case...

  x = (x & 0x55555555) + ((x & 0xAAAAAAAA) >> 1)

Next, get the totals for each group of four bits...

  x = (x & 0x33333333) + ((x & 0xCCCCCCCC) >> 2)

And so on, until the value contains a single total of the number of
bits in the whole original value.

  x = (x & 0x0F0F0F0F) + ((x & 0xF0F0F0F0) >> 2)
  x = (x & 0x00FF00FF) + ((x & 0xFF00FF00) >> 4)
  x = (x & 0x0000FFFF) + ((x & 0xFFFF0000) >> 8)

A simple looping method would need 32 shifts, 32 bitwise ands and an
average of 16 increments. This takes 5 shifts, 10 bitwise ands and 5
additions. Probably a worthwhile gain.



LuLu points out that you can compress the bitset even further by
storing bits only for odd numbers. This compression can be extended by
looking at larger repeating patterns. For example, if we look at a
each value modulo 2*3 = 6, we can exclude multiples of both two and
three from the bitset.

  0 1 2 3 4 5  - value % 6
  * . * . * .  - multiples of two
  * . . * . .  - multiples of three
  * . * * * .  - multiples of two or three

Basically, you only need to store bits where (value % 6) is one of the
values that gets a dot in the bottom line. Then to determine which bit
you need, you set up something like...

  lookup = [None, 0, None, None, None, 1]

  if lookup [value % 6] is None :
    if (value == 2) or (value == 3) :
      value is prime
    else :
      value is not prime
  else
    lookup_bit ((value / 3) + lookup [value % 6])

So you only need to hold one bit in three (those with value % 6 equal
to either 1 or 5) with both 2 and 3 being special-cased.

Use modulo 2*3*5 == 30 and I think the bitset reduces to 8/30 == just
under 27%, which tells me that along with a rapidly increasing
pattern-checking lookup table size you get rapidly decreasing returns.
No way am I going to work out the figures for modulo 2*3*5*7 = 210!

The modulo 6 case is certainly worthwhile, and maybe the module 30 if
you really need those extra few % points (which are still worth a
couple of MB at a guess), but I doubt it's worth going further than
that.


-- 
Steve Horne

steve at ninereeds dot fsnet dot co dot uk