Skip to content

Sieve Methods

Prime Sieve

Introduction

If we want to know how many primes there are less than or equal to \(n\)?

A natural idea is to perform a primality test for each number less than or equal to \(n\). This brute-force approach obviously cannot achieve optimal complexity.

Sieve of Eratosthenes

Process

Consider this fact: For any positive integer \(n > 1\), its \(x\) times multiple is a composite number (\(x > 1\)). Using this conclusion, we can avoid many unnecessary checks.

If we consider each number from small to large, and simultaneously mark all (larger) multiples of the current number as composite, then the numbers not marked at the end of the run are primes.

Implementation

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
vector<int> prime;
bool is_prime[N];

void Eratosthenes(int n) {
  is_prime[0] = is_prime[1] = false;
  for (int i = 2; i <= n; ++i) is_prime[i] = true;
  for (int i = 2; i <= n; ++i) {
    if (is_prime[i]) {
      prime.push_back(i);
      if ((long long)i * i > n) continue;
      for (int j = i * i; j <= n; j += i)
        // Because multiples from 2 to i - 1 have been sieved before,
        // here we start directly from multiples of i, improving speed.
        is_prime[j] = false;  // Multiples of i are not primes.
    }
  }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
prime = []
is_prime = [False] * N


def Eratosthenes(n):
    is_prime[0] = is_prime[1] = False
    for i in range(2, n + 1):
        is_prime[i] = True
    for i in range(2, n + 1):
        if is_prime[i]:
            prime.append(i)
            if i * i > n:
                continue
            for j in range(i * i, n + 1, i):
                is_prime[j] = False

The above is the Sieve of Eratosthenes, with a time complexity of \(O(n\log\log n)\).

Proof

Now let's look at the derivation process:

If each operation on the array costs 1 unit of time, then the time complexity is:

\[ O\left(\sum_{k=1}^{\pi(n)}{\frac{n}{p_k}}\right)=O\left(n\sum_{k=1}^{\pi(n)}{\frac{1}{p_k}}\right) \]

where \(p_k\) denotes the \(k\)-th smallest prime, and \(\pi(n)\) denotes the number of primes \(\le n\). \(\sum_{k=1}^{\pi(n)}\) represents the first for loop, where the cumulative upper bound \(\pi(n)\) is the number of times if (prime[i]) enters the true branch; \(\frac{n}{p_k}\) represents the number of executions of the second for loop.

According to Mertens' second theorem, there exists a constant \(B_1\) such that:

\[ \sum_{k=1}^{\pi(n)}{\frac{1}{p_k}}=\log\log n+B_1+O\left(\frac{1}{\log n}\right) \]

Therefore, the time complexity of the Sieve of Eratosthenes is \(O(n\log\log n)\). Next, we prove a weaker version of Mertens' second theorem: \(\sum_{k\le\pi(n)}1/p_k=O(\log\log n)\):

Since \(\pi(n)=\Theta(n/\log n)\), the size of the \(n\)-th prime is \(\Theta(n\log n)\). Thus we have

\[ \begin{aligned} \sum_{k=1}^{\pi(n)}{\frac{1}{p_k}} &=O\left(\sum_{k=2}^{\pi(n)}{\frac{1}{k\log k}}\right) \\ &=O\left(\int_2^{\pi(n)}{\frac{\mathrm dx}{x\log x}}\right) \\ &=O(\log\log\pi(n))=O(\log\log n) \end{aligned} \]

Of course, the above approach is still not efficient enough. Applying the following methods can slightly improve the algorithm's performance.

Sieving up to the Square Root

Obviously, to find all primes up to \(n\), it is sufficient to sieve only primes not exceeding \(\sqrt n\).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
vector<int> prime;
bool is_prime[N];

void Eratosthenes(int n) {
  is_prime[0] = is_prime[1] = false;
  for (int i = 2; i <= n; ++i) is_prime[i] = true;
  // i * i <= n means i <= sqrt(n)
  for (int i = 2; i * i <= n; ++i) {
    if (is_prime[i])
      for (int j = i * i; j <= n; j += i) is_prime[j] = false;
  }
  for (int i = 2; i <= n; ++i)
    if (is_prime[i]) prime.push_back(i);
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
prime = []
is_prime = [False] * N


def Eratosthenes(n):
    is_prime[0] = is_prime[1] = False
    for i in range(2, n + 1):
        is_prime[i] = True
    # Let i loop up to <= sqrt(n)
    for i in range(2, isqrt(n) + 1):  # `isqrt` is a new function in Python 3.8
        if is_prime[i]:
            for j in range(i * i, n + 1, i):
                is_prime[j] = False
    for i in range(2, n + 1):
        if is_prime[i]:
            prime.append(i)

This optimization does not affect the asymptotic time complexity. Actually, repeating the above proof, we get \(n \ln \ln \sqrt n + o(n)\). According to logarithmic properties, they are asymptotically the same, but the number of operations is significantly reduced.

Sieving Only Odd Numbers

Since all even numbers except 2 are composite, we can skip them directly and only care about odd numbers.

First, this halves our memory requirement; second, the required operations are approximately halved.

Reducing Memory Usage

We note that only a bool array is needed for sieving. An element of a bool array generally occupies \(1\) byte (i.e., \(8\) bits), but storing a Boolean value requires only \(1\) bit.

We can use knowledge of bit operations to compress each Boolean value into a single bit, so we only need \(n\) bits (i.e., \(\dfrac n 8\) bytes) instead of \(n\) bytes, significantly reducing memory usage. This method is called "bit-level compression."

It is worth mentioning that there are data structures that automatically perform bit-level compression, such as vector<bool> and bitset<> in C++.

Additionally, vector<bool> and bitset<> provide constant optimizations to the program. The Sieve of Eratosthenes with time complexity \(O(n \log \log n)\), when optimized with bitset<> or vector<bool>, can even outperform the Euler sieve with time complexity \(O(n)\).

See bitset: Combined with the Sieve of Eratosthenes.

Segmented Sieving

From the optimization "Sieve up to the Square Root," we know that it is not necessary to keep the entire is_prime[1...n] array all the time. For sieving, it is sufficient to keep primes up to \(\sqrt n\), i.e., prime[1...sqrt(n)]. Then divide the entire range into blocks and sieve each block separately. This way, we don't need to keep multiple blocks in memory simultaneously, and the CPU can handle caching better.

Let \(s\) be a constant determining the block size. Then we have \(\lceil {\frac n s} \rceil\) blocks, and block \(k\) (\(k = 0 \dots \lfloor {\frac n s} \rfloor\)) contains numbers in the interval \([ks, ks + s - 1]\). We can process blocks sequentially, i.e., for each block \(k\), we traverse all primes (from \(1\) to \(\sqrt n\)) and use them for sieving.

It is worth noting that we need to slightly modify the strategy when processing the first block: First, all primes in \([1, \sqrt n]\) should be retained; second, numbers \(0\) and \(1\) should be marked as non-prime. When processing the last block, we should not forget that the last number \(n\) is not necessarily at the end of the block.

The following implementation uses segmented sieving to count the number of primes less than or equal to \(n\).

Implementation 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
int count_primes(int n) {
  constexpr static int S = 10000;
  vector<int> primes;
  int nsqrt = sqrt(n);
  vector<char> is_prime(nsqrt + 1, true);
  for (int i = 2; i <= nsqrt; i++) {
    if (is_prime[i]) {
      primes.push_back(i);
      for (int j = i * i; j <= nsqrt; j += i) is_prime[j] = false;
    }
  }
  int result = 0;
  vector<char> block(S);
  for (int k = 0; k * S <= n; k++) {
    fill(block.begin(), block.end(), true);
    int start = k * S;
    for (int p : primes) {
      int start_idx = (start + p - 1) / p;
      int j = max(start_idx, p) * p - start;
      for (; j < S; j += p) block[j] = false;
    }
    if (k == 0) block[0] = block[1] = false;
    for (int i = 0; i < S && start + i <= n; i++) {
      if (block[i]) result++;
    }
  }
  return result;
}

The asymptotic time complexity of segmented sieving is the same as that of the Sieve of Eratosthenes (unless blocks are very small), but the required memory is reduced to \(O(\sqrt{n} + S)\), and it has better cache performance. On the other hand, division is performed for each pair of block and prime in \([1, \sqrt{n}]\), which is much worse for smaller blocks. Therefore, a balance must be maintained when choosing the constant \(S\).

Choosing block size \(S\) between \(10^4\) and \(10^5\) yields the best speed.

Linear Sieve

The Sieve of Eratosthenes still has room for optimization; it marks a composite number multiple times. Is there a way to eliminate unnecessary steps? The answer is yes.

If each composite number can be marked only once, then the time complexity can be reduced to \(O(n)\).

Implementation 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
vector<int> pri;
bool not_prime[N];

void pre(int n) {
  for (int i = 2; i <= n; ++i) {
    if (!not_prime[i]) {
      pri.push_back(i);
    }
    for (int pri_j : pri) {
      if (i * pri_j > n) break;
      not_prime[i * pri_j] = true;
      if (i % pri_j == 0) {
        // i % pri_j == 0
        // In other words, i has been sieved by pri_j before.
        // Since primes in pri are from small to large, the result of i multiplied by other primes will definitely be
        // sieved by multiples of pri_j, so there is no need to sieve here first, so we can directly break
        // here.
        break;
      }
    }
  }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
pri = []
not_prime = [False] * N


def pre(n):
    for i in range(2, n + 1):
        if not not_prime[i]:
            pri.append(i)
        for pri_j in pri:
            if i * pri_j > n:
                break
            not_prime[i * pri_j] = True
            if i % pri_j == 0:
                """
                i % pri_j == 0
                In other words, i has been sieved by pri_j before.
                Since primes in pri are from small to large, the result of i multiplied by other primes will definitely be
                sieved by multiples of pri_j, so there is no need to sieve here first, so we can directly break
                here.
                """
                break

The above linear sieve is also known as the Euler sieve.

Note

Note that the sieve method for primes also obtains the smallest prime factor of each number.

Sieve for Euler's Totient Function

Note that in the linear sieve, every composite number is sieved by its smallest prime factor. For example, let \(p_1\) be the smallest prime factor of \(n\), \(n' = \frac{n}{p_1}\), then \(n\) is sieved by \(n' \times p_1\) during the linear sieve.

Observing the linear sieve process, we also need to handle two parts. Below we discuss the case of \(n' \bmod p_1\).

If \(n' \bmod p_1 = 0\), then \(n'\) contains all prime factors of \(n\).

\[ \begin{aligned} \varphi(n) & = n \times \prod_{i = 1}^s{\frac{p_i - 1}{p_i}} \\\\ & = p_1 \times n' \times \prod_{i = 1}^s{\frac{p_i - 1}{p_i}} \\\\ & = p_1 \times \varphi(n') \end{aligned} \]

What if \(n' \bmod p_1 \neq 0\)? In this case, \(n'\) and \(p_1\) are coprime. According to the property of Euler's totient function, we have:

\[ \begin{aligned} \varphi(n) & = \varphi(p_1) \times \varphi(n') \\\\ & = (p_1 - 1) \times \varphi(n') \end{aligned} \]

Implementation

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
vector<int> pri;
bool not_prime[N];
int phi[N];

void pre(int n) {
  phi[1] = 1;
  for (int i = 2; i <= n; i++) {
    if (!not_prime[i]) {
      pri.push_back(i);
      phi[i] = i - 1;
    }
    for (int pri_j : pri) {
      if (i * pri_j > n) break;
      not_prime[i * pri_j] = true;
      if (i % pri_j == 0) {
        phi[i * pri_j] = phi[i] * pri_j;
        break;
      }
      phi[i * pri_j] = phi[i] * phi[pri_j];
    }
  }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
pri = []
not_prime = [False] * N
phi = [0] * N


def pre(n):
    phi[1] = 1
    for i in range(2, n + 1):
        if not not_prime[i]:
            pri.append(i)
            phi[i] = i - 1
        for pri_j in pri:
            if i * pri_j > n:
                break
            not_prime[i * pri_j] = True
            if i % pri_j == 0:
                phi[i * pri_j] = phi[i] * pri_j
                break
            phi[i * pri_j] = phi[i] * phi[pri_j]

Sieve for Möbius Function

Definition

According to the definition of the Möbius function, let \(n\) be a composite number, \(p_1\) be the smallest prime factor of \(n\), \(n'=\frac{n}{p_1}\), then:

\[ \mu(n)= \begin{cases} 0 & n' \bmod p_1 = 0\\\\ -\mu(n') & \text{otherwise} \end{cases} \]

If \(n\) is prime, then \(\mu(n)=-1\).

Implementation

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
vector<int> pri;
bool not_prime[N];
int mu[N];

void pre(int n) {
  mu[1] = 1;
  for (int i = 2; i <= n; ++i) {
    if (!not_prime[i]) {
      mu[i] = -1;
      pri.push_back(i);
    }
    for (int pri_j : pri) {
      if (i * pri_j > n) break;
      not_prime[i * pri_j] = true;
      if (i % pri_j == 0) {
        mu[i * pri_j] = 0;
        break;
      }
      mu[i * pri_j] = -mu[i];
    }
  }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
pri = []
not_prime = [False] * N
mu = [0] * N


def pre(n):
    mu[1] = 1
    for i in range(2, n + 1):
        if not not_prime[i]:
            pri.append(i)
            mu[i] = -1
        for pri_j in pri:
            if i * pri_j > n:
                break
            not_prime[i * pri_j] = True
            if i % pri_j == 0:
                mu[i * pri_j] = 0
                break
            mu[i * pri_j] = -mu[i]

Sieve for Number of Divisors

Let \(d_i\) denote the number of divisors of \(i\), and \(num_i\) denote the count of the smallest prime factor of \(i\).

Theorem on Number of Divisors

Theorem: If \(n=\prod_{i=1}^m p_i^{c_i}\), then \(d_i=\prod_{i=1}^m (c_i+1)\).

Proof: We know that \(p_i^{c_i}\) has divisors \(p_i^0,p_i^1,\dots ,p_i^{c_i}\), totaling \(c_i+1\). By the multiplication principle, the number of divisors of \(n\) is \(\prod_{i=1}^m (c_i+1)\).

Implementation

Since \(d_i\) is a multiplicative function, the linear sieve can be used.

Here is a brief introduction to the principle of linear sieve implementation.

  1. When \(i\) is prime, \(\textit{num}_i \gets 1,\textit{d}_i \gets 2\), and let \(q = \left\lfloor \dfrac {i}{p} \right\rfloor\), where \(p\) is the smallest prime factor of \(i\).
  2. When \(p\) is a prime factor of \(q\), \(\textit{num}_i \gets \textit{num}_q + 1,\textit{d}_i \gets \dfrac{\textit{d}_q}{\textit{num}_i} \times (\textit{num}_i + 1)\).
  3. When \(p,q\) are coprime, \(\textit{num}_i \gets 1,\textit{d}_i \gets \textit{d}_q \times (\textit{num}_i+1)\).
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
vector<int> pri;
bool not_prime[N];
int d[N], num[N];

void pre(int n) {
  d[1] = 1;
  for (int i = 2; i <= n; ++i) {
    if (!not_prime[i]) {
      pri.push_back(i);
      d[i] = 2;
      num[i] = 1;
    }
    for (int pri_j : pri) {
      if (i * pri_j > n) break;
      not_prime[i * pri_j] = true;
      if (i % pri_j == 0) {
        num[i * pri_j] = num[i] + 1;
        d[i * pri_j] = d[i] / num[i * pri_j] * (num[i * pri_j] + 1);
        break;
      }
      num[i * pri_j] = 1;
      d[i * pri_j] = d[i] * 2;
    }
  }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
pri = []
not_prime = [False] * N
d = [0] * N
num = [0] * N


def pre(n):
    d[1] = 1
    for i in range(2, n + 1):
        if not not_prime[i]:
            pri.append(i)
            d[i] = 2
            num[i] = 1
        for pri_j in pri:
            if i * pri_j > n:
                break
            not_prime[i * pri_j] = True
            if i % pri_j == 0:
                num[i * pri_j] = num[i] + 1
                d[i * pri_j] = d[i] // num[i * pri_j] * (num[i * pri_j] + 1)
                break
            num[i * pri_j] = 1
            d[i * pri_j] = d[i] * 2

Sieve for Sum of Divisors

\(f_i\) denotes the sum of divisors of \(i\), \(g_i\) denotes \(p^0+p^1+p^2+\dots p^k\) of the smallest prime factor of \(i\).

Implementation

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
vector<int> pri;
bool not_prime[N];
int g[N], f[N];

void pre(int n) {
  g[1] = f[1] = 1;
  for (int i = 2; i <= n; ++i) {
    if (!not_prime[i]) {
      pri.push_back(i);
      g[i] = i + 1;
      f[i] = i + 1;
    }
    for (int pri_j : pri) {
      if (i * pri_j > n) break;
      not_prime[i * pri_j] = true;
      if (i % pri_j == 0) {
        g[i * pri_j] = g[i] * pri_j + 1;
        f[i * pri_j] = f[i] / g[i] * g[i * pri_j];
        break;
      }
      f[i * pri_j] = f[i] * f[pri_j];
      g[i * pri_j] = 1 + pri_j;
    }
  }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
pri = []
not_prime = [False] * N
f = [0] * N
g = [0] * N


def pre(n):
    g[1] = f[1] = 1
    for i in range(2, n + 1):
        if not not_prime[i]:
            pri.append(i)
            g[i] = i + 1
            f[i] = i + 1
        for pri_j in pri:
            if i * pri_j > n:
                break
            not_prime[i * pri_j] = True
            if i % pri_j == 0:
                g[i * pri_j] = g[i] * pri_j + 1
                f[i * pri_j] = f[i] // g[i] * g[i * pri_j]
                break
            f[i * pri_j] = f[i] * f[pri_j]
            g[i * pri_j] = 1 + pri_j

General Multiplicative Functions

Suppose a multiplicative function \(f\) satisfies: For any prime \(p\) and positive integer \(k\), \(f(p^k)\) can be computed in polynomial time with respect to \(k\) of low degree. Then \(f(1),f(2),\dots,f(n)\) can be sieved in \(O(n)\) time.

Let the prime factorization of composite number \(n\) be \(\prod_{i=1}^k p_i^{\alpha_i}\), where \(p_1<p_2<\dots<p_k\) are primes. In the linear sieve, we record \(g_n=p_1^{\alpha_1}\). If \(n\) is sieved by \(x\cdot p\) (\(p\) is prime), then \(g\) satisfies the following recurrence:

\[ g_n= \begin{cases} g_x\cdot p & x\bmod p=0\\\\ p & \text{otherwise} \end{cases} \]

If \(n=g_n\), it means \(n\) is a prime power, and \(f(n)\) can be computed in \(O(1)\); otherwise, \(f(n)=f(\frac{n}{g_n})\cdot f(g_n)\).

Part of this section is translated from the blog post Решето Эратосфена and its English translation Sieve of Eratosthenes. The Russian version is licensed under Public Domain + Leave a Link; the English version is licensed under CC-BY-SA 4.0.

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用