Skip to content

Min_25 Sieve

Definition

From the perspective of its algorithmic idea, this sieve method is also known as the "Extended Eratosthenes Sieve".

Since it was invented and first used by Min_25, it is called the "Min_25 sieve".

Properties

It can solve the prefix sum problem for a class of multiplicative functions in \(O\left(\frac{n^{\frac{3}{4}}}{\log{n}}\right)\) or \(\Theta\left(n^{1 - \epsilon}\right)\) time complexity.

Requirements: \(f(p)\) is a sum of completely multiplicative functions that can be quickly evaluated with respect to \(p\) (e.g., polynomials); \(f(p^{c})\) can be quickly evaluated.

Notation

  • Unless otherwise specified, all variables denoted by \(p\) in this section take values from the set of all prime numbers.
  • \(x / y := \left\lfloor\frac{x}{y}\right\rfloor\)
  • \(\operatorname{isprime}(n) := [ |\{d : d \mid n\}| = 2 ]\), i.e., its value is \(1\) when \(n\) is prime, and \(0\) otherwise.
  • \(p_{k}\): the \(k\)-th smallest prime among all primes (e.g., \(p_{1} = 2, p_{2} = 3\)). Specifically, let \(p_{0} = 1\).
  • \(\operatorname{lpf}(n) := [1 < n] \min\{p : p \mid n\} + [1 = n]\), i.e., the smallest prime factor of \(n\). Specifically, when \(n=1\), its value is \(1\).
  • \(F_{\mathrm{prime}}(n) := \sum_{2 \le p \le n} f(p)\)
  • \(F_{k}(n) := \sum_{i = 2}^{n} [p_{k} \le \operatorname{lpf}(i)] f(i)\)

Explanation

Observing the definition of \(F_{k}(n)\), we can see that the answer is \(F_{1}(n) + f(1) = F_{1}(n) + 1\).

Consider how to compute \(F_{k}(n)\). By enumerating the smallest prime factor and its exponent of each \(i\), we obtain the recurrence:

\[ \begin{aligned} F_{k}(n) &= \sum_{i = 2}^{n} [p_{k} \le \operatorname{lpf}(i)] f(i) \\ &= \sum_{\substack{k \le i \\ p_{i}^{2} \le n}} \sum_{\substack{c \ge 1 \\ p_{i}^{c} \le n}} f\left(p_{i}^{c}\right) ([c > 1] + F_{i + 1}\left(n / p_{i}^{c}\right)) + \sum_{\substack{k \le i \\ p_{i} \le n}} f(p_{i}) \\ &= \sum_{\substack{k \le i \\ p_{i}^{2} \le n}} \sum_{\substack{c \ge 1 \\ p_{i}^{c} \le n}} f\left(p_{i}^{c}\right) ([c > 1] + F_{i + 1}\left(n / p_{i}^{c}\right)) + F_{\mathrm{prime}}(n) - F_{\mathrm{prime}}(p_{k - 1}) \\ &= \sum_{\substack{k \le i \\ p_{i}^{2} \le n}} \sum_{\substack{c \ge 1 \\ p_{i}^{c + 1} \le n}} \left(f\left(p_{i}^{c}\right) F_{i + 1}\left(n / p_{i}^{c}\right) + f\left(p_{i}^{c + 1}\right)\right) + F_{\mathrm{prime}}(n) - F_{\mathrm{prime}}(p_{k - 1}) \end{aligned} \]

The last step of derivation is based on the fact that for \(c\) satisfying \(p_{i}^{c} \le n < p_{i}^{c + 1}\), we have \(p_{i}^{c + 1} > n \iff n / p_{i}^{c} < p_{i} < p_{i + 1}\), so \(F_{i + 1}\left(n / p_{i}^{c}\right) = 0\).
The boundary value is \(F_{k}(n) = 0 (p_{k} > n)\).

Assuming all \(F_{\mathrm{prime}}(n)\) have been computed, there are two ways to compute all \(F_{k}(n)\):

  1. Directly compute according to the recurrence.
  2. Enumerate \(p\) from large to small for transitions. The increment of the transition is nonzero only when \(p^{2} < n\), so we can optimize using suffix sums according to the recurrence.

