Intitutive method to find prime numbers in a range

Question

While trying to find prime numbers in a range (see problem description), I came across the following code:

_{(Code taken from here)}

// For each prime in sqrt(N) we need to use it in the segmented sieve process.
for (i = 0; i < cnt; i++) {
    p = myPrimes[i]; // Store the prime.
    s = M / p;
    s = s * p; // The closest number less than M that is composite number for this prime p.

    for (int j = s; j <= N; j = j + p) {
        if (j < M) continue; // Because composite numbers less than M are of no concern.

        /* j - M = index in the array primesNow, this is as max index allowed in the array
           is not N, it is DIFF_SIZE so we are storing the numbers offset from.
           while printing we will add M and print to get the actual number. */
        primesNow[j - M] = false;
    }
}

// In this loop the first prime numbers for example say 2, 3 are also set to false.
for (int i = 0; i < cnt; i++) { // Hence we need to print them in case they're in range.
    if (myPrimes[i] >= M && myPrimes[i] <= N) // Without this loop you will see that for a
                                              // range (1, 30), 2 & 3 doesn't get printed.
        cout << myPrimes[i] << endl;
}

// primesNow[] = false for all composite numbers, primes found by checking with true.
for (int i = 0; i < N - M + 1; ++i) {
    // i + M != 1 to ensure that for i = 0 and M = 1, 1 is not considered a prime number.
    if (primesNow[i] == true && (i + M) != 1)
        cout << i + M << endl; // Print our prime numbers in the range.
}

However, I didn't find this code intuitive and it was not easy to understand.

Can someone explain the general idea behind the above algorithm?
What alternative algorithms are there to mark non-prime numbers in a range?

The basic principle is simple; once you have found a prime, all multiples of that number are not primes. — sp2danny, Sep 26 '15 at 21:59

score 0 · Answer 1 · edited May 23 '17 at 12:22

That's overly complicated. Let's start with a basic Sieve of Eratosthenes, in pseudocode, that outputs all the primes less than or equal to n:

function primes(n)
    sieve := makeArray(2..n, True)
    for p from 2 to n
        if sieve[p]
            output(p)
            for i from p*p to n step p
                sieve[p] := False

This function calls output on each prime p; output can print the primes, or sum the primes, or count them, or do whatever you want to do with them. The outer for loop considers each candidate prime in turn; The sieving occurs in the inner for loop where multiples of the current prime p are removed from the sieve.

Once you understand how that works, go here for a discussion of the segmented Sieve of Eratosthenes over a range.

user3715912 · Answer 2 · 2015-10-14T03:25:56.227

Have you considered the sieve on a bit level, it can provide a bit larger number of primes, and with the buffer, you could modify it to find for example the primes between 2 and 2^60 using 64 bit ints, by reusing the same buffer, while preserving the offsets of the primes already discovered. The following will use an array of integers.

Declerations

#include <math.h>         // sqrt(), the upper limit need to eliminate
#include <stdio.h>        // for printing, could use <iostream>

Macros to manipulate bit, the following will use 32bit ints

#define BIT_SET(d, n)   (d[n>>5]|=1<<(n-((n>>5)<<5)))      
#define BIT_GET(d, n)   (d[n>>5]&1<<(n-((n>>5)<<5)))
#define BIT_FLIP(d, n)  (d[n>>5]&=~(1<<(n-((n>>5)<<5))))

unsigned int n = 0x80000;   // the upper limit 1/2 mb, with 32 bits each 
                            // will get the 1st primes upto 16 mb    
int *data = new int[n];     // allocate
unsigned int r = n * 0x20;  // the actual number of bits avalible

Could use zeros to save time but, on (1) for prime, is a bit more intuitive

for(int i=0;i<n;i++)
    data[i] = 0xFFFFFFFF;

unsigned int seed = 2;           // the seed starts at 2  
unsigned int uLimit = sqrt(r);   // the upper limit for checking off the sieve

BIT_FLIP(data, 1);      // one is not prime

Time to discover the primes this took under a half second

// untill uLimit is reached
while(seed < uLimit) {
    // don't include itself when eliminating canidates
    for(int i=seed+seed;i<r;i+=seed) 
        BIT_FLIP(data, i);

    // find the next bit still active (set to 1), don't include the current seed
    for(int i=seed+1;i<r;i++) {
        if (BIT_GET(data, i)) {
            seed = i;
            break;
        }
    }
}

Now for the output this will consume the most time

unsigned long bit_index = 0;      // the current bit
int w = 8;                        // the width of a column
unsigned pc = 0;                  // prime, count, to assist in creating columns

for(int i=0;i<n;i++) {
    unsigned long long int b = 1;  // double width, so there is no overflow

    // if a bit is still set, include that as a result

    while(b < 0xFFFFFFFF) {
        if (data[i]&b) {
            printf("%8.u ", bit_index);
            if(((pc++) % w) == 0)
                putchar('\n');   // add a new row                
        }

        bit_index++;
        b<<=1;                       // multiply by 2, to check the next bit    
    }

}

clean up

delete [] data;

Intitutive method to find prime numbers in a range

2 Answers2