Codeforces Round #548 (Div. 2) D. Steps to One_codeforces round 548 (div. 2)-CSDN博客

本文链接：https://blog.csdn.net/qq_40538328/article/details/88858478

题目链接：Codeforces Round #548 (Div. 2) D. Steps to One

题目大意：一个数列，一开始是空的，每次往他最后一个位置随机的加上一个[1,m]范围内的数字，然后对当前数列所有的数字求gcd，如果gcd不为1，那么就继续添加，否则停止。问最后这个数列期望的长度。

思路：拿到题目首先明白一定是道期望dp。我们考虑dp[x]为此时序列gcd为x，到序列gcd为1时(结束)的长度。由全期望公式 $E[x]=\sum p[y]\cdot E[x\left | y].$ 得

$dp[x]=1+\frac{\sum_{1}^{m}dp[gcd(x,i)]}{m}.$

由此我们得到了一个 $O(n^{2})$ 的解法，显然复杂度时不够的 (1≤m≤100000).

因此我们考虑对每一个 $d=gcd(x,i)$ 来说 $d\left|x$ ，因此我们可以将枚举的复杂度降为 $O(logm)$ ,令f(d)表示 $\sum_{1}^{m}\left [ gcd(x,i)==d \right ]$

则上述状态转移可化为

$dp[x]=1+\frac{\sum dp[d]\cdot f(d) }{m}$ ,其中 $d\left|x$

(可对f(d)反演，稍后再说)

现在来计算f(d):

对于集合 $set=$ $\left \{ i \right|i\in \left [ 1,m\right ],gcd(x,i)=d}$ ,令 $\left\{\begin{matrix} & x=k_{x}\cdot d & \\ & i=k_{i}\cdot d \end{matrix}\right.$ ，即 $gcd(k_{x},k_{i})=1$ ，得到 $m^{\\'}=\frac{m}{d},x^{'}=\frac{m}{d}$ ,集合 $set^{'}$ 元素个数与原集合元素个数相同，该子问题转化为求在 $set^{'}$ 内有多少与x互质的元素。设n为 $x^{'}$ 的质因子数量，则子问题复杂度为 $O(2^{n}\cdot n)$ .

总时间复杂度为 $O(m\cdot log_{2}m\cdot 2^{n}\cdot n)$

代码(注释有很多细节)：

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <iostream>
#include <vector>
#include <map>
#include <unordered_map>
#include <set>
#include <queue>
typedef long long ll;
using namespace std;
#define INF 0x3f3f3f3f
const int mod=1e9+7;
const int maxn=1e5+10;
ll pow_mod(ll a,ll b,ll p){
    ll ret=1;
    while(b){
        if(b & 1) ret=(ret*a)%p;
        a=(a*a)%p;
        b>>=1;
    }
    return ret;
}
unordered_map<int,int>prime[maxn];//prime[j][i]表示j的素因子i的幂次(个数)
vector<int>fact[maxn];
bool isPrime[maxn];
ll dp[maxn];//dp[i]表示当前gcd和为i时到gcd和为1时的步长
ll m;
ll findCount(ll x,ll n){
	ll currCnt=m/x;
	vector<ll> arr;
	for(auto &it:prime[n])
		arr.push_back(it.first);
	
	ll sz=arr.size();
	ll lim=m/x;
	for(ll i=1;i<(1<<sz);i++){
		ll sign=1,prod=1;
		for(ll j=0;j<sz;j++)
			if(i&(1ll<<j)){
				prod*=arr[j];
				sign*=-1;
			}
		currCnt+=sign*(lim/prod);//sign表示容斥
	}
	return currCnt;
}
ll invm;
int main()
{
	memset(isPrime,true,sizeof(isPrime));
	//处理maxn内的所有素数及i的素因子
	for(ll i=2;i<maxn;i++){
		if(isPrime[i]){
			for(ll j=i;j<maxn;j+=i){
				ll curr=j,cnt=0;
				while(curr%i==0){
					curr/=i;
					cnt++;
				}
				prime[j][i]=cnt;
				isPrime[j]=false;
			}
		}
	}
	//处理maxn内的所有约数
	for(ll i=1;i<maxn;i++)
		for(ll j=2*i;j<maxn;j+=i)
			fact[j].push_back(i);
	cin>>m;
	invm=pow_mod(m,mod-2,mod);//m的逆元
	dp[1]=0;
	for(ll i=2;i<=m;i++){
		ll rhs=1;
		for(ll x:fact[i]){//遍历i的所的约数
			ll cnt=findCount(x,i/x);//gcd(i,j)==x的个数，j取1~m
			rhs+=(dp[x]*cnt)%mod*invm%mod;
			if(rhs>=mod)rhs-=mod;//常数优化
		}
		ll cnt=m/i;
		dp[i]=rhs*m%mod*pow_mod(m-cnt,mod-2,mod)%mod;//移相结果
	}
	ll res=1;
	for(ll i=1;i<=m;i++){
		res+=dp[i]*invm%mod;
		if(res>=mod)
			res-=mod;
	}
	cout<<res<<endl;
	return 0;
}

这道题的代码有很多值得学习的地方：

在线性筛的同时处理每个数的质因子。
因数分解的预处理。
取模的常数优化。
在计算f(n)的时候应用的容斥原理(原因是会有gcd(x,i)==x( $i$ $\geq$ $x$ ))),同时移项的时候也要注意！！！！！！(这一条其他题目也应用的非常多)

		ll cnt=m/i;
		dp[i]=rhs*m%mod*pow_mod(m-cnt,mod-2,mod)%mod;//移相结果

该部分代码为移项后的化简结果。

反演的方法日后再补。