Min_25 sieve

Definition

From the perspective of the ideas behind this sieve, it is also called the “Extended Eratosthenes Sieve.”

Because it was invented and first used by Min\_25, it is called the “Min\_25 sieve.”

Properties

It can compute prefix sums of a class of multiplicative functions in $O\left(\frac{n^{\frac{3}{4}}}{\log{n}}\right)$ or $\Theta\left(n^{1 - \epsilon}\right)$ time.

The requirements are that $f(p)$ is a sum of completely multiplicative functions that can be evaluated quickly as a function of $p$ (for example, polynomials), and that $f(p^{c})$ can be evaluated quickly.

Notation

• Unless otherwise stated, every variable denoted by $p$ in this section ranges over all primes.
• $x / y := \left\lfloor\frac{x}{y}\right\rfloor$
• $\operatorname{isprime}(n) := [ |\{d : d \mid n\}| = 2 ]$, namely, when $n$ is prime, its value is $1$, and otherwise it is $0$.
• $p_{k}$: the $k$-th smallest prime among all primes (for example, $p_{1} = 2, p_{2} = 3$). In particular, let $p_{0} = 1$.
• $\operatorname{lpf}(n) := [1 < n] \min\{p : p \mid n\} + [1 = n]$, that is, the smallest prime factor of $n$. In particular, 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

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

Let us consider how to compute $F_{k}(n)$. Enumerating the smallest prime factor of each $i$ and its exponent gives the recurrence:

The last derivation is based on the following fact: when $p_{i}^{c} \le n < p_{i}^{c + 1}$ holds for $c$, 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$. Its boundary value is $F_{k}(n) = 0 (p_{k} > n)$.

Suppose that all $F_{\mathrm{prime}}(n)$ values have already been computed. There are two ways to compute all $F_{k}(n)$ values:

1. Compute them directly according to the recurrence.
2. Enumerate $p$ from large to small and perform the transitions. The increment is nonzero only when $p^{2} < n$, so the recurrence can be optimized with suffix sums.

Now consider how to compute $F_{\mathrm{prime}}{(n)}$.
Looking at the process of computing $F_{k}(n)$, it is easy to see that $F_{\mathrm{prime}}$ has useful values only at the points $1, 2, \dots, \left\lfloor\sqrt{n}\right\rfloor, n / \sqrt{n}, \dots, n / 2, n$, of which there are $O(\sqrt{n})$.
In general, $f(p)$ is a low-degree polynomial in $p$, which can be written as $f(p) = \sum a_{i} p^{c_{i}}$.
Thus, the contribution of each $p^{c_{i}}$ 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 becomes: given $n, s, g(p) = p^{s}$, for every $m = n / i$, find $\sum_{p \le m} g(p)$.

Now define $G_{k}(n) := \sum_{i = 2}^{n} \left[p_{k} < \operatorname{lpf}(i) \lor \operatorname{isprime}(i)\right] g(i)$, namely, after the $k$-th round of the Eratosthenes sieve, the sum of the $g$ values of the numbers that remain.
For a composite number $x \le n$, we necessarily 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)$, meaning that after $\ell$ rounds of the Eratosthenes sieve, only primes remain. Considering the boundary value of $G$, it is clearly $G_{0}(n) = \sum_{i = 2}^{n} g(i)$. (Recall that we specially defined $p_{0} = 1$.)
For the transition, consider the Eratosthenes sieve process and separate the contributions of each part:

1. For the part with $n < p_{k}^{2}$, the value of $G$ does not change, so $G_{k}(n) = G_{k - 1}(n)$.
2. For the part with $p_{k}^{2} \le n$, every removed number has $p_{k}$ as a prime factor, contributing $-g(p_{k}) G_{k - 1}(n / p_{k})$.
3. In the second part, because $p_{k}^{2} \le n \iff p_{k} \le n / p_{k}$, the condition $\operatorname{lpf}(i) < p_{k}$ is met by numbers $i$. This part must be added back, giving $g(p_{k}) G_{k - 1}(p_{k - 1})$.

