upd：256 天后的更新，修复了错误和公式格式，充实了内容

under 题解 P4450

蒟蒻第一道莫反题

本题解适合新手学习基础，dalao 请跳过，不要嘲讽 qwq

---

前置知识：

单位函数（$\epsilon$）：

$$\epsilon(i) = \begin{cases} 1 &, i = 1 \\ 0 &, i \ne 1 \end{cases}$$

即当且仅当 i = 1 时取值取值为 1，否则为 0。很简单对不对？

Dirichlet 卷积（$\ast$）：

Dirichlet 卷积是数论函数之间的运算。即

$$h = f \ast g$$

其中 $f, g, h$ 均为数论函数。则

$$h(n) = \sum_{d | n}{f(d) g \left(\frac nd \right)}$$

也可以写作

$$h(n) = \sum_{i · j = n}{f(i) g(j)}$$

我们发现：$\epsilon$ 函数是 Dirichlet 卷积的单位元，即 $f = f \ast \epsilon$。

Mobius 函数（$\mu$）：

$$\mu(n) = \begin{cases} 1 &, n = 1 \\ (-1)^m &, n = p_1p_2 \cdots p_m \\ 0 &, \mathrm{otherwise} \end{cases}$$

其中 $p_1,p_2, \cdots, p_m$ 是不同的质数。$\mu(n)$ 恰在 $n$ 无平方因子时取值非零。显然 $\mu$ 是积性函数。

$$\sum_{d | n}{\mu(d)} = \epsilon(n) \\ \epsilon = 1 \ast \mu$$

Mobius 变换

设 $f$ 是数论函数，若数论函数 $g$ 满足

$$\begin{align} g &= f * \mathbf1 \\ \textrm{That is} \quad g(n) &= \sum_{d | n}{f(d)} \end{align}$$

则称 $g$ 是 $f$ 的 Mobius 变换。

推一下式子。

$$g = f \ast 1 \Rightarrow f = f \ast \epsilon = f \ast 1 \ast \mu = g \ast \mu$$

即

$$g = f \ast 1 \Leftrightarrow f = g \ast \mu$$

至此就可以来做这道题了。

---

回到了原题。题目要求的是

$$\sum_{i = 1}^A{\sum_{j = 1}^B{[\gcd(i, j) = d]}}$$

对式子进行套路性的变形

$$\sum_{i = 1}^{\lfloor A / d \rfloor}{\sum_{j = 1}^{\lfloor B / d \rfloor}{[\gcd(i, j) = 1]}}$$

出现了 $[ \gcd(i, j) = 1]$，也即 $\epsilon( \gcd(i, j))$。把 $\epsilon = 1 \ast \mu$ 代入，得到了

$$\sum_{i = 1}^{\lfloor A / d \rfloor}{\sum_{j = 1}^{\lfloor B / d \rfloor}{\sum_{g | \gcd(i, j)}{\mu(g)}}}$$

（似乎越来越丑了）此时可以把 $d | \gcd(i, j)$ 拆开：

$$\begin{align} & \sum_{i = 1}^{\lfloor A / d \rfloor}{\sum_{j = 1}^{\lfloor B / d \rfloor}{\sum_{g | i, g | j}{\mu(g)}}} \\ =& \sum_{g = 1}^{\min(A / d, B / d)}{\mu(g) \left\lfloor \frac{A}{dg} \right\rfloor \left\lfloor \frac{B}{dg} \right\rfloor} \end{align}$$

至此，计算答案可以 $O(n)$ 解决。

算法还可以继续优化。发现式子里还有两个形如 $\lfloor n / d \rfloor$ 的因子，因此可以整除分块优化成 $O(\sqrt{n})$ 的。配合 $O(n)$ 的线性筛 $\mu$ 函数，总复杂度 $O(n)$，可以轻松通过本题。

---

code 1

std::bitset<1000003> v;
int p[78501], mu[1000003];

void sieve(int N) {
  int m = 0;

  v.set(1), mu[1] = 1;
  for (int i = 2; i <= N; ++i) {
    if (!v[i]) p[++m] = i, mu[i] = -1;
    for (int j = 1; j <= m && i * p[j] <= N; ++j) {
      v.set(i * p[j]);
      if (i % p[j] == 0) break;
      mu[i * p[j]] = -mu[i];
    }
    mu[i] += mu[i - 1];
  }
}

int main() {
  int A, B, d;
  long long ans = 0;

  std::cin >> A >> B >> d;
  A /= d, B /= d;
  if (A > B) std::swap(A, B);

  sieve(A);
  for (int l = 1, r, a, b; l <= A; l = r + 1) {
    a = A / l, b = B / l;
    r = std::min(A / a, B / b);
    ans += 1ll * a * b * (mu[r] - mu[l - 1]);
  }
  std::cout << ans << '\n';
  return 0;
}

code 2

如果你像我一样闲得难受，可以学来杜教筛写写这题，复杂度好像还更高了。

不会杜教筛但想学的朋友可以来看我的这篇题解。

这个多次询问的题目让我对杜教筛及空间优化有了更深入的理解（

std::bitset<10003> v, vm;
int MX, base, p[1231], mu[10003], mu_[10003];

inline int to(int x) {
  return MX / x;
}

void sieve(int N); // 同 code 1

int getMu(int N) {
  if (N <= base) return mu[N];
  if (vm[to(N)]) return mu_[to(N)];
  int ans = 0;
  for (int l = 2, r, L; l <= N; l = r + 1) {
    L = N / l, r = N / L;
    ans += (r - l + 1) * getMu(L);
  }
  vm.set(to(N));
  return mu_[to(N)] = 1 - ans;
}

int main() {
  int A, B, d;
  long long ans = 0;

  std::cin >> A >> B >> d;
  A /= d, B /= d;
  if (A > B) std::swap(A, B);

  MX = A, base = pow(A, 2. / 3);
  sieve(base);
  for (int l = 1, r, a, b, gR, gL; l <= A; l = r + 1) {
    a = A / l, b = B / l;
    r = std::min(A / a, B / b);
    vm.reset(), gR = getMu(r);
    vm.reset(), gL = getMu(l - 1);
    ans += 1ll * a * b * (gR - gL);
  }
  std::cout << ans << '\n';
  return 0;
}