Sieves for primes

Why this matters when you build with AI

Almost every number theory task starts with a list of primes. An AI tool will write a sieve for you in a few seconds, and the code often looks right. The mistakes are at the edges: n = 1, a loop bound that overflows int, a memory use of 125 MB for a sieve you only needed for 10^9 bits, or an index formula that is off by one. These bugs give a wrong count, not a crash, so they pass a quick read.

You can catch them with two checks. First, compare against brute force on small n, where brute force means testing each number by trial division. Second, compare the total with a value you know. Two known values are pi(10^7) = 664579 and pi(10^9) = 50847534, where pi(n) is the number of primes up to n. These come from OEIS sequence A006880, and the programs below reproduce both.

This lesson builds four sieves in C++17, shows where each one breaks, and then uses a segmented sieve on a window of numbers near 10^12. A sieve is a method that marks the numbers it rules out and keeps the ones left unmarked.

The idea: the Sieve of Eratosthenes

Take the numbers from 2 to n and cross out every multiple of 2 after 2. The next number left is 3, which is prime, so cross out its multiples after 3. Continue with each number that is still not crossed out. Whatever stays uncrossed is prime.

Worked example for n = 30. Start with 2 to 30. The first unmarked number is 2, so cross out 4, 6, 8, ..., 30. The next unmarked is 3, so cross out 9, 12, 15, ..., 27. The next is 5, which crosses out 25. The next is 7, and 7 * 7 = 49 is above 30, so the work is done. The unmarked numbers are 2, 3, 5, 7, 11, 13, 17, 19, 23 and 29. That is ten primes, and the count matches the table of pi(n) values.

Two details matter. You start crossing out at i * i, not at 2 * i, because every smaller multiple k * i with k < i has a prime factor below i and was already crossed out. You only need to go while i * i <= n, because a composite number up to n has a prime factor no larger than its square root.

The work done is the sum of n / p over the primes p up to the square root of n. Mertens' theorem gives that sum as about n log log n, which is close to linear for every size you will meet in a contest.

Here is the Python version, which is fast enough for n up to a few million. It is input free, so the site can run it:

Python · editable

Real output of that program:

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
count: 25

The classic sieve in C++

The C++ version stores one byte per number. It uses char instead of bool because vector<bool> packs bits and is slower to index in a tight loop. The ll type is used for i * i so the product cannot overflow when i is large.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
 
// Classic Sieve of Eratosthenes: is_prime[i] for 0 <= i <= n.
vector<char> classic_sieve(int n) {
    vector<char> is_prime(n + 1, 1);
    is_prime[0] = 0;
    if (n >= 1) is_prime[1] = 0;     // 0 and 1 are not prime
    for (ll i = 2; i * i <= n; i++) {
        if (!is_prime[i]) continue;
        for (ll j = i * i; j <= n; j += i) is_prime[j] = 0;
    }
    return is_prime;
}

For n = 10^8 this uses 100 MB, which is more than many judges allow. The next two variants cut the memory down.

Odd-only sieve

Every even number above 2 is composite, so the odd-only sieve never stores them. Index i stands for the odd number 2i + 1. Index 0 is the number 1, which is not prime. The number n maps to index (n - 1) / 2, the largest odd number <= n.

Two index facts make the loop work. The first odd multiple of an odd prime p to cross out is p * p, and its index is (p * p - 1) / 2. The next odd multiple is p * p + 2p, which is p steps further along in index form. So the inner loop is j += p.

// Odd-only sieve: index i stands for the odd number 2i + 1. Returns the number of primes <= n.
ll count_primes_odd_only(ll n) {
    if (n < 2) return 0;
    ll half = (n - 1) / 2;                 // odd numbers 3, 5, ..., up to n are indices 1..half
    vector<char> composite(half + 1, 0);   // composite[i] = 1 means 2i+1 is not prime
    composite[0] = 1;                      // the number 1
    for (ll i = 1; (2 * i + 1) * (2 * i + 1) <= n; i++) {
        if (composite[i]) continue;
        ll p = 2 * i + 1;
        for (ll j = (p * p - 1) / 2; j <= half; j += p) composite[j] = 1;
    }
    ll count = 1;                          // the prime 2
    for (ll i = 1; i <= half; i++)
        if (!composite[i]) count++;
    return count;
}