Thus,

Complexity analysis

For computing $F_{k}(n)$, the time complexity of the first method has been proved to be $O\left(n^{1 - \epsilon}\right)$ (see Section 2.3 of Zhen-Ting Zhu's training-team paper 《Some Special Summation Problems for Number-Theoretic Functions》);
for the second method, its essence is the second part of the Du sieve, which is also mentioned in Ren Zhizhou's paper 《Several Methods for Summing Multiplicative Functions》 (Section 6.5.4), and its time complexity has been proved to be $O\left(\frac{n^{\frac{3}{4}}}{\log{n}}\right)$.

For computing $F_{\mathrm{prime}}(n)$, the implementation is in fact the same as the first part of the Du sieve.
For each $m = n / i$, only enumerating a value satisfying $p_{k}^{2} \le m$ and then processing the prime $p_{k}$ contributes to the time complexity, so the time complexity can be estimated as follows:

For the space complexity, observe that both $F_{k}$ and $F_{\mathrm{prime}}$ take effective values only at $n / i$, for a total of $O(\sqrt{n})$ points. Recording only the effective values reduces the space complexity to $O(\sqrt{n})$.

First, a single number-theoretic block decomposition gives all effective values; record them in an array of size $O(\sqrt{n})$ called $\text{lis}$. For an effective value $v$, let $\text{id}(v)$ be the index of $v$ in $\text{lis}$. It is easy to see that for every effective value $v$, $\text{id}(v) \le \sqrt{n}$.

Then consider separately effective values no greater than $\sqrt{n}$ and those greater than $\sqrt{n}$: for effective values no greater than $\sqrt{n}$ with value $v$, use an array $\text{le}$ to record its $\text{id}(v)$, namely $\text{le}_v = \text{id}(v)$; for effective values greater than $\sqrt{n}$ with value $v$, use an array $\text{ge}$ to record $\text{id}(v)$. Since $v$ is too large, use $v' = n / v < \sqrt{n}$ to record $\text{id}(v)$, namely $\text{ge}_{v'} = \text{id}(v)$.

Thus, two arrays of size $O(\sqrt{n})$ can record the $\text{id}$ of every effective value and support $O(1)$ queries. When computing $F_{k}$ or $F_{\mathrm{prime}}$, use the $\text{id}$ of an effective value as the index instead of the value itself, thereby reducing the space complexity to $O(\sqrt{n})$.

Procedure

For computing $F_{k}(n)$, in implementations we generally choose the first method because it is easier to implement and often performs better than the second method when the data size is small.

For computing $F_{\mathrm{prime}}(n)$, implement the recurrence directly.

For $p_{k}^{2} \le n$, a linear sieve can preprocess $s_{k} := F_{\mathrm{prime}}(p_{k})$ to replace $F_{k}$ in the recurrence for $F_{\mathrm{prime}}(p_{k - 1})$.
Correspondingly, in the recurrence for $G$, the value $G_{k - 1}(p_{k - 1}) = \sum_{i = 1}^{k - 1} g(p_{i})$ can also be preprocessed in this way.

When using the Extended Eratosthenes Sieve to compute the prefix sum of a multiplicative function $f$, the following points must be clear:

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

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

1. Sieve the primes in $[1, \sqrt{n}]$ and the first $\sqrt{n}$ values of $f$;
2. For every term in the polynomial representation of $f(p)$, sieve out the corresponding $G$ and combine them to obtain all useful values of $F_{\mathrm{prime}}$ numbered $O(\sqrt{n})$;
3. Implement the recurrence for $F_{k}$ recursively to obtain $F_{1}(n)$.

Examples

/* "LOJ #6053" Simple Function */
#include <cmath>
#include <iostream>

constexpr int MAXS = 200000;  // 2sqrt(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;
}

// all of the above are to avoid negative numbers and handle modular reduction
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 the prime
      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 bound

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

  return 0;
}