Extending Sieve of Eratosthenes beyond a billion

Question

For MAX value of 1000000000 (\$10^9\$), it takes about 45 seconds to reach the line where it prints "done". How can I speed this up?

Is printing to screen always going to take so much time? I know it's hardware.

On increasing the size of MAX beyond \$10^9\$ gives me an exception on the line memset(a, true, MAX);. Is there any limit to this function? I should be able to use all of the RAM. Running this program uses around 954 MB of memory.

void sieve_of_eratosthenes(){
    bool* a;
    a = (bool*)malloc(MAX * sizeof(bool));
    memset(a, true, MAX);
    unsigned long int i = 1;
    while (i < MAX){

        while (((++i)<MAX)&&(!a[i]));

        if (2 * i >= MAX)
            break;

        for (unsigned long int k = 2 * i; k < MAX; k += i)
            if (a[k])
                a[k] = false;
    }
    std::cout << "done\n";
    for (unsigned long int i = 2; i < MAX; i++)
        if (a[i])
            std::cout << i << "\n";
    getchar(); getchar();
}

Answer 1

There is a huge number of ways to make this run faster.

First you save an awful lot of space by storing only odd numbers in the sieve. The only even prime is the number 2. Note that on a modern computer the time your algorithm takes is roughly equivalent to the amount of data that it reads and writes, so halving the space needed will half the execution time. So the bool at index \$i\$ represent the number \$2_i+1\$.

Now you only need to remove odd numbers. When you set everything to true initially, set array[0] to false since 1 is not prime. Now when you remove multiples of a prime p, you know that all multiples less than \$p^2\$ have already been removed, so you start the loop for removing primes at \$p^2\$. And you increase by 2p at a time, because p*p, (p+2)*p, (p+4)*p etc. are the odd multiples of p. And since the first prime you remove is \$p^2\$, you can stop the loop looking for primes when \$p^2\$ <= MAX.

Now things are a bit tricky: You look for the next prime p. For that you check a [i] until you find one that is true. The index i maps to the prime p = 2i + 1. \$p^2\$ is an odd number and the number \$p^2\$ is stored at the index j = (p^2 - 1) / 2. So you clear the numbers a [j], a [j + p] etc.

So far it's simple, now we get a bit more clever. Your MAX is big, say \$10^9\$. When you remove the multiples of 3, you get a repeating pattern. For the numbers p = 1, 3, 5, 7, 9, 11, 13, 15, 17, you get the pattern (true, false, true), (true, false, true), (true, false, true) and so on. Three values repeating forever. Fill the first 15 numbers with that pattern, then remove the multiples of five: From (T,F,T, T,F,T, T,F,T, T,F,T, T,F,T) after removing the multiples of five you keep (T,F,F, T,F,T, T,F,T, T,F,T, F,F,T). This pattern of fifteen bools will repeat forever. You make 7 copies giving 105 numbers and remove the multiples of 7. You make 11 copies giving 1155 numbers and remove the multiples of 11. You make 13 copies giving 15,015 numbers and remove the multiples of 13. You may make 17 copies giving 255,255 numbers and remove the multiplies of 17. This is very quick because you didn't have a billion numbers, only a few hundred thousand. When you're done you duplicate the data into the whole array with memcpy. For the last part you need to watch out not to overwrite the end of the array. That was only six numbers, but these 6 numbers do a very significant part of the work!

For the other primes, you would remove all their odd multiples. We can do that faster. Take p = 101. You would remove 101p, 103p, 105p, 107p and so on. But 105p is divisible by 3, so it has been removed already. Same for 111p and so on. So here is what you do: You remove \$p^2\$. If p+2 is not divisible by 3, you also remove (p+2)*p. The next number you would try to remove would be a multiple of 3. So instead of increasing the number you remove by 2p each time, you increase by 4p (avoiding the multiple of 3), then by 2p, then again by 4p, 2p, 4p, 2p and so on. That saves one third of the work.

One more improvements: Instead of using one byte to store each number, you use only one bit. This means you need to perform shifting operations to access each bit, but the amount of memory used shrinks by a factor 1. That lets you handle much bigger numbers, and less memory will work quicker.

And the final, big one: Optimise memory accesses. Let's say you have MAX = 1 billion, and use one bit per odd number = 62.5 megabyte. That's more than fits into your processor cache. Let's say you have 2MB L3 cache = 32 million numbers. In that case you perform the sieve operation completely for the first 32 million numbers. This will run a lot faster because your data is all within the L3 cache. Let's say you have 256KB L2 cache = 4 million numbers. In that case you perform the sieve operation completely for the first 4 million numbers, then the next four million, and so on. This is even faster because now all the data you use is within the L2 cache and can be read / written very quickly.

Answer 2

A few things you can do:

• Use a bit vector (with bit manipulation) instead of a bool array. This gives you 8 times more memory. If you use C++ vector< bool> you get this optimization for free.

• Store only odd numbers in the array. This will save additional factor of 2 in memory. You need to adapt the logic of the program a little bit, but this is not extremely difficult. This saves also some time because in the first iteration the algorithms just marks half of the numbers.

• The outer loop (the while (i < MAX)) only needs to run from 0 to ceil(sqrt(MAX)).

• The inner loop can start at i*i not just 2*i.

Summarizing the first two will give you factor of 16 more memory (or factor of 16 more numbers to sieve with given memory) and the last two will significantly improve the run-time (not sure by how much exactly).

The segfault is most likely due to the fact that a == NULL, i.e. the malloc routine fails to allocate the memory. You need to check for NULL after the malloc call and before memset.

Answer 3

Although it's best to use container classes in C++ instead of allocating memory manually, I'll still point this out since it should still be said:

a = (bool*)malloc(MAX * sizeof(bool));

In C++, it's best to use new instead of malloc(). This will also avoid the need for such a cast, which is required in C++ not not in C.

You can also initialize a instead of declaring it first and then assigning to it.

bool* a = new bool[MAX];

But you should still use a container class, as mentioned by @Andreas H. Plus, you never used free() here, so this function will cause a memory leak.

Answer 4

Is there some reason you have to use the Sieve? I haven't tried it with the packing approaches suggested to cut memory use but the simple test everything up to the square root route blows away the implementation of the Sieve you're using. As soon as you start causing cache misses the performance of the Sieve nosedives and even with the bit-packing there's no way you'll keep it in the L1 cache.