NC113552 Steps to One

题目地址：

https://ac.nowcoder.com/acm/problem/113552

基本题意：

给定一个空集合 $a$ ，和一个数 $m$ ;
每次等概率的从 $[1,m]$ 中选择一个数加入集合;
如果集合里所有数的 $gcd$ 不等于 $1$ 就一直加入直到 $gcd$ 为 $1$ 为止;
求结束时集合 $a$ 的期望大小。

基本思路：

有一说一我数论确实一窍不通，真的是鼓起了很大的勇气来挑战这道题，这里记录一下自己一步步分析的过程。

首先我们 $dp$ ，设 $dp[i]$ 记录集合 $gcd$ 为 $i$ 时的集合期望大小(也就是操作步数)，那么我们从 $gcd$ 小的往 $gcd$ 大的方向转移，就应该能得到这样的转移方程 : $dp[i] = 1 + \frac{1}{m} * \sum_{j=1}^{m}dp[gcd(i,j)]$ ，我们思考 $gcd(i,j)<=i$ 一定是恒成立的，而由于 $gcd(i,j) = i$ 时是自己转移自己，所以要分 $gcd(i,j) < i$ 的情况，和 $gcd(i,j) = i$ 两种来讨论。

我们先看 $gcd(i,j) < i$ 的情况，这种情况看上去是能直接转移了但是我们还要考虑复杂度的问题，如果不做处理复杂度是 $O(m^2)$ 肯定是不行的，那么我们考虑能不能预处理出 $[1,m]$ 中的数和 $i$ 进行 $gcd$ 运算时出现的每种可能答案与每种答案的出现次数。

首先这个每种可能答案很明显只会是 $i$ 的约数。那么上面的结论转化一下也就是说我们要找 $i$ 的每个约数会同样是 $[1,m]$ 中多少数的约数。

那么我们这样思考，假如数 $i$ 存在约数 $x$ ,那么在 $[1,m]$ 中的 $x,2*x,3*x,4*x...k*x <= m$ 这些数就也一定会存在 $x$ 这个约数，那么初步看上去 $i$ 的每个约数 $x$ 的出现次数就是 $floor(\frac{m}{x})$ 。

但是再思考一下就会发现有重复的情况，比如说 $i=4$ 的时候，那么算 $4$ 的时候会把实际 $gcd$ 出来应该是 $2$ 的部分重复的给计算了，那么考虑再容斥一下去重。这样我们求得了每种可能结果和对应出现次数，原转移方程就能做到如下的转换： $dp[i] = 1 + \frac{1}{m}*\sum_{j|i}dp[j] * cnt_{j}$ 其中 $j|i$ 表示 $j$ 能整除 $i$ 。

我们再考虑一下之前所说的 $gcd(i,j) = i$ 这个情况，这个情况很明显，当且仅当 $i|j$ 时才会出现而 $j\epsilon[1,m]$ ，根据上面推的那部分，同原理能知道这时候还要加 $floor(\frac{m}{i})$ 个 $dp[i]$ 。

最后我们综上结合两种情况整理我们的转移方程：
我们设 $x = \sum_{j|i\&\&gcd(i,j)<i}dp[j] * cnt_{j}$ ;
那么就能得到 $dp[i] = 1 + \frac{1}{m}*(x + dp[i] * floor(\frac{m}{i}))$
我们近一步化简将 $dp[i]$ 移到一边最终就有： $dp[i] = \frac{m + x}{m - floor(\frac{m}{i})}$ 。
那么我们终于就得到了这个真正的 $dp$ 转移方程！
最终答案就是： $ans = \frac{1}{m} * \sum_{i = 1}^{m} (dp[i] + 1)$ ;
就是枚举第一次放进去的数，然后第一次放每个数的概率同样是 $\frac{1}{m}$ 。
ps.好像还有莫比乌斯反演的做法，实在不会了QAQ。

参考代码：

#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>
using namespace std;
#define IO std::ios::sync_with_stdio(false)
#define int long long
#define rep(i, l, r) for (int i = l; i <= r; i++)
#define per(i, l, r) for (int i = l; i >= r; i--)
#define mset(s, _) memset(s, _, sizeof(s))
#define pb push_back
#define pii pair <int, int>
#define mp(a, b) make_pair(a, b)
#define INF (int)1e18

inline int read() {
  int x = 0, neg = 1; char op = getchar();
  while (!isdigit(op)) { if (op == '-') neg = -1; op = getchar(); }
  while (isdigit(op)) { x = 10 * x + op - '0'; op = getchar(); }
  return neg * x;
}
inline void print(int x) {
  if (x < 0) { putchar('-'); x = -x; }
  if (x >= 10) print(x / 10);
  putchar(x % 10 + '0');
}

const int maxn = 1e5 + 10;
const int p = 1e9 + 7;
int qsm(int x,int n) {
  int res = 1;
  while (n) {
    if (n & 1) res = res * x % p;
    x = x * x % p;
    n >>= 1;
  }
  return res;
}
int m;
vector<int> fac[maxn],cnt[maxn];
int dp[maxn];
signed main() {
  IO;
  cin >> m;
  //预处理过程;
  rep(i, 1, m) {
    for (int j = i; j <= m; j += i) //找到[1,m]中每个数的因数;
      fac[j].push_back(i), cnt[j].push_back(0);
  }
  for (int i = 1; i <= m; i++) {
    for (int j = (int) fac[i].size() - 1; j >= 0; j--) {
      cnt[i][j] = m / fac[i][j];
      for (int k = j + 1; k < fac[i].size(); k++) { //容斥去重;
        if (fac[i][k] % fac[i][j] == 0) cnt[i][j] -= cnt[i][k];
      }
    }
  }

  dp[1] = 0;
  for (int i = 2; i <= m; i++) {
    int tmp = 0;
    for (int j = 0; j < fac[i].size(); j++) {
      if (fac[i][j] != i) tmp = (tmp + dp[fac[i][j]] * cnt[i][j] % p) % p;
    }
    int q1 = (m + tmp) % p;
    int k = m / i;
    int q2 = (m - k + p) % p;
    dp[i] = q1 * qsm(q2, p - 2) % p;
  }
  int ans = 0;
  rep(i, 1, m) {
    ans = (ans + dp[i] + 1) % p;
  }
  ans = ans * qsm(m, p - 2) % p;
  cout << ans % p << '\n';
  return 0;
}