The bound (n - 1) / 2 is the one line where AI tools often slip. Writing n / 2 looks natural, and it counts the odd number n + 1 when n is even. That number is composite for the values that were tested first, so the check passed for 10^7 and 10^8 by luck. Here is what the wrong bound gives on small inputs, from a program that uses n / 2 and nothing else changed:

odd_only_bad(4) = 3 (correct: 2)
odd_only_bad(6) = 4 (correct: 3)
odd_only_bad(10) = 5 (correct: 4)

The fix was found with an exhaustive check for every n from 0 to 3000 against a brute force, which reported 0 mismatches after the fix. The lesson is to test every small n, not only a few round values.

Segmented sieve

The classic sieve needs memory for the whole range. A segmented sieve splits the range into blocks, usually a few tens of thousands of numbers each, and sieves one block at a time. The primes needed to cross out composites in a block are the base primes, which are all primes up to the square root of the top of the range. There are only a few thousand of them for ranges up to 10^9.

For each block [lo, hi] and each base prime p, find the first multiple of p that is at least lo. Start there, but never below p * p, so that p itself is not crossed out. Then step by p and mark. Counting the unmarked numbers gives the primes in that block.

// Count primes in [L, R]. base must hold every prime p with p * p <= R, in increasing order.
ll count_primes_segmented(ll L, ll R, const vector<ll> &base) {
    if (L < 2) L = 2;
    if (R < L) return 0;
    const ll BLOCK = 1 << 16;
    ll count = 0;
    vector<char> mark;
    for (ll lo = L; lo <= R; lo += BLOCK) {
        ll hi = min(R, lo + BLOCK - 1);
        mark.assign(hi - lo + 1, 1);
        for (ll p : base) {
            if (p * p > hi) break;
            // First multiple of p inside the block, but never p itself, so start at p*p at the earliest.
            ll start = max(p * p, (lo + p - 1) / p * p);
            for (ll m = start; m <= hi; m += p) mark[m - lo] = 0;
        }
        for (ll i = 0; i <= hi - lo; i++) count += mark[i];
    }
    return count;
}

The base primes come from a small sieve up to sqrt(R). The helper below uses the classic sieve, and isqrt_ll returns the exact integer square root for 64-bit values, which sqrtl alone does not guarantee.

// Primes up to r, by the classic sieve. Used to build base primes.
vector<ll> primes_upto(ll r) {
    vector<ll> out;
    if (r < 2) return out;
    vector<char> is_prime = classic_sieve((int)r);
    for (ll i = 2; i <= r; i++)
        if (is_prime[i]) out.push_back(i);
    return out;
}
 
// Integer square root, exact for 64-bit values.
ll isqrt_ll(ll n) {
    ll r = (ll)sqrtl((long double)n);
    while (r * r > n) r--;
    while ((r + 1) * (r + 1) <= n) r++;
    return r;
}

Measured on this machine, one segmented pass over [2, 10^9] returned 50847534 in 1939 ms. That matches the known value of pi(10^9). The timing will differ on your machine.

Linear sieve

The classic sieve crosses out a composite once for each of its prime factors. The number 30 is crossed out by 2, by 3 and by 5. The linear sieve crosses out each composite exactly once, and it records its smallest prime factor while doing it. Gries and Misra published this method in 1978 in the Communications of the ACM. Its running time is linear in n.

The rule is this. For each i from 2 to n, if i has no recorded smallest prime factor, i is prime. Then for each prime p in the list, in increasing order, stop when p is larger than the smallest prime factor of i, or when p * i is above n. Otherwise, record p as the smallest prime factor of p * i.

Why this crosses out each composite once: every composite c has exactly one pair (p, i) with c = p * i, p prime, and p no larger than the smallest prime factor of i. That pair is p = spf(c) and i = c / p. The stop condition is what rules out every other pair.

Trace for n = 30. The first 12 lines of the real output of the trace mode of the program below are:

