反演套 DP 的好题（不用反演貌似也能做
Description
Vivek initially has an empty array $a$ and some integer constant $m$.
He performs the following algorithm:
- Select a random integer $x$ uniformly in range from $1$ to $m$ and append it to the end of $a$.
- Compute the greatest common divisor of integers in $a$.
- In case it equals to $1$, break
- Otherwise, return to step $1$.
Find the expected length of $a$. It can be shown that it can be represented as $\frac PQ$ where $P$ and $Q$ are coprime integers and $Q\neq 0\pmod{10^9+7}$. Print the value of $P\cdot Q^{-1}\pmod{10^9+7}$.
Input
The first and only line contains a single integer $m$($1\le m\le 100000$).
Output
Print a single integer — the expected length of the array $a$ written as $P\cdot Q^{-1}\pmod{10^9+7}$.
Examples
input
1output
1input
2output
2input
4output
333333338Note
In the first example, since Vivek can choose only integers from $1$ to $1$, he will have $a=[1]$ after the first append operation, and after that quit the algorithm. Hence the length of $a$ is always $1$, so its expected value is $1$ as well.
In the second example, Vivek each time will append either $1$ or $2$, so after finishing the algorithm he will end up having some number of $2$'s (possibly zero), and a single $1$ in the end. The expected length of the list is $1⋅\frac 12+2⋅\frac 1{2^2}+3⋅\frac 1{2^3}+\dots =2$.
题意
每一步在 $1\sim m$ 中任选一个整数，问期望多少步后选出的数的最大公约数是 $1$。答案对 $1\ 000\ 000\ 007$ 取模。
题解
因为每一步不会让已经选了的元素的 $\gcd$ 和变大，因此认为是一个除自环外的有向无环图。对于自环我们很好处理，所以把它看成是一道期望 DP。
令 $f[i]$ 表示当前的 $\gcd$ 和为 $i$，到 $\gcd$ 和为 $1$ 的状态的期望步数。因此把状态转移方程写出来 $$ f[i]=1+\frac{\sum_{j=1}^m f[\gcd(i,j)]}m $$ 这样的转移是 $O(m^2)$ 的。但是我们发现，对于很多 $j$，$\gcd(i,j)$ 都是相等的，因此我们把这样的数整合到一起。
令 $F(n)$ 表示 $1\sim m$ 中有多少个数 $i$ 满足 $\gcd(x,i)=n$，其中视 $x$ 为常数。
则计算 $f[i]$ 就转化为了 $$ f[i]=1+\frac{\sum_{d|i}f[d]\times F(d)}n,x=i $$ 这样差不多就把枚举优化到了 $\log n$ 的 $d|i$。
考虑怎么计算 $F(n)$ $$ F(n)=\sum_{i=1}^m[\gcd(x,i)=n] $$ 令 $G(n)=\sum_{n|d}F(d)$，则 $$ \begin{aligned} G(n)&=\sum_{n|d}F(d)\ &=\sum_{n|d}\sum_{i=1}^m[\gcd(x,i)=d]\ &=\sum_{i=1}^m[n|\gcd(x,i)]\ ?&=\sum_{i=1}^{\left\lfloor\frac mn\right\rfloor}\left[1|\gcd\left(\frac xn,i\right)\right] \end{aligned} $$ 实际上这样是有问题的，因为（在后面）无法保证 $n|x$，此时 $G(n)$ 就一定为 $0$ 了。
我们再多化一步： $$ \begin{aligned} G(n)&=\sum_{i=1}^{\left\lfloor\frac mn\right\rfloor}\left[1|\gcd\left(\frac xn,i\right)\right][n|x]\ &=\left\lfloor\frac mn\right\rfloor\cdot[n|x] \end{aligned} $$ 根据 $G(n)=\sum_{n|d}F(d)$，我们反演到 $F$，得 $$ \begin{aligned} F(n)&=\sum_{n|d}\mu\left(\frac dn\right)G(d)\ &=\sum_{i=1}^{\left\lfloor\frac mn\right\rfloor}\mu(i)G(ni)\ &=\sum_{i=1}^{\left\lfloor\frac mn\right\rfloor}\mu(i)\left\lfloor\frac{m}{ni}\right\rfloor[(ni)|x] \end{aligned} $$ 我们发现后面的布尔表达式可以当作条件。原本的条件本来就是 $\to +\infty$ 的，只不过超过了 $\left\lfloor\frac mn\right\rfloor$ 没有意义。因此直接把条件换成 $[(ni)|x]$ 即可。又因为 $n|x$ 在上面的枚举过程中是成立的，同时可以转化为 $\left[i|\frac xn\right]$。 $$ F(n)=\sum_{i|\frac xn}\mu(i)\left\lfloor\frac{m}{ni}\right\rfloor $$ 这样的一次枚举是 $O\left(d\left(\frac xn\right)\right)$ 的，由于 $1\sim m$ 的约数个数和均摊是 $O(\log m)$ 的，其中最多的有 $128$ 个约数，但是这样的数肯定不是很多，并且其中很多被枚举到的数都是质因数，迭代一下并不会造成很大的复杂度。
然后我们需要再把状态转移方程稍微转化一下，把 $f[i]$ 移到左边 $$ \begin{aligned} f[i]&=1+\frac{\sum_{d|i,d<i}f[d]\times F(d)+f[i]\times F(i)}n,x=i\ \frac{n-F(i)}{n}\cdot f[i]&=1+\frac{\sum_{d|i,d<i}f[d]\times F(d)}n,x=i\ f[i]&=\frac{n+\sum_{d|i,d<i}f[d]\times F(d)}{n-F(i)},x=i \end{aligned} $$ 就得到了真正的转移方程。
时间复杂度 $O(m\log^2 m)$。
Code：
1#include<cstdio> 2#include<cstring> 3#include<vector> 4#define ll long long 5#define p 1000000007 6using std::vector; 7vector<int> v[100100];//约数用 vector 存一下，每次 √m 枚举不是很稳 8ll qpow(ll x,ll y) 9{ 10 ll ans=1; 11 while(y) 12 { 13 if(y&1) 14 ans=ans*x%p; 15 x=x*x%p; 16 y>>=1; 17 } 18 return ans; 19} 20bool is[100100]; 21int pri[100100],mu[100100],cnt=0; 22ll f[100100]; 23int n; 24int calc(int x,int y)//1~n 中 gcd(x,i)=y 的数的个数 25{ 26 int g=x/y,ans=0; 27 for(int i=0;i<v[g].size();++i) 28 ans+=mu[v[g][i]]*(n/v[g][i]/y); 29 return ans; 30} 31int main() 32{ 33 scanf("%d",&n); 34 f[1]=1; 35 mu[1]=1; 36 for(int i=2;i<=n;++i) 37 { 38 if(!is[i]) 39 { 40 pri[++cnt]=i; 41 mu[i]=-1; 42 } 43 for(int j=1;j<=cnt&&i*pri[j]<=n;++j) 44 { 45 is[i*pri[j]]=1; 46 if(i%pri[j]) 47 mu[i*pri[j]]=-mu[i]; 48 else 49 { 50 mu[i*pri[j]]=0; 51 break; 52 } 53 } 54 } 55 for(int i=1;i<=n;++i) 56 for(int j=i;j<=n;j+=i) 57 v[j].push_back(i); 58 ll ans=1,inv=qpow(n,p-2); 59 for(int i=2;i<=n;++i) 60 { 61 for(int j=0;j<v[i].size()-1;++j) 62 f[i]=(f[i]+calc(i,v[i][j])*f[v[i][j]]%p)%p; 63 f[i]=(f[i]*inv+1)%p; 64 ll g=n-calc(i,i); 65 f[i]=f[i]*n%p*qpow(g,p-2)%p; 66 ans=(ans+f[i])%p; 67 } 68 printf("%lld\n",ans*qpow(n,p-2)%p); 69 return 0; 70}