杜教筛

前置知识：Mobius 反演

令 \(f\) 为数论函数，杜教筛可以在 \(O(n^{\frac{2}{3}})\) 的时间内求出 \(S(n) = \sum\limits_{i = 1}^n f(i)\)。

求解过程

模板：洛谷 P4213

因为 \(f\) 不好求，首先我们要找到好求的 \(g\) 和 \(f * g\)，最常见的有 \(\varphi * 1 = \epsilon,\mu * 1 = id\)，其中 \(\epsilon\) 是单位函数，\(id = id(n) = n\)。

\[\begin{aligned} \sum\limits_{i = 1}^n f * g &= \sum\limits_{i= 1}^n \sum\limits_{d | i} g(d)\cdot f(\frac{i}{d}) \\ &= \sum\limits_{d = 1}^n g(d) \sum\limits_{i = 1}^{\frac{n}{d}} f(i) \\ &= g(1)S(n) + \sum\limits_{d = 2}^n g(d)S(\lfloor \frac{n}{d} \rfloor) \end{aligned} \]

由此可知 \(S(n) = \sum\limits_{i = 1} ^ n f * g - \sum\limits_{d =2} g(d)S(\lfloor \frac{n}{d} \rfloor)\)。因为 \(\sum g, \sum\limits_{i= 1}^n f * g\) 易求，只需递归求 \(S(\lfloor \frac{n}{d} \rfloor)\) 即可。

而 \(S(\lfloor \frac{n}{d} \rfloor)\) 可以使用整除分块，这样的时间复杂度为 \(O(n^{\frac{3}{4}})\)，通过欧拉筛求出 \(S(1)\sim S(10^6)\) 可将时间复杂度降至 \(O(n^{\frac{2}{3}})\)。有时可以使用 unordered_map 记忆化搜索降低常数。

时间复杂度证明见 oi-wiki

#include <bits/stdc++.h>

using namespace std;
using ll = long long;
using pli = pair<ll, int>;

const int MAXN = 4e6 + 3;

bool vis[MAXN];
unordered_map<int, pli> pos;
int T, n, k, p[MAXN], mu[MAXN], phi[MAXN], s2[MAXN];
ll s1[MAXN];

void build() { // 预处理
  s1[1] = s2[1] = 1;
  for (int i = 2; i < MAXN; i++) {
    if (!vis[i]) {
      p[++k] = i, mu[i] = -1, phi[i] = i - 1;
    }
    s1[i] = s1[i - 1] + phi[i], s2[i] = s2[i - 1] + mu[i];
    for (int j = 1; j <= k && p[j] * i < MAXN; j++) {
      int x = i * p[j]; vis[x] = 1;
      phi[x] = phi[i] * p[j];
      if (i % p[j] == 0) break;
      phi[x] -= phi[i], mu[x] = -mu[i];
    }
  }
}

pli solve(int n) {
  if (n < MAXN) return {s1[n], s2[n]};
  if (pos.count(n)) return pos[n];
  ll sphi = n * (n + 1ll) / 2, smu = 1;
  for (ll l = 2, r; l <= n; l = r + 1ll) { // 整除分块
    int d = n / l; r = n / d;
    auto [x, y] = solve(d);
    sphi -= x * (r - l + 1), smu -= y * (r - l + 1);
  }
  return pos[n] = {sphi, smu};
}

int main() {
  ios::sync_with_stdio(0), cin.tie(0);
  for (cin >> T, build(); T--; ) {
    cin >> n;
    auto [x, y] = solve(n);
    cout << x << ' ' << y << '\n';
  }
  return 0;
}

题目1：洛谷 P3172

给定 \(n, k, l, r(1 \le n, k, l, r \le 10^9, r- l \le 10^5)\)，求 \(\sum\limits_{i_1 = l}^r\sum\limits_{i_2 = l}^r \cdots \sum\limits_{i_n = l}^r [\gcd(i_1, i_2, \dots,i_n) == k]\)。

令 \(l = \lceil \frac{l}{k} \rceil, r = \lfloor \frac{r}{k} \rfloor\)，则有

\[\begin{aligned} ans &= \sum\limits_{i_1 = l}^r\sum\limits_{i_2 = l}^r \cdots \sum\limits_{i_n = l}^r [\gcd(i_1, i_2, \dots,i_n) == 1] \\ &= \sum\limits_{x = 1} \mu(x)(\lfloor \frac{r}{x} \rfloor - \lceil \frac{l}{x} \rceil) ^ n \end{aligned} \]

杜教筛求 \(\sum \mu(x)\)，然后整出分块即可 （注意有一个是上取整，有一个是下取整）。

#include <bits/stdc++.h>

using namespace std;
using ll = int;

const int MAXN = 3e6 + 3, Mod = 1e9 + 7;

bool vis[MAXN];
unordered_map<int, ll> pos;
int n, k, L, R, ans, p[MAXN], mu[MAXN], f[MAXN], s[MAXN];

int qpow(int x, int k) {
  int ret = 1;
  for (; k; k >>= 1, x = 1ll * x * x % Mod) {
    if (k & 1) ret = 1ll * ret * x % Mod;
  }
  return ret;
}

void build() {
  k = 0, s[1] = mu[1] = f[1] = 1;
  for (int i = 2; i < MAXN; i++) {
    if (!vis[i]) {
      p[++k] = i, mu[i] = -1;
      f[i] = qpow(i, n);
    }
    s[i] = s[i - 1] + mu[i];
    for (int j = 1; j <= k && p[j] * i < MAXN; j++) {
      vis[i * p[j]] = 1, f[i * p[j]] = 1ll * f[i] * f[p[j]] % Mod;
      if (i % p[j] == 0) break;
      mu[i * p[j]] = -mu[i];
    }
  }
}

ll solve(int n) {
  if (n < MAXN) return s[n];
  if (pos.count(n)) return pos[n];
  ll ret = 1;
  for (int l = 2, r; l <= n; l = r + 1) {
    r = n / (n / l);
    ret -= (r - l + 1ll) * solve(n / l);
  }
  return pos[n] = ret;
}

int main() {
  ios::sync_with_stdio(0),cin.tie(0);
  cin >> n >> k >> L >> R;
  L = (L + k - 1) / k, R /= k, build();
  for (int i = 1; i < MAXN; i++) {
    (ans += mu[i] * f[R / i - (L + i - 1) / i + 1]) %= Mod;
  }
  for (int l = MAXN, r; l <= R; l = r + 1) {
    int d1 = (L + l - 1) / l, d2 = R / l; r = R / d2;
    if (d1 > 1) r = min((L + d1 - 2) / (d1 - 1) - 1, r);
    ans = (ans + (0ll + solve(r) - solve(l - 1)) * f[d2 - d1 + 1]) % Mod;
  }
  cout << (ans + Mod) % Mod;
  return 0;
}

习题：洛谷 P3768