i=2: prime, record lp[2]=2
i=2: mark 4 = 2 * 2, lp=2
i=3: prime, record lp[3]=3
i=3: mark 6 = 2 * 3, lp=2
i=3: mark 9 = 3 * 3, lp=3
i=4: mark 8 = 2 * 4, lp=2
i=5: prime, record lp[5]=5
i=5: mark 10 = 2 * 5, lp=2
i=5: mark 15 = 3 * 5, lp=3
i=5: mark 25 = 5 * 5, lp=5
i=6: mark 12 = 2 * 6, lp=2
i=7: prime, record lp[7]=7

The full trace has 40 lines, and the last line is primes up to 30: 2 3 5 7 11 13 17 19 23 29. Look at i = 6: it writes 12 with p = 2, then stops before p = 3, because 3 is larger than lp[6] = 2. So 18 is not written from i = 6. It is written once, from i = 9 with p = 2.

// Linear sieve: lp[i] is the smallest prime factor of i, for 2 <= i <= n.
// primes receives every prime <= n in increasing order.
void linear_sieve(int n, vector<int> &lp, vector<int> &primes) {
    lp.assign(n + 1, 0);
    primes.clear();
    for (int i = 2; i <= n; i++) {
        if (lp[i] == 0) {          // no smaller prime divides i, so i is prime
            lp[i] = i;
            primes.push_back(i);
        }
        for (int p : primes) {
            // Stop when p is larger than the smallest prime factor of i, or when p * i > n.
            if (p > lp[i] || (ll)p * i > n) break;
            lp[p * i] = p;         // p is the smallest prime factor of p * i
        }
    }
}

The check compares this sieve with the other three on n = 10^7, and the four counts agree. Real output of the check, compiled with g++ -std=c++17 -O2:

pi(1e7): classic=664579 odd_only=664579 linear=664579 segmented=664579
pi(1e9) segmented = 50847534

Counting the writes made by the loop shows the difference from the broken version that drops the stop condition. For n = 10^7 the linear sieve makes 9,335,420 writes. That is exactly the number of composites from 2 to 10^7, since there are 9,999,999 numbers in that range and 664,579 of them are prime. Without the stop, the same program makes 22,850,051 writes.

Timings on this machine, from two separate runs with g++ -O2, for n = 10^8:

run 1: classic 1876 ms, odd_only 872 ms, linear 1506 ms
run 2: classic 2546 ms, odd_only 1308 ms, linear 2530 ms

The two runs disagree by up to 40 percent, so the numbers show the order of the methods on this machine and nothing more. The odd-only version was fastest in both runs, probably because it touches half as many bytes.

Worked problem: primes in a window near 10^12

The exercise asks for the count of primes in [L, R] for R up to 10^12. Sieving all of [2, 10^12] is not practical. At the speed of the 1939 ms run above, a linear extrapolation gives roughly half an hour for 10^12 numbers on this machine, and that is far beyond any contest time limit. A window of 10^7 numbers is practical, and that is what the solution accepts.

Input: a line with T (at most 10), then T lines each with L and R. Limits: 1 <= L <= R <= 10^12 and R - L <= 10^7. Output: one count per line.

The base primes go up to 10^6, because sqrt(10^12) = 10^6. There are 78,498 of them. They are built once and reused for every query.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
 
vector<ll> primes_upto(ll r) {
    vector<ll> out;
    if (r < 2) return out;
    vector<char> is_prime(r + 1, 1);
    is_prime[0] = 0;
    is_prime[1] = 0;
    for (ll i = 2; i * i <= r; i++)
        if (is_prime[i])
            for (ll j = i * i; j <= r; j += i) is_prime[j] = 0;
    for (ll i = 2; i <= r; i++)
        if (is_prime[i]) out.push_back(i);
    return out;
}
 