Now consider how to compute \(F_{\mathrm{prime}}{(n)}\).
Observing the process of computing \(F_{k}(n)\), it is easy to see that \(F_{\mathrm{prime}}\) only has useful point values at \(1, 2, \dots, \left\lfloor\sqrt{n}\right\rfloor, n / \sqrt{n}, \dots, n / 2, n\), totaling \(O(\sqrt{n})\) points.
Generally, \(f(p)\) is a low-degree polynomial in \(p\), which can be expressed as \(f(p) = \sum a_{i} p^{c_{i}}\).
Then for each \(p^{c_{i}}\), its contribution to \(F_{\mathrm{prime}}(n)\) is \(a_{i} \sum_{2 \le p \le n} p^{c_{i}}\).
Considering the contribution of each \(p^{c_{i}}\) separately, the problem transforms into: given \(n, s, g(p) = p^{s}\), for all \(m = n / i\), compute \(\sum_{p \le m} g(p)\).

Note

\(g(p) = p^{s}\) is a completely multiplicative function!

Thus, define \(G_{k}(n) := \sum_{i = 2}^{n} \left[p_{k} < \operatorname{lpf}(i) \lor \operatorname{isprime}(i)\right] g(i)\), i.e., the sum of \(g\) values of the numbers remaining after the \(k\)-th round of the sieve of Eratosthenes.
For a composite number \(x \le n\), we must have \(\operatorname{lpf}(x) \le \sqrt{x} \le \sqrt{n}\). Let \(p_{\ell(n)}\) be the largest prime not exceeding \(\sqrt{n}\), then \(\sum_{2\le p\le n}g(p) = G_{\ell(n)}(n)\), i.e., after \(\ell\) rounds of the sieve, the remaining numbers are all primes. Consider the boundary value of \(G\), which is clearly \(G_{0}(n) = \sum_{i = 2}^{n} g(i)\). (Remember? We specially defined \(p_{0} = 1\).)
For the transition, consider the process of the sieve of Eratosthenes, discussing the contribution of each part separately:

  1. For the part where \(n < p_{k}^{2}\), the \(G\) value remains unchanged, i.e., \(G_{k}(n) = G_{k - 1}(n)\).
  2. For the part where \(p_{k}^{2} \le n\), the numbers being sieved must have the prime factor \(p_{k}\), i.e., \(-g(p_{k}) G_{k - 1}(n / p_{k})\).
  3. For the second part, since \(p_{k}^{2} \le n \iff p_{k} \le n / p_{k}\), the \(i\) satisfying \(\operatorname{lpf}(i) < p_{k}\) will be subtracted extra. This part should be added back, i.e., \(g(p_{k}) G_{k - 1}(p_{k - 1})\).

Thus:

\[ G_{k}(n) = G_{k - 1}(n) - \left[p_{k}^{2} \le n\right] g(p_{k}) (G_{k - 1}(n / p_{k}) - G_{k - 1}(p_{k - 1})) \]

Complexity Analysis

For the computation of \(F_{k}(n)\), the time complexity of its first method is proven to be \(O\left(n^{1 - \epsilon}\right)\) (see zzt's national team paper 2.3);
For the second method, its essence is the second part of the Zhouge sieve, also mentioned in the Zhouge paper (6.5.4), and its time complexity is proven to be \(O\left(\frac{n^{\frac{3}{4}}}{\log{n}}\right)\).

For the computation of \(F_{\mathrm{prime}}(n)\), in fact, its implementation is the same as the first part of the Zhouge sieve.
Considering that for each \(m = n / i\), only when enumerating \(p_{k}\) satisfying \(p_{k}^{2} \le m\) for transitions will contribute to the time complexity, the time complexity can be estimated as:

\[ \begin{aligned} T(n) &= \sum_{i^{2} \le n} O\left(\pi\left(\sqrt{i}\right)\right) + \sum_{i^{2} \le n} O\left(\pi\left(\sqrt{\frac{n}{i}}\right)\right) \\ &= \sum_{i^{2} \le n} O\left(\frac{\sqrt{i}}{\ln{\sqrt{i}}}\right) + \sum_{i^{2} \le n} O\left(\frac{\sqrt{\frac{n}{i}}}{\ln{\sqrt{\frac{n}{i}}}}\right) \\ &= O\left(\int_{1}^{\sqrt{n}} \frac{\sqrt{\frac{n}{x}}}{\log{\sqrt{\frac{n}{x}}}} \mathrm{d} x\right) \\ &= O\left(\frac{n^{\frac{3}{4}}}{\log{n}}\right) \end{aligned} \]

For space complexity, it can be observed that both \(F_{k}\) and \(F_{\mathrm{prime}}\) only take effective point values at \(n / i\), totaling \(O(\sqrt{n})\) points. Only recording these effective values can optimize the space complexity to \(O(\sqrt{n})\).

First, all effective values can be obtained through a single number-theoretic blocking, recorded in an array \(\text{lis}\) of size \(O(\sqrt{n})\). For an effective value \(v\), let \(\text{id}(v)\) be the index of \(v\) in \(\text{lis}\). It is easy to see that for all effective values \(v\), \(\text{id}(v) \le \sqrt{n}\).

Then consider effective values less than or equal to \(\sqrt{n}\) and those greater than \(\sqrt{n}\) separately: for effective values \(v \le \sqrt{n}\), use an array \(\text{le}\) to record \(\text{id}(v)\), i.e., \(\text{le}_v = \text{id}(v)\); for effective values \(v > \sqrt{n}\), use an array \(\text{ge}\) to record \(\text{id}(v)\). Since \(v\) is too large, we use \(v' = n / v < \sqrt{n}\) to record \(\text{id}(v)\), i.e., \(\text{ge}_{v'} = \text{id}(v)\).

