Sieve Methods

Prime Sieves

Introduction

How many primes are there less than or equal to $n$?

A natural idea is to test each number less than or equal to $n$ for primality once. This brute-force approach clearly cannot achieve optimal complexity.

Sieve of Eratosthenes

Process

Consider the following fact: for any positive integer greater than $1$, let it be $n$; its $x$-multiple is composite when $x > 1$. Using this fact, we can avoid many unnecessary tests.

If we process the numbers from small to large and, for each current number, mark all of its multiples greater than itself as composite, then the numbers that remain unmarked when the process ends are exactly the primes.

Implementation

C++

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)
        // We have already sieved multiples of 2 through i - 1,
        // so start directly from multiples of i here to improve the running time.
        is_prime[j] = false;  // None of the multiples of i is prime.
    }
  }
}

Python

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, also known as the Eratosthenes sieve or the Eulerian abbreviation in Chinese usage. Its time complexity is $O(n\log\log n)$.

Sieve Only up to the Square Root

Clearly, to find all primes up to $n$, it is enough to sieve only primes no greater than $\sqrt n$.

C++

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

Python

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 run only while i <= sqrt(n).
    for i in range(2, isqrt(n) + 1):  # `isqrt` was added 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 change the asymptotic time complexity. Repeating the proof above gives $n \ln \ln \sqrt n + o(n)$, which is asymptotically equivalent by the properties of logarithms, but the number of operations is significantly smaller.

Sieve Only Odd Numbers

Since every even number other than 2 is composite, we can skip even numbers and consider only odd numbers.

This halves the memory requirement and also roughly halves the number of required operations.

Reducing Memory Usage

During sieving, we only need an array of type bool. One element of a bool array generally occupies $1$ byte (that is, $8$ bits), although storing a Boolean value requires only $1$ bit.

Using the relevant knowledge of bit operations, we can pack each Boolean value into one bit. We then need only $n$ bits (that is, $\dfrac n 8$ bytes), rather than $n$ bytes, which can significantly reduce memory usage. This method is called bit-level compression.

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

In addition, vector<bool> and bitset<> offer constant-factor optimizations. After optimization with bitset<> or vector<bool>, an Eratosthenes sieve with time complexity $O(n \log \log n)$ can even outperform an Euler sieve with time complexity $O(n)$.

See bitset: Combining It with the Sieve of Eratosthenes.

Segmented Sieve

The optimization of sieving only up to the square root shows that we do not need to keep the entire is_prime[1...n] array at all times. To perform the sieve, it is enough to retain the primes up to $\sqrt n$, namely prime[1...sqrt(n)], and divide the whole range into blocks, sieving each block separately. This avoids keeping multiple blocks in memory at the same time and also lets the CPU make better use of its cache.

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

The first block requires a slight adjustment. First, all primes in $[1, \sqrt n]$ should be retained; second, the numbers $0$ and $1$ should be marked as non-prime. When processing the last block, do not forget that the final number $n$ need not be at the end of the block.

The following implementation uses a segmented sieve to count the primes less than or equal to $n$.

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 a segmented sieve is the same as that of the sieve of Eratosthenes (unless the blocks are extremely small), but its memory usage is reduced to $O(\sqrt{n} + S)$ and it has better cache behavior. On the other hand, it performs a division for every pair consisting of a block and a prime in $[1, \sqrt{n}]$, which is particularly unfavorable for small blocks. Therefore, the constant $S$ must be chosen as a balance between these factors.

A block size $S$ between $10^4$ and $10^5$ generally gives the best speed.

Linear Sieve

The sieve of Eratosthenes can still be optimized: it marks a composite number multiple times. Is there a way to eliminate these meaningless repeated operations? Yes.

If every composite number can be marked exactly once, the time complexity can be reduced to $O(n)$.

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 already been sieved by pri_j.
        // Since the primes in pri are in increasing order, multiplying i by any other prime would produce a number that is certainly
        // sieved by a multiple of pri_j, so it need not be sieved here; simply break
        // here.
        break;
      }
    }
  }
}

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 already been sieved by pri_j.
                Because the primes in pri are ordered increasingly, products of i with later primes will certainly be
                removed as multiples of pri_j. There is no need to sieve them here first, so break
                immediately.
                """
                break

The above linear sieve is also called the Euler sieve.

Computing Euler's Totient Function with a Sieve

In the linear sieve, every composite number is eliminated by its smallest prime factor. For example, let $p_1$ be the smallest prime factor of $n$ and let $n' = \frac{n}{p_1}$. During the linear sieve, $n$ is eliminated through $n' \times p_1$.

Examining the linear sieve process, we need to handle two cases. We distinguish them according to $n' \bmod p_1$.

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

If $n' \bmod p_1 \neq 0$, then $n'$ and $p_1$ are coprime. By the properties of Euler's totient function,

Implementation

C++

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];
    }
  }
}

Python

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]

Computing the Möbius Function with a Sieve

Definition

By the definition of the Möbius function, let $n$ be composite, let $p_1$ be the smallest prime factor of $n$, and let $n'=\frac{n}{p_1}$. Then

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

Implementation

C++

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];
    }
  }
}

Python

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]

Computing the Number of Divisors with a Sieve

Let $d_i$ denote the number of divisors of $i$, and let $num_i$ denote the number of occurrences of the smallest prime factor of $i$.

Divisor-Count Theorem

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

Proof: the divisors of $p_i^{c_i}$ are $p_i^0,p_i^1,\dots ,p_i^{c_i}$, for a total of $c_i+1$. By the multiplication principle, the number of divisors of $n$ is therefore $\prod_{i=1}^m (c_i+1)$.

Implementation

Since $d_i$ is a multiplicative function, we can use a linear sieve.

Here is a brief description of the principle behind the implementation.

1. When $i$ is prime, set $\textit{num}_i \gets 1,\textit{d}_i \gets 2$. Also 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$, set $\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, set $\textit{num}_i \gets 1,\textit{d}_i \gets \textit{d}_q \times (\textit{num}_i+1)$.

C++

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

Python

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

Computing the Sum of Divisors with a Sieve

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

Implementation

C++

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

Python

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$ has the property that, for every prime $p$ and positive integer $k$, in time polynomial in $k$ we can compute $f(p^k)$; then in $O(n)$ time we can sieve the values $f(1),f(2),\dots,f(n)$.

Suppose the prime factorization of a composite number $n$ is $\prod_{i=1}^k p_i^{\alpha_i}$, where $p_1<p_2<\dots<p_k$ are primes. In the linear sieve, record $g_n=p_1^{\alpha_1}$. If $n$ is eliminated by $x\cdot p$, where $p$ is prime, then $g$ satisfies

If $n=g_n$, then $n$ is a power of a prime, so $f(n)$ can be computed in $O(1)$ time. Otherwise, $f(n)=f(\frac{n}{g_n})\cdot f(g_n)$.

Part of this section is adapted from the blog post Решето Эратосфена and its English-language edition 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.