# Extending Sieve of Eratosthenes beyond a billion

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();
}

• to support large numbers, use segmented sieve of Eratosthenes, see primesieve. – jfs Dec 20 '15 at 5:11
• How can you speed up the search of primes on a site without mathjax! Haha GL. – Alec Teal Dec 20 '15 at 18:48
• I have rolled back the last edit. Please see what you may and may not do after receiving answers. If you'd like the new version to be reviewed, you could ask a new question and just link back to this original version. – SuperBiasedMan Dec 21 '15 at 9:31
• @AlecTeal actually... code review does have mathjax enabled :) $ does the trick – Vogel612 Dec 21 '15 at 9:36 • Are you compiling for 32-bit mode? ~950MB is in the ballpark of size limit I've seen for memory allocation in a 32-bit process. One of the major reasons for 64-bit ISAs is to address more memory. 32-bit x86 has been obsolete for a long time. Also, use calloc to get already-zeroed memory, and have your sieve set it to true. Then you can avoid a costly memset. – Peter Cordes Dec 30 '19 at 11:07 ## 4 Answers 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.

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.

• It's probably worth noting that C++ has a std::bitset class that is basically equivalent to std::vector<bool>. – wil93 Dec 20 '15 at 23:33

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.

• bool* a = new bool[MAX]; doesn't allow me to have array size more than 10^8. I used vector<bool>, it slowed me down. Where to use free? – piepi Dec 19 '15 at 22:03
• @piepi: What about with a container class, such as std::vector (with int, not bool)? free() should be called at the end of the function, before returning. – Jamal Dec 19 '15 at 22:05
• Don't you think int will take larger space than bool? – piepi Dec 19 '15 at 22:29
• @piepi: Maybe, but there are actually issues with std::vector<bool> specifically. – Jamal Dec 19 '15 at 22:35
• For an efficient use of memory, I was thinking of using variables with 1 bit only. Is that possible? I read on SO that it might not be efficient to do that cause processors can easily process the packed data rather than when packing and unpacking is involved? They were more pro-int or pro-bool than bit vars. – piepi Dec 19 '15 at 22:40

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.

• Fast segmented sieves pretty much always use segments that fit into the L1 cache; with bigger ranges (far beyond 2^32) the difficulty lies in remembering/managing the working offsets (phase) for each prime from segment to segment. For any given segment, the number of primes for which a multiple occurs in that segments is only a tiny fraction of the primes up to the square root of the segment's upper limit, and hence a lot of time is spent on iteration over 'uninteresting' primes. A 32KB segment close to 2^64 still contains only 32K odd numbers but 203 million prime factors potentially apply. – DarthGizka Dec 20 '15 at 9:39
• Hence the trick is to use bucketing so that a given working offset can be 'posted' to the segment where it will be used again, or at least close to it, such that this offset (or rather its prime) will not have to be considered when working on intervening segments. The gory details can be read on primesieve.org. – DarthGizka Dec 20 '15 at 9:43
• @DarthGizka Huh? "primesieve generates primes about 50 times faster (single-threaded) than an ordinary C/C++ sieve of Eratosthenes implementation and about 10,000 times faster than trial-division." Which would say that trial division is 200x slower than a naive sieve approach--but when I tried it I found trial division blew a naive sieve out of the water at least over the range I tried it on. – Loren Pechtel Dec 20 '15 at 22:18
• Sieves are for answering repeated queries or for generating primes in bulk; if you only need to test a single prime - or a tiny handful - then Trial Division must always be faster (especially as more than 85% of composites are uncovered by the first dozen small primes alone). The reason is that sieves must process all prime factors up to the square root before they can start answering queries whereas Trial Division can stop as soon as the compositeness of the candidate has been proved. Anyway, I was merely addressing your concerns regarding caches and telling you where the real problem lies. – DarthGizka Dec 21 '15 at 6:44
• @DarthGizka I was doing them in bulk, same as the sieve. I forget how high I went but the sieve lost horribly. – Loren Pechtel Dec 21 '15 at 19:37