// Count primes in [L, R] with one segmented pass. base holds all primes p with p * p <= R.
ll count_primes_segmented(ll L, ll R, const vector<ll> &base) {
    if (L < 2) L = 2;
    if (R < L) return 0;
    const ll BLOCK = 1 << 16;
    ll count = 0;
    vector<char> mark;
    for (ll lo = L; lo <= R; lo += BLOCK) {
        ll hi = min(R, lo + BLOCK - 1);
        mark.assign(hi - lo + 1, 1);
        for (ll p : base) {
            if (p * p > hi) break;
            ll start = max(p * p, (lo + p - 1) / p * p);
            for (ll m = start; m <= hi; m += p) mark[m - lo] = 0;
        }
        for (ll i = 0; i <= hi - lo; i++) count += mark[i];
    }
    return count;
}
 
int main() {
    int T;
    cin >> T;
    // Base primes up to 1e6 cover every R <= 1e12, since sqrt(1e12) = 1e6.
    vector<ll> base = primes_upto(1000000);
    while (T--) {
        ll L, R;
        cin >> L >> R;
        // Allowed input: 1 <= L <= R <= 1e12 and R - L <= 1e7.
        cout << count_primes_segmented(L, R, base) << "\n";
    }
    return 0;
}

Real output on a sample of three queries:

input:
3
1000000000000 1000000000100
1 100
999999999989 999999999989
 
output:
4
25
1

The first window holds four primes. The second is the primes up to 100, which is 25. The third checks a single number, 999999999989, which is prime and gives 1.

How this was tested

The solution was compared with two other methods. On 200 random windows with L up to 5000 and widths up to 3000, the counts matched trial division every time, 0 mismatches. On six windows of 20,001 numbers each, placed near 10^12 - 10^6, the counts matched a deterministic Miller-Rabin test with the first twelve prime bases. Miller-Rabin is a primality test covered in a later lesson, and these bases are known to be correct for all 64-bit inputs. A single window of 10^7 numbers took 0.06 s, including building the base primes, on this machine. Its count was not checked by an independent method, because checking 10^7 numbers that way would take too long in Python.

Complexity

Method	Time	Memory
Classic sieve	`O(n log log n)`	one byte per number, `n` bytes
Odd-only sieve	`O(n log log n)`	half the classic memory
Segmented sieve on `[L, R]`	about `(R - L) log log R + sqrt(R)`	one block, plus base primes
Linear sieve	`O(n)`	`lp` array of `4n` bytes for `int`

The linear sieve is the only one that records a factor for every number, which is what the later lessons on factorisation and divisor counts need. Its memory is four times the classic sieve when the values are int, so for n = 10^8 expect about 400 MB, which may be more than a judge allows.

Common mistakes

Treating 1 as prime. If the sieve starts with is_prime[1] = 1, the count for n = 1 is 1 instead of 0, and every count is off by one. The fix is the line is_prime[1] = 0.

Overflow in i * i with int. The product 46341 * 46341 is 2147488281, which is above the largest int value of 2147483647. Signed overflow is undefined behaviour in C++, so the loop test can misbehave once i reaches 46341. Use long long for the product, as the code above does.

Odd-only bound n / 2. The output of this bug is shown above: odd_only_bad(4) = 3 where the answer is 2. The fix is (n - 1) / 2.

Linear sieve without the stop. The result stays correct, but the loop writes 22,850,051 times for n = 10^7 instead of 9,335,420. The program runs slower and loses the linear time bound. The stop condition is the whole point of the method.

Segmented sieve starting at the multiple p. If the start is ceil(lo / p) * p with no max(p * p, ...), the base prime p itself is crossed out when it lies in the block. The count then loses one prime for each base prime in the window.

vector<bool> for large n. The sieve for n = 10^9 needs 10^9 bits, which is 125 MB, and bit-level access is slower than bytes in the inner loop. For n = 10^9, use the segmented sieve or the odd-only version with bytes for the whole range.

Checklist to keep

Mark 0 and 1 as not prime before the loop.
Start crossing out at i * i, and stop when i * i > n.
Use long long for every product of two indices.
In the odd-only sieve, the largest odd index is (n - 1) / 2.
In the segmented sieve, start each prime's crossing-out at max(p * p, first multiple).
Build base primes up to sqrt(R) with an exact integer square root.
Check the count against a known value, such as pi(10^7) = 664579.
Test every small n against trial division, not only round values.

