CF1139D Steps to One

题目

打\(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)=\lfloor\frac{n}{k}\rfloor\times[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;
}