打(div2)的时候被这道题锤自闭了
确实是非常神仙的的一道题
我们发现不能和结尾的数互质，结尾的数也不会和前面的数互质，所以我们发现整个序列都是不互质的
那么我们可以考虑把这个关系抽象成一张拓扑图，每一个数(x)向其约数(k)连边，表示我们通过在末尾添加一个数可以使得整个序列的(gcd)从(x)变成(k)
于是我们得到了一张(DAG)，我们再搞一个虚点，向所有数连边，表示起点
每一条边都有一个走的概率，我们已经把问题转化成了求从起点到(1)的期望路径长度
从起点连向各点的边概率显然是(frac{1}{m})
从(x)连向其约数(k)的边的概率显然是(frac{sum_{i=1}^m[(x,i)=k]}{m})
现在的问题转化成了如何求
[f(x,k)=sum_{i=1}^m[(x,i)=k]
]
考虑反演
设
[F(x,k)=sum_{i=1}^m[k|(x,i)]
]
显然存在
[F(x,k)=sum_{k|d}f(x,d)=lfloorfrac{n}{k}
floor imes[k|x]
]
直接反演可得
[f(x,k)=sum_{k|d}mu(frac{d}{k})F(x,d)
]
注意这里的(d)仅能取(x)的约数
我们求(f(x,k))直接暴力就好了，复杂度基于(sum_{i=1}^md^2(i))，不是很会算，大概(O(nlog^2n))吧
还有一个问题，就是每一个点自己能像自己转移，看起来没有什么办法直接(dp)了呀
照样做就好了
设(dp_x)表示从(1)到(x)的期望路径长度
显然
[dp_x=frac{f(x,x)}{m}dp_x+1+sum_{k|x,k!=x}frac{f(x,k)}{m}dp_k
]
移项可得
[(1-frac{f(x,x)}{m})dp_x=1+sum_{k|x,k!=x}frac{f(x,k)}{m}dp_k
]
求出右边之后直接解方程就好
代码
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<vector>
#define re register
#define LL long long
const int maxn=200005;
const int M=2500005;
const LL mod=1e9+7;
int n,num;
std::vector<int> v[maxn];
int f[maxn],mu[maxn],p[maxn>>1],head[maxn],vis[maxn];
LL a[maxn],g[maxn],dp[maxn];
struct E{int v,nxt;LL w;}e[M];
inline void add(int x,int y,LL w) {
e[++num].v=y;e[num].nxt=head[x];
head[x]=num;e[num].w=w;
}
inline LL ksm(LL a,LL b) {LL S=1;while(b) {if(b&1) S=S*a%mod;b>>=1;a=a*a%mod;}return S;}
inline LL dfs(int x) {
if(vis[x]) return dp[x];
vis[x]=1;
for(re int i=head[x];i;i=e[i].nxt)
dp[x]=(dp[x]+e[i].w*(dfs(e[i].v)+1ll)%mod)%mod;
if(a[x]) {
dp[x]=(dp[x]+a[x])%mod;
dp[x]=dp[x]*ksm((1-a[x]+mod)%mod,mod-2);
dp[x]%=mod;
}
return dp[x];
}
int main() {
scanf("%d",&n);
for(re int i=1;i<=n;i++)
for(re int j=i;j<=n;j+=i) v[j].push_back(i);
f[1]=mu[1]=1;
for(re int i=2;i<=n;i++) {
if(!f[i]) p[++p[0]]=i,mu[i]=-1;
for(re int j=1;j<=p[0]&&p[j]*i<=n;j++) {
f[p[j]*i]=1;
if(i%p[j]==0) break;
mu[p[j]*i]=-1*mu[i];
}
}
LL inv=ksm(n,mod-2);
for(re int i=1;i<=n;i++) {
for(re int j=0;j<v[i].size();j++) {
int now=0,x=v[i][j];
if(v[i][j]==i) continue;
for(re int k=j;k<v[i].size();k++) {
int t=v[i][k];
if(t%x==0) now+=mu[t/x]*(n/t);
}
add(i,x,now*inv%mod);
}
a[i]=n/i;
a[i]=(a[i]*inv)%mod;
}
for(re int i=1;i<=n;i++) add(0,i,inv);
vis[1]=0;
std::cout<<dfs(0);
return 0;
}