Codeforces Round #548 (Div. 2) D. Steps to One

链接

https://codeforces.com/contest/1139/problem/D

题解

我好菜啊
本来想回紫，结果在掉分的路上渐行渐远
首先可以直接写出期望$dp$方程
$f_i$表示我写的序列的最大公约数已经等于$i$的情况下，写到结束为止的期望长度
那么$$f_i=1+\frac{1}{m}\sum _{j=1}^n{f}_{gcd(i,j)}$$
移项整理得
$$\frac{m-\lfloor \frac{m}{i}\rfloor }{m}{f}_i=1+\frac{1}{m}\sum _{j=1}^n{f}_{gcd(i,j)}[(i,j)̸=n]$$
直接$dp$的复杂度是$O(m^2)$
会超时
那我就统计一下$i$的每种约数对$i$的贡献
即先枚举$x$（作为约数），再去枚举$x$的倍数$n$
现在就是要求一个
$$\sum _{i=1}^m[(n,i)=x]$$
反演得$$\sum _{d|n}\mu (d)\lfloor \frac{m}{d}\rfloor$$
预处理$n$的约数就行了

代码

#include <bits/stdc++.h>
#define eps 1e-8
#define iinf 0x3f3f3f3f
#define linf (1ll<<60)
#define maxn 200010
#define base 2000
#define cl(x) memset(x,0,sizeof(x))
#define mod 1000000007ll
using namespace std;
typedef long long ll;
ll read(ll x=0)
{
	ll c, f=1;
	for(c=getchar();!isdigit(c);c=getchar())if(c=='-')f=-f;
	for(;isdigit(c);c=getchar())x=x*10+c-48;
	return f*x;
}
ll fastpow(ll a, ll b)
{
	ll ans=1, t=a;
	for(;b;b>>=1,t=t*t%mod)if(b&1)ans=ans*t%mod;
	return ans;
}
struct fs
{
	ll a, b;
	fs(ll a, ll b):a(a),b(b){}
	fs(){}
}f[maxn], ans;
fs operator+(fs x, fs y)
{
	return fs((x.a*y.b+x.b*y.a)%mod,x.b*y.b%mod);
}
fs operator*(fs x, fs y)
{
	return fs((x.a*y.a)%mod,(x.b*y.b)%mod);
}
ll M, phi[maxn], mark[maxn], prime[maxn], mu[maxn];
vector<ll> ys[maxn];
void init()
{
	ll i, j;
	phi[1]=1;
	mu[1]=1;
	for(i=2;i<maxn;i++)
	{
		if(!mark[i])
		{
			prime[++*prime]=i;
			phi[i]=i-1;
			mu[i]=-1;
		}
		for(j=1;j<=*prime and i*prime[j]<maxn;j++)
		{
			mark[i*prime[j]]=1;
			if(i%prime[j]==0)
			{
				phi[i*prime[j]]=phi[i]*prime[j];
				break;
			}
			else phi[i*prime[j]]=phi[i]*phi[prime[j]], mu[i*prime[j]]=-mu[i];
		}
	}
	for(i=1;i<maxn;i++)for(j=i;j<maxn;j+=i)
	{
		ys[j].push_back(i);
	}
}
ll calc(ll m, ll n)
{
	ll ans=0, i, d;
	for(i=0;i<ys[n].size();i++)
	{
		d=ys[n][i];
		ans=(ans+mu[d]*(m/d))%mod;
	}
	return (ans%mod+mod)%mod;
}
int main()
{
	ll i, j, k, cnt;
	M=read();
	init();
	for(i=1;i<=M;i++)f[i]=fs(0ll,1ll);
	fs ans(0ll,1ll);
	for(i=1;i<=M;i++)
	{
		if(i>1)
		{
			f[i]=fs(1ll,M)*f[i]+fs(1ll,1ll);
			f[i]=f[i]*fs(M,M-M/i);
		}
		for(j=i+i;j<=M;j+=i)
		{
			cnt=calc(M/i,j/i);
			f[j]=f[j]+fs(cnt,1ll)*f[i];
		}
	}
	for(i=1;i<=M;i++)
	{
		ans=ans+(f[i]+fs(1ll,1ll))*fs(1ll,M);
	}
	cout<<ans.a*fastpow(ans.b,mod-2)%mod;
	return 0;
}