Prime Sieving in C++ (part 1)

Part 1 | Part 2

If you are a competitive programmer, you must know how to generate prime numbers quickly. As we all know, prime numbers are positive integers n > 1 that is only divisible by 1 and n. In this post, I'll be discussing various (adequately fast) sieve techniques that can be useful, not only in generating primes, but also in some other number theory stuff like generating divisors, the totient function, and more (see part 2).

What is Sieving?

According to Wikipedia, sieving is a method to filter wanted elements from unwanted ones. In our case, we are dealing with prime numbers. Below, we'll see the classical prime sieve and its extensions.

1. Trial Division O(n^2)

bool isPrime[n] = {false, false};
for (int i = 2; i < n; ++i) {
 isPrime[i] = true;
 for (int j = 2; j < i; ++j) {
  if (i % j == 0) {
   isPrime[j] = false;
   break;
  }
 }
}

Everyone started with the brute force sieve by trial division. This is essentially checking if a number is divisible by everything smaller than it using the modulo (%) operator. This method straightforward, but slow.


This might be O(n^2), but this sieve can run ~1.3 seconds for n = 100,000 on normal machines, which is already doable for programming problems with small primes.

2. Optimal Trial Division O(n * √n)

bool isPrime[n] = {false, false};
for (int i = 2; i < n; ++i) {
 isPrime[i] = true;
 for (int j = 2; j*j <= i; ++j) { // same as j <= sqrt(i)
  if (i % j == 0) {
   isPrime[j] = false;
   break;
  }
 }
}

The Optimal Trial Division Sieve is an optimization of the brute force sieve. To check if n is prime, instead of checking all numbers less than n, we can just check all numbers up to √n.

Why is this still correct? Short proof, we can represent n as n = a * b. Then it follows that either one of them will always be less than or equal to √n. Fancy, right? We can now modify our code above by changing one line.

This fancy optimization can run in just 0.2s for n = 1 million and ~5.4 seconds for n = 10 million (10^7). Can we get any faster?

3. Eratosthenes Sieve O(n * log log n)

bool isPrime[n] = {false, false};
for (int p = 2; p < n; ++p)
 isPrime[p] = true;
for (int p = 2; p < n; ++p)
 if (isPrime[p])
  for (int i = 2*p; i < n; i += p)
   isPrime[i] = false;

The most famous sieve used by competitive programmers is Eratosthenes prime sieve. As you may guess, it's named after the Greek dude who made the algorithm. Wikipedia has this rather fancy animation of how it works.

What does the algorithm do? Basically, it's the reversed version of trial division. Instead of checking if each number is prime through the numbers below it, the sieve starts with the primes you have and crosses out all its multiples above it. First, you start with 2, cross out 4, 6, 8, etc., find the next prime 3 and cross out 6, 9, 12, etc., until you get to the end. Complicated as it may sound, the code itself is rather short.

This code runs in less than 0.2s for n = 10 million (10^7) and ~3.2s for n = 100 million (10^8). WAIT WHAT??? How is this one faster than trial division when it basically does the same thing? Moreover, how the heck is it O(n log log n)?

Don't worry, we have a proof. We might need a little knowledge about sequences and series. Let's look at  how many numbers are being checked per iteration in the second loop. It is easy to show that we don't iterate if the number is composite. But what if the number is prime? If the number p is prime, the number of iterations will be ⌊n / p⌋. Sum up all the iterations, we will have:

But the right part has a prime harmonic series! Such series asymptotically approaches log log n. Done n times, our proof is now complete, it is indeed O(n log log n).

Improving Eratosthenes Sieve O(n * log log n)

Though the complexity roughly remains the same, there have been some improvements done to the sieve. For example, we can just check odd numbers in the outer loop or provide a modulus wheel that only checks 6n + 1 and 6n + 5, or similar. In addition, we can merge use the square root idea optimization with Eratosthenes by crossing out composites starting from p^2 instead of 2*p because a composite number below p^2 must have already been crossed out by a prime below p=√p^2.
There are various space and time optimizations for Eratosthenes sieve (e.g. the Segmented Sieve which only uses √n space). Well, Google is your friend.

4. Sieve of Atkin O(n)

code ommited cuz it's not practical

The Sieve of Atkin is (by far) the fastest prime sieve ever known. It exposes the usage of the modulus wheel of primes and properties of remainders. It does the same crossing out method of the Eratosthenes sieve, but instead of crossing out multiples of a prime, it crosses out multiples of squares of primes with rather complicated remainder pruning techniques. More recently, mathematicians have even found an O(n / log log n) version of Atkin's sieve, which blew the minds of modern number theorists.

The sad thing is, the sieve of Atkin has very complicated rules and is very hassle to code for the sake of competitive programming. Moreover, even though it should be theoretically faster, it is actually slower in practice compared to Eratosthenes' sieve because of frequent cache misses due to skipping remainders. Even the prominent primesieve.org uses Eratosthenes' and not Atkin's because they concluded that the former is indeed faster in practice. For this sake, I won't dwell on the actual code and proof because it is not essential for competitive programming anyway. Well, if you're curious, you can browse online for more formal research papers on Atkin's sieve.

What sieve should I use?

The most practical one is undoubtedly Eratosthenes sieve. It is both fast to run and code, and it can generate primes less than 10^8 in less than 4 seconds on a normal computer (what more if it was a Russian computer run by an online judge!). Not only is it fast, it can also be extended to generate other number theory stuff like number of divisors, sum of divisors, euler totient, prime factorization, and more!

If you want to know how to extend Eratosthenes sieve, you may want to see part 2.