详细题解

1. 问题转化

设 $E[x]$ 表示当前数组的 GCD 为 $x$ 时，还需继续添加的数的期望个数（不包含当前已经存在的数）。特别地，$E[1]=0$，因为 GCD 已经为 $1$ 则立即停止。

初始时数组为空，我们首先随机选择一个数 $x$（$1\le x\le m$），然后进入状态 $x$，因此总期望长度为：

$$\text{ans} = 1 + \frac{1}{m}\sum_{x=1}^{m} E[x].$$

对于 $x>1$，我们考虑下一步：随机选择 $y\in[1,m]$，新的 GCD 变为 $g = \gcd(x,y)$。因此：

$$E[x] = 1 + \frac{1}{m}\sum_{y=1}^{m} E[\gcd(x,y)]. $$

边界 $E[1]=0$。

2. 化简递推式

定义 $c(x,d) = \#\{ y\in[1,m] \mid \gcd(x,y)=d \}$。则上式化为：

$$E[x] = 1 + \frac{1}{m}\sum_{d\mid x} E[d]\cdot c(x,d). $$

注意 $d$ 必须整除 $x$，因为 $\gcd(x,y)$ 一定是 $x$ 的约数。

当 $d=x$ 时，$c(x,x) = \left\lfloor\frac{m}{x}\right\rfloor$（即所有 $x$ 的倍数）。此时方程中包含 $E[x]$ 自身项，移项得：

$$E[x]\left(1 - \frac{c(x,x)}{m}\right) = 1 + \frac{1}{m}\sum_{d\mid x, d<x} E[d]\cdot c(x,d). $$

于是可以解出：

$$E[x] = \frac{1 + \frac{1}{m}\sum_{d\mid x, d<x} E[d]\cdot c(x,d)}{1 - \frac{c(x,x)}{m}}. $$

由于 $x>1$ 时 $c(x,x)=\lfloor m/x\rfloor < m$，分母非零。

3. 计算 $c(x,d)$

利用莫比乌斯函数：

$$c(x,d) = \#\{ y: d\mid y,\ \gcd(x/d,\ y/d)=1 \} = \sum_{k\mid (x/d)} \mu(k) \left\lfloor\frac{m}{d\cdot k}\right\rfloor. $$

记 $t = x/d$，则 $c(x,d) = \sum_{k\mid t} \mu(k) \left\lfloor\frac{m}{d k}\right\rfloor$。

4. 算法步骤

• 预处理 $1$ 到 $m$ 的所有因子列表（使用筛法）。
• 线性筛预处理莫比乌斯函数 $\mu$。
• 按 $x$ 从小到大的顺序计算 $E[x]$（因为 $d<x$ 时 $E[d]$ 已经求出）：
  • 对每个 $x$，枚举其所有因子 $d$（包括 $x$）：
    • 令 $t = x/d$，枚举 $t$ 的所有因子 $k$，累加 $\mu(k) \cdot \lfloor m/(d k)\rfloor$ 得到 $c(x,d)$。
  • 将 $d<x$ 的贡献累加到 $\text{sum}$ 中：$\text{sum} = \sum_{d\mid x, d<x} E[d]\cdot c(x,d)$。
  • 计算 $E[x] = \frac{1 + \text{sum}/m}{1 - c(x,x)/m}$，所有运算在模 $10^9+7$ 下进行（除法用模逆）。
• 最后答案 $\text{ans} = 1 + \frac{1}{m}\sum_{x=1}^{m} E[x] \pmod{MOD}$。

5. 复杂度

• 预处理因子：$O(m\log m)$。
• 对于每个 $x$，需要枚举所有因子对 $(d,k)$，总次数为 $\sum_{x=1}^m \tau(x)^2$，其中 $\tau(x)$ 为 $x$ 的因子个数。对于 $m=10^5$，$\max\tau(x)=128$，平均约 $100$，总操作量约 $10^7$，完全可以接受。

6. 代码实现

#include <bits/stdc++.h>
using namespace std;

const int MOD = 1e9 + 7;
const int MAXM = 100005;

int mu[MAXM];
vector<int> factors[MAXM];
bool is_prime[MAXM];
vector<int> primes;

void sieve_mu(int n) {
    fill(is_prime, is_prime + n + 1, true);
    mu[1] = 1;
    for (int i = 2; i <= n; ++i) {
        if (is_prime[i]) {
            primes.push_back(i);
            mu[i] = -1;
        }
        for (int p : primes) {
            if (i * p > n) break;
            is_prime[i * p] = false;
            if (i % p == 0) {
                mu[i * p] = 0;
                break;
            } else {
                mu[i * p] = -mu[i];
            }
        }
    }
}

int modpow(int a, int e) {
    int res = 1;
    while (e) {
        if (e & 1) res = 1LL * res * a % MOD;
        a = 1LL * a * a % MOD;
        e >>= 1;
    }
    return res;
}

int inv(int x) {
    return modpow(x, MOD - 2);
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int m;
    cin >> m;
    sieve_mu(m);
    // 预处理因子
    for (int i = 1; i <= m; ++i) {
        for (int j = i; j <= m; j += i) {
            factors[j].push_back(i);
        }
    }
    vector<int> E(m + 1, 0);
    E[1] = 0;
    int inv_m = inv(m);
    for (int x = 2; x <= m; ++x) {
        long long sum = 0;
        int cxx = m / x; // c(x,x)
        for (int d : factors[x]) {
            int t = x / d;
            long long cnt = 0;
            for (int k : factors[t]) {
                cnt += 1LL * mu[k] * (m / (d * k));
            }
            cnt %= MOD;
            if (d == x) {
                // 跳过，用于分母
                cxx = cnt; // 实际上cxx已经等于floor(m/x)，但用公式算出来应该一样，这里保留原值也可
            } else {
                sum = (sum + 1LL * E[d] * cnt) % MOD;
            }
        }
        // 注意：cxx 实际就是 floor(m/x)，但为了避免重复计算，我们直接用 m/x
        cxx = m / x;
        long long denom = (1 - 1LL * cxx * inv_m % MOD + MOD) % MOD;
        long long num = (1 + 1LL * sum * inv_m % MOD) % MOD;
        E[x] = num * inv(denom) % MOD;
    }
    long long total = 0;
    for (int x = 1; x <= m; ++x) total = (total + E[x]) % MOD;
    int ans = (1 + 1LL * total * inv_m) % MOD;
    cout << ans << '\n';
    return 0;
}