Thus, two arrays of size \(O(\sqrt{n})\) can be used to record the \(\text{id}\) of all effective values and support \(O(1)\) queries. When computing \(F_{k}\) or \(F_{\mathrm{prime}}\), using the \(\text{id}\) of the effective value as the subscript instead of the effective value itself can optimize the space complexity to \(O(\sqrt{n})\).

Procedure

For the computation of \(F_{k}(n)\), we usually implement the first method due to its lower implementation difficulty, and it often performs better than the second method for small data sizes;

For the computation of \(F_{\mathrm{prime}}(n)\), it can be directly implemented according to the recurrence.

For \(p_{k}^{2} \le n\), linear sieve can be used to precompute \(s_{k} := F_{\mathrm{prime}}(p_{k})\) to replace \(F_{\mathrm{prime}}(p_{k - 1})\) in the recurrence of \(F_{k}\).
Correspondingly, \(G_{k - 1}(p_{k - 1}) = \sum_{i = 1}^{k - 1} g(p_{i})\) in the recurrence of \(G\) can also be precomputed using this method.

When using the Extended Eratosthenes Sieve to compute the prefix sum of a multiplicative function \(f\), the following points should be clarified:

  • How to quickly (usually in linear time) sieve the first \(\sqrt{n}\) values of \(f\);
  • The polynomial representation of \(f(p)\);
  • How to quickly compute \(f(p^{c})\).

After clarifying the above points, implement the following parts in order:

  1. Sieve the primes within \([1, \sqrt{n}]\) and the first \(\sqrt{n}\) values of \(f\);
  2. For each term in the polynomial representation of \(f(p)\), sieve the corresponding \(G\), and combine them to obtain all \(O(\sqrt{n})\) useful point values of \(F_{\mathrm{prime}}\);
  3. Implement the recursion according to the recurrence of \(F_{k}\) to compute \(F_{1}(n)\).

Examples

Luogu P4213【Template】Du's Sieve

Compute \(\displaystyle \sum_{i = 1}^{n} \varphi(i)\) and \(\displaystyle \sum_{i = 1}^{n} \mu(i)\).

Solution

For computing the prefix sum of \(\varphi(i)\), first note that \(f(p) = p - 1\). For the linear term \((p)\) of \(f(p)\), we have \(g(p) = p, G_{0}(n) = \sum_{i = 2}^{n} g(i) = \frac{(n + 2) (n - 1)}{2}\); for the constant term \((-1)\) of \(f(p)\), we have \(g(p) = -1, G_{0}(n) = \sum_{i = 2}^{n} g(i) = -n + 1\). Sieve twice and add them up to obtain all \(O(\sqrt{n})\) required point values for \(F_{\mathrm{prime}}\).

For computing the prefix sum of \(\mu(i)\), note that \(f(p) = -1\). Then \(g(p) = -1, G_{0}(n) = \sum_{i = 2}^{n} g(i) = -n + 1\). Directly sieve to obtain all \(O(\sqrt{n})\) required point values for \(F_{\mathrm{prime}}\).

LOJ 6053 Simple Function

Given \(f(n)\):

\[ f(n) = \begin{cases} 1 & n = 1 \\ p \operatorname{xor} c & n = p^{c} \\ f(a)f(b) & n = ab \land a \perp b \end{cases} \]

