题目链接：Codeforces 1139D

Vivek 最初有一个空数组 $a$ 和一个整数 $m$。接下来他会进行如下操作：

1. 随机选择一个 $1$ 到 $m$ 之间的整数 $x$ 并将它加入到数组 $a$ 的最后。
2. 计算出数组 $a$ 中元素的最大公约数。
3. 如果最大公约数等于 $1$ 那么退出操作。
4. 否则回到步骤 $1$。

请你计算出数组 $a$ 的期望长度，对 $10 ^ 9 + 7$ 取模。

数据范围：$1\le m \le 10 ^ 5$。

Solution

首先我们定义 $\text{DP}$ 状态：$f(i)$ 表示在 $\gcd$ 为 $i$ 时，将其变成 $1$ 的期望操作次数。可以发现答案就是

$$ 1 + \sum_{i = 1} ^ m \frac{f(i)}{m} $$

接下来考虑转移方程。我们可以枚举下一个数字，得到方程：

$$ f(n) = 1 + \sum_{i = 1} ^ m \frac{f(\gcd(i, n))}{m} $$

直接转移的复杂度是 $O(m ^ 2)$ 的。发现 $\gcd$ 可能会出现很多重复的，并且 $\gcd(i, n)$ 一定是 $n$ 的因子，我们考虑用 $g(n, d)$ 表示满足 $1\le i \le m, \gcd(i, n) = d$ 的 $i$ 的数量。这样转移方程化为：

$$ f(n) = 1 + \sum_{d\mid n} \frac{f(d)\cdot g(n,d)}{m} $$

注意这里的一个细节：$d = n$ 时等式两侧都有 $f(n)$，因此我们将这个情况特殊考虑。先将等式两侧乘上 $m$ 变成：

$$ m\cdot f(n) = m + \sum_{d\mid n} f(d)\cdot g(n,d) $$

将 $d = n$ 的情况提出来：

$$ m\cdot f(n) = m + f(n)\cdot g(n, n) + \sum_{d \mid n, d\neq n} f(d)\cdot g(n, d) $$

其中 $g(n, n) = \left\lfloor\frac{m}{n}\right\rfloor​$，移项得到：

$$ \left(m - \left\lfloor\frac{m}{n}\right\rfloor\right)\cdot f(n) = m + \sum_{d \mid n, d\neq n} f(d)\cdot g(n, d) $$

现在问题变成求 $g(n, d)​$，我们将其求和式列出来：

$$ \begin{align*} g(n, d) & = \sum_{i = 1} ^ m [\gcd(i, n) = d] \\ & = \sum_{i = 1} ^ {\left\lfloor\frac{m}{d}\right\rfloor} \left[\gcd\left(i, \frac{n}{d}\right) = 1\right] \\ \end{align*} $$

我们将其莫比乌斯反演，并变换求和顺序：

$$ \begin{align*} g(n, d) & = \sum_{i = 1} ^ {\left\lfloor\frac{m}{d}\right\rfloor} \sum_{t\mid i, t\mid \frac{n}{d}} \mu(t) \\ & = \sum_{t\mid\frac{n}{d}} \sum_{i = 1} ^ {\left\lfloor\frac{m}{d}\right\rfloor}[t\mid i] \\ & = \sum_{t\mid\frac{n}{d}} \left\lfloor\frac{m}{dt}\right\rfloor \end{align*} $$

将它重新代入原式：

$$ \left(m - \left\lfloor\frac{m}{n}\right\rfloor\right)\cdot f(n) = m + \sum_{d \mid n, d\neq n} f(d)\sum_{t\mid\frac{n}{d}} \left\lfloor\frac{m}{dt}\right\rfloor $$

其实对于这个式子，我们预处理每个数的因子然后暴力计算复杂度是正确的。

最后证明一下复杂度：

$$ \begin{align*} & \sum_{i = 1} ^ m \sum_{j \mid i} d(j) \\ = & \sum_{i = 1} ^ m d(i)\cdot\left\lfloor\frac{m}{i}\right\rfloor \\ = & \mathcal O(m ^ {\frac{4}{3}}\log m) \end{align*} $$

时间复杂度：$\mathcal O(m ^ {\frac{4}{3}}\log m)$。

Code

#include <cstdio>
#include <vector>

const int N = 1e5 + 5;
const int P = 1e9 + 7;

int n, tot, p[N], mu[N], inv[N], f[N];
bool flg[N];
std::vector<int> v[N];

void add(int &x, int y) {
    (x += y) >= P && (x -= P);
}
void sieve(int n) {
    mu[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!flg[i]) p[++tot] = i, mu[i] = -1;
        for (int j = 1; j <= tot && i * p[j] <= n; j++) {
            flg[i * p[j]] = true;
            if (i % p[j] == 0) {
                mu[i * p[j]] = 0;
                break;
            } else {
                mu[i * p[j]] = -mu[i];
            }
        }
    }
    for (int i = 1; i <= n; i++) {
        add(mu[i], P);
        for (int j = i + i; j <= n; j += i) v[j].push_back(i);
    }
    inv[1] = 1;
    for (int i = 2; i <= n; i++) {
        inv[i] = 1LL * (P - P / i) * inv[P % i] % P;
    }
}
int main() {
    scanf("%d", &n);
    sieve(n);
    for (int i = 2; i <= n; i++) {
        f[i] = n;
        for (int d: v[i]) {
            int s = 1LL * n / i * mu[i / d] % P;
            for(int t: v[i / d]) {
                add(s, 1LL * n / d / t * mu[t] % P);
            }
            add(f[i], 1LL * f[d] * s % P);
        }
        f[i] = 1LL * f[i] * inv[n - n / i] % P;
    }
    int ans = 0;
    for (int i = 1; i <= n; i++) add(ans, f[i]);
    ans = 1LL * ans * inv[n] % P;
    add(ans, 1);
    printf("%d\n", ans);
    return 0;
}