1. 程式人生 > 實用技巧 >【CF1139D】Steps to One(期望+莫比烏斯反演)

【CF1139D】Steps to One(期望+莫比烏斯反演)

點此看題面

大致題意: 一個空數列,每次隨機加入一個\(1\sim m\)的元素,直至數列中所有元素\(gcd=1\)。求期望長度。

期望

關於期望有一個著名的公式:

\[E(X)=\sum_{i\ge1}P(X\ge i) \]

這裡的\(P(X\ge i)\)即為最終長度大於等於\(i\)的概率。

接下來的做法都要以這一公式為基礎。

推式子

考慮最終長度大於等於\(i\),等價於長度為\(i-1\)的序列\(gcd\not=1\)的概率。

總方案數是非常好計算的,就是\(m^{i-1}\),因此我們只要求出長度為\(i-1\)的序列\(gcd=1\)的方案數就可以了。

這一步可以使用容斥,容斥係數就是莫比烏斯函式

:(可參考此題:【BZOJ2440】[中山市選2011] 完全平方數

\[\sum_{j=1}^m\mu(j)\lfloor\frac mj\rfloor^i \]

整個式子就是:

\[1+\sum_{i\ge2}\frac{m^{i-1}-\sum_{j=1}^m\mu(j)\lfloor\frac mj\rfloor^{i-1}}{m^{i-1}} \]

觀察到式子中的\(i\)都以\(i-1\)的形式出現,方便起見我們將所有\(i\)\(1\)得到:

\[1+\sum_{i\ge1}\frac{m^i-\sum_{j=1}^m\mu(j)\lfloor\frac mj\rfloor^i}{m^i} \]

\(j=1\)\(\mu(j)\lfloor \frac mj\rfloor^i=m^i\),因此我們單獨提出這個式子,與前面的\(m^i\)抵消,得到:

\[1-\sum_{i\ge 1}\frac{\sum_{j=2}^m\mu(j)\lfloor\frac mj\rfloor^i}{m^i} \]

接下來就是喜聞樂見地調整列舉順序:

\[1-\sum_{j=2}^m\mu(j)\sum_{i\ge1}(\frac{\lfloor\frac mj\rfloor}{m})^i \]

由於\(\frac{\lfloor\frac mj\rfloor}{m}\)顯然是一個小於\(1\)的數,所以適用於公式\(\sum_{i\ge1} x^i=\frac{x}{1-x}\)

,因此得到:

\[1-\sum_{j=2}^m\mu(j)\frac{\frac{\lfloor\frac mj\rfloor}m}{1-\frac{\lfloor\frac mj\rfloor}m} \]

\[1-\sum_{j=2}^m\mu(j)\frac{\lfloor\frac mj\rfloor}{m-\lfloor\frac mj\rfloor} \]

這個式子可以除法分塊\(O(\sqrt n)\)求解。(但好像沒啥意義?除非出成多組資料)

程式碼

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 100000
#define X 1000000007
using namespace std;
int n,Pt,P[N+5],mu[N+5],s[N+5];
I int QP(RI x,RI y) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
I void Sieve()//線性篩預處理
{
	RI i,j;for(mu[1]=1,i=2;i<=n;++i)
		for(!P[i]&&(mu[P[++Pt]=i]=X-1),j=1;j<=Pt&&i*P[j]<=n;++j)
			if(P[i*P[j]]=1,i%P[j]) mu[i*P[j]]=X-mu[i];else break;
	for(i=1;i<=n;++i) s[i]=(s[i-1]+mu[i])%X;//預處理莫比烏斯函式字首和
}
int main()
{
	RI l,r,t=1;for(scanf("%d",&n),Sieve(),l=2;l<=n;l=r+1)//除法分塊
		r=n/(n/l),t=(1LL*(X-1)*(s[r]-s[l-1]+X)%X*(n/l)%X*QP(n-n/l,X-2)+t)%X;
	return printf("%d\n",t),0;
}