Hands-on exercise

Build the four sieves in your own scratch folder. Check that the classic, odd-only and linear sieves agree on every n from 0 to 3000 against trial division. Then compile the segmented window program and run three windows of your choice near 10^12, each with width 20,001. Check each count with a Miller-Rabin test that you write yourself, using the bases 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31 and 37.

A good answer shows zero mismatches on every check and reports which method you used for each. Expect roughly 700 to 800 primes in each window of 20,001 numbers near 10^12, because primes thin out as numbers grow. Those counts are not in this lesson, so compare only with your own Miller-Rabin results.

Prompt to try with an AI tool

Use this prompt, then check the answer with the steps above:

Write a C++17 linear sieve for n up to 10^7 that stores the smallest prime factor of every
number from 2 to n. Then add a test that compares it with trial division for every number
from 2 to 3000 and prints the first wrong value, or 0 mismatches if there is none.

A good answer has the stop condition p > lp[i] and the check p * i > n cast to long long. Its test covers every value from 2 to 3000, not a sample. Run the test before trusting the sieve, and if the output reports a wrong value, check the stop condition first.

Practice

Try each one on paper or in your scratch folder before you open the answer.

Practice 1: primes up to 10^6

Count the primes up to 10^6 with any of the sieves.

Hint: The count is a known value from the table of pi(n), and the Python sieve in this lesson can handle this size.

Show the answer

There are 78,498 primes up to 10^6. A Python run of the sieve from this lesson returned 78498. This is also the number of base primes used by the window solution.

Practice 2: primes in a short window

Count the primes in [90, 110].

Hint: List the numbers from 90 to 110 and cross out the even ones and the multiples of 3, 5 and 7.

Show the answer

The primes are 97, 101, 103, 107 and 109, so the count is 5.

Practice 3: smallest prime factors

Use the table of smallest prime factors for 1 to 30. What is lp[25], and what is lp[27]?

Hint: Use the trace above, or the array printed below.

Show the answer

lp[25] = 5 and lp[27] = 3. The full array printed by the linear sieve for 2 to 30 was:

2:2 3:3 4:2 5:5 6:2 7:7 8:2 9:3 10:2 11:11 12:2 13:13 14:2 15:3 16:2 17:17 18:2 19:19 20:2 21:3 22:2 23:23 24:2 25:5 26:2 27:3 28:2 29:29 30:2

Practice 4: odd-only count

Use the odd-only function to count the primes up to 100. Do the index calculation for the number 97 by hand.

Hint: 97 is odd, so its index is (97 - 1) / 2.

Show the answer

The count is 25. The index for 97 is (97 - 1) / 2 = 48. The sieve for n = 3000 was compared with trial division on every value, and the odd-only count matched at each n, so the value 25 at n = 100 is covered by that check.

Practice 5: base primes for a window

Explain how many base primes the window solution needs for R = 10^12, and what the largest one is.

Hint: The base primes are all primes up to the square root of R.

Show the answer

The square root of 10^12 is 10^6. The base primes are all primes up to 10^6, which is 78,498 primes. The largest one is 999983, the largest prime below 10^6.

Summary

Four sieves cover most prime-list needs in a contest. The classic sieve is the simplest, the odd-only sieve halves its memory, the segmented sieve handles windows far above the memory limit, and the linear sieve records a smallest prime factor for every number in linear time. Each one has an edge case that AI tools get wrong, and each is checked best with brute force on small values and a known count on large ones.

Sources and scope

Gries and Misra, "A linear sieve algorithm for finding prime numbers", Communications of the ACM, volume 21, number 12, pages 999 to 1003, 1978, DOI 10.1145/359657.359660. Its method is the one used for the linear sieve above. The known values pi(10^7) = 664579, pi(10^8) = 5761455 and pi(10^9) = 50847534 come from OEIS A006880, and the program reproduced each of them. The statement that there are 78,498 primes up to 10^6 is also from the same OEIS list and matched a Python run. The Mertens estimate is a known result, and the timings come from g++ 15.2.0 on Linux with -std=c++17 -O2 on one machine, so they will differ on yours.