Compute \(\displaystyle \sum_{i = 1}^{n} f(i)\).

Solution

Note that \(f(p) = p - 1 + 2[p = 2]\). Then sieve similarly to the method for \(\varphi\), with a special case for \(2\).

Reference Code
  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
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
/* 「LOJ #6053」Simple Function */
#include <cmath>
#include <iostream>

constexpr int MAXS = 200000;  // 2 * sqrt(n)
constexpr int mod = 1000000007;

template <typename x_t, typename y_t>
void inc(x_t &x, const y_t &y) {
  x += y;
  (mod <= x) && (x -= mod);
}

template <typename x_t, typename y_t>
void dec(x_t &x, const y_t &y) {
  x -= y;
  (x < 0) && (x += mod);
}

template <typename x_t, typename y_t>
int sum(const x_t &x, const y_t &y) {
  return x + y < mod ? x + y : (x + y - mod);
}

template <typename x_t, typename y_t>
int sub(const x_t &x, const y_t &y) {
  return x < y ? x - y + mod : (x - y);
}

template <typename _Tp>
int div2(const _Tp &x) {
  return ((x & 1) ? x + mod : x) >> 1;
}

// The above functions are for preventing negative numbers and modulo operations.
template <typename _Tp>
long long sqrll(const _Tp &x) {  // Square function
  return (long long)x * x;
}

int pri[MAXS / 7], lpf[MAXS + 1], spri[MAXS + 1], pcnt;

void sieve(const int &n) {
  for (int i = 2; i <= n; ++i) {
    if (lpf[i] == 0) {  // Record prime numbers
      lpf[i] = ++pcnt;
      pri[lpf[i]] = i;
      spri[pcnt] = sum(spri[pcnt - 1], i);  // Prefix sum
    }
    for (int j = 1, v; j <= lpf[i] && (v = i * pri[j]) <= n; ++j) lpf[v] = j;
  }
}

long long global_n;
int lim;
int le[MAXS + 1],  // x <= sqrt(n)
    ge[MAXS + 1];  // x > sqrt(n)
#define idx(v) (v <= lim ? le[v] : ge[global_n / v])

int G[MAXS + 1][2], Fprime[MAXS + 1];
long long lis[MAXS + 1];
int cnt;

void init(const long long &n) {
  for (long long i = 1, j, v; i <= n; i = n / j + 1) {
    j = n / i;
    v = j % mod;
    lis[++cnt] = j;
    (j <= lim ? le[j] : ge[global_n / j]) = cnt;
    G[cnt][0] = sub(v, 1ll);
    G[cnt][1] = div2((long long)(v + 2ll) * (v - 1ll) % mod);
  }
}

void calcFprime() {
  for (int k = 1; k <= pcnt; ++k) {
    const int p = pri[k];
    const long long sqrp = sqrll(p);
    for (int i = 1; lis[i] >= sqrp; ++i) {
      const long long v = lis[i] / p;
      const int id = idx(v);
      dec(G[i][0], sub(G[id][0], k - 1));
      dec(G[i][1], (long long)p * sub(G[id][1], spri[k - 1]) % mod);
    }
  }
  /* F_prime = G_1 - G_0 */
  for (int i = 1; i <= cnt; ++i) Fprime[i] = sub(G[i][1], G[i][0]);
}

int f_p(const int &p, const int &c) {
  /* f(p^{c}) = p xor c */
  return p ^ c;
}

int F(const int &k, const long long &n) {
  if (n < pri[k] || n <= 1) return 0;
  const int id = idx(n);
  long long ans = Fprime[id] - (spri[k - 1] - (k - 1));
  if (k == 1) ans += 2;
  for (int i = k; i <= pcnt && sqrll(pri[i]) <= n; ++i) {
    long long pw = pri[i], pw2 = sqrll(pw);
    for (int c = 1; pw2 <= n; ++c, pw = pw2, pw2 *= pri[i])
      ans +=
          ((long long)f_p(pri[i], c) * F(i + 1, n / pw) + f_p(pri[i], c + 1)) %
          mod;
  }
  return ans % mod;
}

using std::cin;
using std::cout;

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  cin >> global_n;
  lim = sqrt(global_n);  // Upper limit

  sieve(lim + 1000);  // Preprocessing
  init(global_n);
  calcFprime();
  cout << (F(1, global_n) + 1ll + mod) % mod << '\n';

  return 0;
}