1

I have this snippet of code that generates the primes on "max" in a sufficient time with Sieve of Eratosthenes.

I want to give the function the posibility to use a starting value to calculate a range of primes. So I wonder at what point in the algoritm I have to hand over the starting value..

e.g.

get_primes(unsigned long from, unsigned long to);
get_primes(200, 5000);

-> Saves the prime numbers from 200 to 5000 in the vector.

Unfortunately I don't understand the algorithm completely. [Especially lines 3 to 5, 7 & 10 are unclear]

I tryed to follow the steps by using a debugger but that also did not make me smarter.

It would be great if anyone can explain me this code better and tell me how to set a start value.

Thank you.

vector<unsigned long long> get_primes(unsigned long max) {
    vector<unsigned long long> primes;
    char *sieve;
    sieve = new char[max / 8 + 1];
    memset(sieve, 0xFF, (max / 8 + 1) * sizeof(char));
    for (unsigned long long x = 2; x <= max; x++)
        if (sieve[x / 8] & (0x01 << (x % 8))) {
            primes.push_back(x);
            for (unsigned long long j = 2 * x; j <= max; j += x)
                sieve[j / 8] &= ~(0x01 << (j % 8));
        }
    delete[] sieve;
    return primes;
}
Jarod42
  • 203,559
  • 14
  • 181
  • 302
Codez
  • 23
  • 1
  • 5
  • 2
    Look at [Sieve_of_Eratosthenes](http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes) – Jarod42 Oct 29 '14 at 13:00
  • `char *sieve` is used as big bit-field. – Jarod42 Oct 29 '14 at 13:07
  • Remove `primes.push_back(x);` from the loop; iterate over the sieve afterwards for your interval. The parts you don't understand use "clever" bit-shifting to store eight flag bits in a `char` (which should be unsigned). – molbdnilo Oct 29 '14 at 13:08
  • You might enjoy my explanation of the segmented Sieve of Eratosthenes: [http://stackoverflow.com/questions/10249378/segmented-sieve-of-eratosthenes](http://stackoverflow.com/questions/10249378/segmented-sieve-of-eratosthenes) – user448810 Oct 29 '14 at 14:19
  • what you describe is an [offset sieve of Eratosthenes](http://stackoverflow.com/a/19641049/849891) (pseudocode and link to C code there) - where you sieve a core up to `sqrt(y)`, and *one* additional segment from `x` to `y` using the primes in the core. "Segmented" sieve refers to a continual discovery of primes, segment after segment. – Will Ness Oct 30 '14 at 13:48

3 Answers3

0

You must start at 2 since the sieve first removes all multiples of 2 to find the next prime as 3. It then removes all multiples of 3 to find the next prime as 5 and so on.

Paul Evans
  • 27,315
  • 3
  • 37
  • 54
0

If you want to generate unsigned long long primes using a version of get_primes() then you are in for a very long wait.

For generating primes in the range lo ... hi (inclusive) you need to consider only factors up to sqrt(hi). Hence you need a small sieve (up to 32 bits) for factors and another small sieve of size (hi - lo + 1) for sieving the target range.

Here's an abbreviated version of a sieve that runs up to 2^64 - 1; it uses a full sieve instead of sieving only odd numbers, because it is the reference code I use to verify optimised implementations. The changes for that are straightforward but add even more pitfalls to the whole shebang. As it is it sieves the 10 million primes between 999560010209 and 999836351599 in about 3 seconds, and those between 18446744073265777349u and 18446744073709551557u (i.e. just below 2^64) in about 20 seconds

The factor sieve is global because it gets reused a lot, and sieving the factors can take a while too. I.e. prepping the factors for a range close to 2^64 means sieving all (or most) of the factors up to 2^32 - 1, and thus it can easily take up to 10 seconds.

I wrapped my bitmap code (moral equivalent to std::bitset<>) and the factor sieve into classes; using raw vectors would make the code inflexible and unreadable. I shortened my code, remove a lot of asserts and other noise, and substituted calls to external functions with inlined code (like the call to std::sqrt()), for the sake of exposition. That way you can cull answers like what to do with the offset (here called lo) directly from verified working code.

The point of having separate number_t and index_t is that number_t can be unsigned long long but index_t must be uint32_t for my current infrastructure. The member functions of bitmap_t use the name of the underlying CPU instructions. BTS ... bit test and set, BT ... bit test. Bitmaps are initialised to 0 and a set bit signifies non-prime.

typedef uint32_t index_t;
sieve_t g_factor_sieve;

template<typename number_t, typename OutputIterator>
index_t generate_primes (number_t lo, number_t hi, OutputIterator sink)
{
   // ...

   index_t max_factor = index_t(std::sqrt(double(hi)));

   g_factor_sieve.extend_to_cover(max_factor);

   number_t range = hi - lo;  assert( range <= index_t(index_t(0) - 1) );
   index_t range32 = index_t(range);
   bitmap_t bm(range32);

   if (lo <  2)  bm.bts(1 - index_t(lo));   // 1 is not a prime
   if (lo == 0)  bm.bts(0);                 // 0 is not a prime

   for (index_t n = 2; n <= max_factor && n > 1; n += 1 + (n & 1))
   {
      if (g_factor_sieve.not_prime(n))  continue;

      number_t start = square(number_t(n));
      index_t stride = n << (int(n) > 2 ? 1 : 0);  // double stride for n > 2

      if (start >= lo)
         start -= lo;
      else
         start = (stride - (lo - start) % stride) % stride;

      // double test because of the possibility of wrapping
      for (index_t i = index_t(start); i <= bm.max_bit; )
      {
         bm.bts(i);

         if ((i += stride) < stride)
         {
            break;
         }
      }
   }

   // output

   for (index_t i = 0; ; ++i)
   {
      if (!bm.bt(i))
      {
         *sink = lo + i;

         ++sink;  
         ++n;
      }

      if (i >= bm.max_bit)  break;
   }

   return n;
}
DarthGizka
  • 4,347
  • 1
  • 24
  • 36
  • I gave the range 999560010209 ... 999836351599 as an example because that's the last full 10-million primes batch that can be downloaded from [primos.mat.br](http://primos.mat.br); the last batch (999836351609 to 999999999989) contains only 5.9 million primes. – DarthGizka Oct 29 '14 at 15:56
-1

Try this one; I used it to set starting and ending numbers

for(int x = m;x<n;x++){
    if(x%2!=0 && x%3!=0 && x%5!=0 && x%7!=0 && x%11!=0)
      // then x is prime
}

where m is starting value, and n is the ending value

kimjb64
  • 1
  • 1