under 题解 UVA10214

---

真心推荐刚学莫反的同学来做一下这道题目，这才是真正的 Mobius 反演入门题啊！

P2158 仪仗队是这题的特殊情况，前者可以用欧拉 $\varphi$ 函数解决，而这题用 Mobius 反演简单得多。

没学过莫反的可以来我的这篇题解的 前置知识 部分学习一下。

我们回到这题。

---

分析题意，我们只需将四个象限之一的答案求出，就可推广到整个坐标系。

简化版题意：给定一个 N * M 的矩阵，求

$$\sum \limits_{i = 1}^N \sum \limits_{j = 1}^M [\gcd(i, j) = 1]$$

尝试对这个式子进行套路性的莫反变形，得到

$$\sum \limits_{i = 1}^N \sum \limits_{j = 1}^M \sum \limits_{d | \gcd(i, j)} \mu(d)$$

可以将 $d | \gcd(i, j)$ 拆开

$$\sum \limits_{i = 1}^N \sum \limits_{j = 1}^M \sum \limits_{d | i, d | j} \mu(d)$$

改变枚举顺序，得到

$$\sum \limits_{d = 1}^{\min(N, M)} \mu(d) \lfloor \tfrac{N}{d} \rfloor \lfloor \tfrac{M}{d} \rfloor$$

就可以单次 $O(n)$ 枚举答案了！

但是我们并不知道数据组数 垃圾UVA，所以优化是必要的。和式含有两个形如 $\lfloor \tfrac{N}{d} \rfloor$ 的因子，意味着可以通过数论分块优化到单次 $O( \sqrt{N})$ 的复杂度，可以通过本题。

TIPS:

1. 求出的 ans 是单个象限中的答案，加上坐标轴四个方向的贡献均为 1，因此答案的分子是 ans * 4 + 4。

2. 原点是没有树的，因此分母是 (N * 2 + 1) * (M * 2 + 1) - 1。

3. 极端数据的分母会爆 int，不要忘记 long long 和强制类型转换。

code(c++):

请手动优化常数

#include<algorithm>
#include<bitset>
#include<iomanip>
#include<iostream>

using namespace std;
typedef long long ll;

bitset<2003> v; // 此处作用等价于 bool v[2003];
int N, M, p[305], mu[2003];
ll ans, tot;

inline void sieve(int N) {
  int m = 0;
  mu[1] = 1;
  for(int i = 2; i <= N; ++i) {
    if(!v[i]) p[++m] = i, mu[i] = -1;
    for(int j = 1, x; j <= m && i * p[j] <= N; ++j) {
      x = i * p[j];
      v.set(x); // v[x] = 1;
      if(i % p[j] == 0) break; // 此时 mu[x] = 0;
      mu[x] = -mu[i];
    }
    mu[i] += mu[i - 1]; // 顺便求出前缀和
  }
}

int main() {
  cout << fixed << setprecision(7);
    // 规定输出小数位数
  sieve(2000); // min(maxa, maxb) = 2000
  while(cin >> N >> M) {
    if(!N) return 0;
    ans = 0; // 多测不清空，爆零两行泪
    if(N > M) swap(N, M);
    for(int i = 1, j, n, m; i <= N; i = j + 1) {
      n = N / i, m = M / i; // 可有可无的小常数优化
      j = min(N / n, M / m);
      ans += (ll)(mu[j] - mu[i - 1]) * n * m;
    }
    ans = ans * 4 + 4;
    tot = (ll)(N * 2 + 1) * (M * 2 + 1) - 1;
    cout << (double)ans / tot << '\n';
  }
}