题意：重复从 $1,n$ 中随机加一个数到序列末，直至序列的 $\gcd$ 为 $1$ 为止，序列初始为空，求序列的期望长度。

假设序列的长度为 $k$，那么不难算出每个数都是 $d$ 的倍数的概率 $g_d$

\[g_d=\left(\frac{\lfloor n/d\rfloor}n\right)^k\]

不妨设 $f_n$ 表示长度为 $k$，整个序列的 $\gcd$ 为 $n$ 的概率，不难看出

\[g_d=\sum_{d\mid n}f_n\]

也就是 $g=1\ast f$，莫反一下，$f=\mu\ast g$

\[f_n=\sum_{n\mid d}\mu(d/n)g_d\]

我们只需要计算 $f_1$，也就是 $\gcd=1$ 的概率

\[f_1=\sum_{d}\mu(d)g_d=\sum_{i=1}^n\mu(i)\left(\frac{\lfloor n/i\rfloor}n\right)^k\]

这里只枚举到 $n$ 是因为 $i>n$ 时 $g_i=0$

因此，设长度为 $k$ 时序列的 $\gcd=1$ 的概率为 $t(k)$，那么

\[t(k)=\sum_{i=1}^n\mu(i)\left(\frac{\lfloor n/i\rfloor}n\right)^k\]

令 $T$ 为停时，那么

\[\begin{aligned} \mathbb E[T]&=\sum_{k\ge1}k(t(k-1)-t(k))\\ &=\sum_{k\ge0}t(k)\\ &=\sum_{k\ge0}\sum_{i=1}^n\mu(i)\left(\frac{\lfloor n/i\rfloor}n\right)^k\\ &=\sum_{i=1}^n\mu(i)\sum_{k\ge0}\left(\frac{\lfloor n/i\rfloor}n\right)^k\\ &=\sum_{i=1}^n\mu(i)\frac1{1-\frac{\lfloor n/i\rfloor}n}\\ &=\sum_{i=1}^n\mu(i)\frac n{n-\lfloor n/i\rfloor}\\ \end{aligned}\]

按 $\lfloor n/i\rfloor$ 整除分块即可，前面 $\mu$ 前缀和采用杜教筛，复杂度 $\mathcal O(n^{2/3})$

但是模数 $P\le10^{12}$ 是来恶心大家的吗？

#include<bits/stdc++.h>
#define siz(x) int((x).size())
#define all(x) std::begin(x),std::end(x)
#define fi first
#define se second
#define int loli
using namespace std;
using unt=unsigned;
using loli=long long;
using lolu=unsigned long long;
using pii=pair<int,int>;
mt19937_64 rng(random_device{}());
constexpr int N=1e7+7;
bool pr[N];
int mu[N],n,P;
loli mv[N];
basic_string<int>pt;
unordered_map<loli,loli>mw;
loli calc_mu(loli n){
	if(n<N)return mv[n];
	if(mw.count(n))return mw[n];
	loli s=1;
	for(loli l=2,r;l<=n;l=r+1){
		r=n/(n/l);
		s-=(r-l+1)*calc_mu(n/l);
	}
	return mw[n]=s;
}
using venti=__int128_t;
loli qpow(venti x,venti y){
	venti z=1;
	for(;y;y>>=1,x=x*x%P)
		if(y&1)z=z*x%P;
	return z;
}
signed main(){
//	freopen(".in","r",stdin);
//	freopen(".out","w",stdout);
	ios::sync_with_stdio(false);cin.tie(nullptr);
	mu[1]=mv[1]=1;
	for(int i=2;i<N;i++){
		if(!pr[i])pt+=i,mu[i]=-1;
		for(int j:pt){
			if(i*j>=N)break;
			pr[i*j]=true;
			if(i%j==0){
				mu[i*j]=0;
				break;
			}
			mu[i*j]=-mu[i];
		}
		mv[i]=mv[i-1]+mu[i];
	}
	cin>>n>>P;
	venti ans=1;
	for(int l=2,r;l<=n;l=r+1){
		r=n/(n/l);
		venti t1=calc_mu(r)-calc_mu(l-1);
		venti t2=venti(n/l)*qpow((n-(n/l)),P-2)%P;
		ans=(ans-t1*t2%P+P)%P;
	}
	cout<<loli(ans);
	return 0;
}