Sieve of Eratosthenes with Wheel Factorization

4
votes

I'm implementing a reasonably fast prime number generator and I obtained some nice results with a few optimizations on the sieve of Eratosthenes. In particular, during the preliminary part of the algorithm, I skip all multiples of 2 and 3 in this way:

template<class Sieve, class SizeT>
void PrimeGenerator<Sieve, SizeT>::factorize()
{
    SizeT c = 2;
    m_sieve[2] = 1;
    m_sieve[3] = 1;

    for (SizeT i=5; i<m_size; i += c, c = 6 - c)
        m_sieve[i] = 1;
}

Here m_sieve is a boolean array according to the sieve of Eratosthenes. I think this is a sort of Wheel factorization only considering primes 2 and 3, incrementing following the pattern 2, 4, 2, 4,..

What I would like to do is to implement a greater wheel, maybe considering primes 2,3 and 5.

I already read a lot of documentation about it, but I didn't see any implementation with the sieve of Eratosthenes... a sample code could help a lot, but also some hints would be nice :) Thanks.

4
c++algorithmprimessieve-of-eratostheneswheel-factorization
What is your question? Are you looking for sieve of erathosthenes code? – Jean-Bernard Pellerin
stackoverflow.com/questions/8341295/… has a Haskell implementation in answers. – zw324
I'm sorry if i wasn't clear: I would like to skip all multiples of 2,3 and 5, like I already do for multiples of 2 and 3. And no, I don't need sieve of eratosthenes code – Crybot
Separate the mapping of prime candidates (numbers not divisible by you wheels factors) to bits in the sieve from the actual sieving algorithm. A wheel with 2, 3 and 5 is especially convenient, because you get 8 candidates for each segment of 30 numbers, which can conveniently be stored in a byte. – starblue

4 Answers

4
votes

You can go even further. Here is some OCaml code I wrote a few years ago:

let eratosthene borne =
  let remove_multiples a lst =
    let rec remmult multa li accu = function
        []         -> rev accu
      | head::tail ->
          if multa = head
          then remmult (a*(hd li)) (tl li)  accu      tail
          else remmult   multa        li (head::accu) tail
    in
    remmult (a * a) lst [] lst
  in
  let rec first_primes accu ll =
    let a = hd ll in 
    if a * a > borne then (rev accu) @ ll 
    else first_primes (a::accu) (remove_multiples a (tl ll))
  in
  let start_list =
(* Hard code of the differences of consecutive numbers that are prime*)
(* with 2 3 5 7 starting with 11... *) 
    let rec lrec = 2 :: 4 :: 2 :: 4 :: 6 :: 2 :: 6 :: 4 :: 2 :: 4 :: 6
      :: 6 :: 2 :: 6 :: 4 :: 2 :: 6 :: 4 :: 6 :: 8 :: 4 :: 2 :: 4 :: 2
      :: 4 :: 8 :: 6 :: 4 :: 6 :: 2 :: 4 :: 6 :: 2 :: 6 :: 6 :: 4 :: 2
      :: 4 :: 6 :: 2 :: 6 :: 4 :: 2 :: 4 :: 2 :: 10 :: 2 :: 10 :: lrec 
    and listPrime2357 a llrec accu =
      if a > borne then rev accu
      else listPrime2357 (a + (num (hd llrec))) (tl llrec) (a::accu)
    in
    listPrime2357 (num 11) lrec []
  in
  first_primes [(num 7);(num 5);(num 3);(num 2)] start_list;;

Note the nice trick that OCaml allows for cyclic linked list.

3
votes

Paul Pritchard, an Australian mathematician working for IBM, developed a series of wheel sieves in the 1980s. I discuss one of them at my blog, including examples worked by hand and an implementation in Scheme. It's too big to discuss here. You should be aware that though the asymptotic complexity is less than the Sieve of Eratosthenes, implementation details typically make it slower in practice.

2
votes

2*3*5 = 30
spokes = 2,3,4,5,6,8,9,10,12,15,16,18,20,24,30
numbers between spokes: 1,7,11,13,17,19,23,25,29 

int[] gaps = [6,4,2,4,2,4,2,4];
int[] primes = [2,3,5];
int max = 9001;
int counter, max_visited;
while(max_visited < max) {
  int jump = gaps[counter];
  counter = counter + 1 % gaps.length;
  max_visited += jump;
}

Does that help?

Alternatively, this might have been what you wanted instead, pseudo-code:

primes = [2,3,5];
product = multiply(primes);//imaginary function that returns 30
wheel = new int[product];
foreach(prime in primes)
  for(int x = 1; x <= product/prime; x++)
    wheel[prime*x] = 1;
return wheel
0
votes

There's a much better optimisation of it (it takes O(n) operations instead of O(n log log n) for your variant): https://stackoverflow.com/a/17550147/2559094 , but it takes more memory (uses n + n/log n memory instead of n).

If you still want to go on with your optimisation and consider primes p1, p2, ..., pn, you should write all numbers from 0 to p1*p2*...*pn (use lcm in case you decide to use not only prime numbers) and check which of them satisfy the following system: x != 0 (mod p1) x != 0 (mod p2) ... x != 0 (mod pn)

Then you should find all the gaps between these numbers and make an array of those gaps to use them.