Creating vector of primes using Sieve of Eratosthenes

Question

Please review my Sieve of Eratosthenes implementation.

std::vector<int> generate(const int& n) {
    std::vector<bool> input(n+1, true);
    for (int i = 2; i * i<= n; ++i) {
        for (int j = i; i*j <= n    ; ++j) {
            input[i*j] = false;
        }
    }
    std::vector<int> primes{};
    for (auto i = 1; i <= n; ++i) {
        if (input[i]) {
            primes.push_back(i);
        }
    }
    return primes;
}

How can I improve and make this faster?

Answer 1

From Wikipedia, the Sieve of Eratosthenes algorithm is:

To find all the prime numbers less than or equal to a given integer n by Eratosthenes' method:

1. Create a list of consecutive integers from 2 through n: (2, 3, 4, ..., n).
2. Initially, let p equal 2, the first prime number.
3. Starting from p, enumerate its multiples by counting to n in increments of p, and mark them in the list (these will be 2p, 3p, 4p, etc.; the p itself should not be marked).
4. Find the first number greater than p in the list that is not marked. If there was no such number, stop. Otherwise, let p now equal this new number (which is the next prime), and repeat from step 3.

Your algorithm does not seem right really, notably it skips the step 4 which is essential for fast operation; and there are some non-C++ idioms too

• Do not return the std::vector; instead pass it as a reference to the function.

• Do not pass integers as const references.

• If n < 2, return empty vector.

• The Eratosthenes sieve goes like: "for all values from 2 to sqrt(n) - if the number is a prime then mark all its multiples not prime. However all multiples less than i * i have been marked already composite by previous iterations.

• Do not multiply the loop variables when a simple addition suffices.

• As a simple optimization, the 2 is always in the primes as we skipped the "no values" already; thus we add it. Then we iterate only odd integers starting from 3.

Thus we get

std::vector<int> generate(int n) {
    std::vector<int> result = std::vector<int>();
    if (n < 2) {
        return result;
    }

    // initialize the vector
    std::vector<bool> input(n + 1, true);

    // calculate the upper limit as the square root of
    // of N. all composite numbers <= N must have a 
    // factor <= sqrt(N)
    int sqrtN = (int)sqrt(n);

    // iterate from 2 up to the square root of N
    for (int i = 2; i <= sqrtN; i ++) {
        if (! input[i]) {
            // i is a proven composite number,
            // all its multiples have been 
            // marked not prime by all its prime factors by now.
            continue;
        }

        // as an optimization, all multiples *less* than 
        // the square of i are already marked, (they have 
        // another prime factor less than i), so we can start
        // from the square which is the smallest composite
        // not yet marked
        for (int j = i * i; j <= n; j += i) {
            input[j] = false;
        }
    }

    // as n >= 2, then add 2 here
    result.push_back(2);

    // and check only odd numbers here,
    // no other even number can be set
    for (int i = 3; i <= n; i += 2) {
        if (input[i]) {
            result.push_back(i);
        }
    }

    return result;
}

Of course you do not need the second loop to just fill in the result vector; if the input[i] is true (the continue was not executed), the number is a prime and you can push_back it to the result vector there:

if (! input[i]) {
    // i is an already proven composite number,
    // all its multiples have been 
    // marked not prime by all its prime factors by now.
    continue;
}

// i is a prime
result.push_back(i);

Answer 2

Taking @Antti Haapala as a starting point. Two small optimizations.

- for (int i = 2; i <= sqrtN; i ++) {
+ for (int i = 5; i <= sqrtN; i ++) {

  result.push_back(2);
+ result.push_back(3);

- for (int i = 3; i <= n; i += 2) {
+ for (int i = 5, jump = 2; i <= n; i += jump, jump = 4 - jump) {

The optimization is that you can ignore all values that are multiples of 2 and 3 if you increment by 2 then 4 then 2 then 4 etc...

Now that you know that you are already ignoring multiples of 2 and 3 in your second loop; then in your first loop you don't need to even check for anything below 5. Just make sure to add 3 to the final result.

Answer 3

Nitpickish stuff:

1. You have no input validation, negative n results in an empty result set.

2. When generating the primes I think using the post-increment would be more obvious

   for (auto i = 2; i <= n; i++) {
   

   I take that you want to benefit from the minimal performance increase, but that's more premature optimization, and thus IMO unncessary..

3. Your second for-loop isn't spaced properly: It should be:

   for (int j = i; i * j <= n; ++j) {
   

Less nitpickish stuff:

• The code's logic for eliminating primes is non-obvious. The first thing I thought when I saw your prime-sieve was. "Wouldn't this also cross out 2?". Your code could use a comment here and there ;)

• 1 is not a prime! Since you initialize your input to all true and start reading out at 1, you'll have 1 in the resulting primes.

• A vector can't have a negative amount of elements. As soon as n < -2 you might get problems when initializing. I wouldn't know since I don't do C++, but eh.

Answer 4

You don't need two loops. If it is set up correctly, every time you test a new number, if it is still true, you can output it immediately. And if it is false, you can skip crossing out its multiples. They have already been set to false earlier.

Answer 5

You mark multiples of pimes in the seive with:

for (int j = i * i; j <= n; j += i) {
    input[j] = false;
}

Say i = 5 then you mark 25, 30, 35, 40, 45, 50, 55, 60 ...

You might notice that every second number is a multiple of 2. No need to mark them since they are never tested for inclusion later on. So use j += 2 * i.

If you use the trick to skip multiples of 2 and 3 then you can use that here too, advancing by 2 * i and 4 * i alternatively. But you have to compute with which one to start. Actually I would advance by 6 * i and set input[j] and input[j - 2 * i]. This time you have to compute the right starting point (i * i or (i + 2) * i) instead of the right step.

---

I also have another idea in mined that i always wanted to try and benchmark. In the current code you start by creating a seive of the full size and fully initialize it. Then you mark off multiples of primes. The above answeres already eliminated a lot of duplicate marking but I think one can do better.

Note: There is no need to store even numbers in the seive but it's easier to explain things with them in it. Compressing the seive to skip even numbers is left as exercise.

Start with a seive of 2 * 3 * 5 = 30 and mark it as before (mark 5 * 5 = 25). The next prime is 7. Grow the seive to 2 * 3 * 5 * 7 = 210. The first 30 entries are copied repeadately to fill the seive. Then mark 7. The next prime is 11. Grow the seive to 2 * 3 * 5 * 7 * 11 = 2310, again by copying the first 210 entries. For every prime you find you grow the seive up to the final size you need, which should happen quickly (max 27 for 32bit, 47 for 64bit systems).

Using memcpy to grow the seive could be faster than initializing it at start and the loops to mark individual multiples of primes. It only helps for the first few primes but those are the ones with the most work.