聽起來很 nb，很有名但比較難學的一個算法類型。然而確實很 nb。
我竟然在學 ymx 一年半前就學過的東西。
1. 反演的本質與第一反演公式
1.1. 什么是反演
反演是通過用 \(f\) 表示 \(g\) 的方法求出如何用 \(g\) 表示 \(f\)。
如果我們已知 \(g(n)\) 以及它如何用 \(f\) 表示：\(g(n)=\sum_{i=1}^nc_{n,i}f(i)\)，從而反推出用 \(g\) 表示 \(f\) 的方法：\(f(n)=\sum_{i=1}^nd_{n,i}g(i)\)，那么就算 \(f\) 無法直接求，也可以通過 \(g(i)\) 的值推出 \(f(n)\)。
上面二式互為反演公式。
1.2. 第一反演公式
這里直接給出定理：如果 \(n\) 次多項式 \(p\)，\(q\) 分別滿足以下公式：\(p_n(x)=\sum_{i=1}^n c_{n,i}q_i(x)\)，\(q_n(x)=\sum_{i=1}^n d_{n,i} p_i(x)\)，且 \(c_{j,j}\neq0,d_{j,j}\neq 0\)，那么
簡單來說，如果將 \(c,d\) 寫成矩陣的形式 \(\mathbf{C},\mathbf{D}\)，那么它顯然是下三角矩陣，同時有 \(\mathbf{CD}=\mathbf{I}\)。因為如果將 \(p,q\) 寫成列向量 \(\mathbf{p},\mathbf{q}\)，那么條件中的兩個公式可以表示為 \(\mathbf{p}=\mathbf{Cq}\)，\(\mathbf{q}=\mathbf{Dp}\)，即 \(\mathbf{p}=\mathbf{CDp}\)，得證。由於 \(\mathbf{C}\) 與 \(\mathbf{D}\) 互逆，所以它們的行列式不為 \(0\)，而因為它們是下三角矩陣，因此一定有對角線元素不為 \(0\)，這也就解釋了為什么 \(c_{j,j}\neq 0\)，\(d_{j,j}\neq 0\)。
通俗地講，只要我們找到任何一組符合條件的 \(p,q\) 與 \(c,d\)，就可以將 \(p,q\) 推廣至任意 \(n\) 次多項式，因為 \(\mathbf{C}\) 與 \(\mathbf{D}\) 永遠是互逆的，不隨 \(p,q\) 的改變而改變。這就意味着為了證明某個反演公式成立，我們只需考慮一種方便證明的特殊情況即可。
2. 二項式反演
2.1. 引入
在玩二項式定理的時候，如果你嘗試把 \(x\) 變成 \(1+(x-1)\)，那么可能會有新的收獲：
耶？和反演公式的形式長得很像誒！那么試試能不能用 \(x\) 表示 \(x-1\)：
很好！如果令 \(p_i=x^i\)，\(q_i=(x-1)^i\)，\(c_{i,j}=\dbinom{i}{j}\)，\(d_{i,j}=\dbinom{i}{j}(-1)^{i-j}\)，那么我們就成功找到了一組 \(p,q,c,d\)！
綜上，我們有：
其中 \(f_n\) 表示至多 \(n\) 個，\(g_n\) 表示恰好 \(n\) 個。這就是二項式反演，上面的引入就是它的證明。
如果將矩陣轉置，那么還可以得到另一種表示方法：
其中 \(f_n\) 表示至少 \(n\) 個，\(g_n\) 表示恰好 \(n\) 個。
2.2. 直觀理解
實際上二項式反演的本質是容斥原理。\(f_i\) 表示至少 / 至多 \(i\) 個的方案數，而 \(g_i\) 表示恰好 \(i\) 個的方案數。
一般情況下，使用容斥原理時，如果固定 \(i\)，那么任選 \(i\) 個元素的貢獻不一定相等。這就意味着我們必須 \(2^n\) 枚舉每一種情況。而相反，如果任選 \(i\) 個數時貢獻只與 \(i\) 有關而與所選擇的數無關，由於固定了 \(i\) 故也即貢獻相等，就可以簡化為二項式反演。這是使用二項式反演的情況之一，另一種情況見應用部分。
2.3. 應用
在 OI 中，二項式反演的應用通常與動態規划相結合（當然也不乏直觀解釋中的情況，如下方的三類問題）：DP 求出至少 / 至多選擇 \(i\) 個元素時的總貢獻是簡單的，那么我們就可以通過二項式反演，得到恰好選擇 \(i\) 個元素時的總貢獻。第四節例題部分很好地體現了這一點。
綜上，我們總結出二項式定理的核心：進行至少 / 至多和恰好的轉換。
2.3.1. 錯排問題
\(n\) 個人排列，求每個人都恰好站錯的方案數。
設 \(g_n\) 表示 \(n\) 個人的錯排方案數，而 \(f_n\) 表示所有排列方案數，即 \(n!\)，那么顯然有 \(f_n=\sum_{i=0}^n\dbinom{n}{i}g_i\)，\(i\) 的意義是站錯的人的個數。從而反推出 \(g_n=\sum_{i=0}^n(-1)^{n-i}\dbinom{n}{i}f_i=\sum_{i=0}^n(-1)^{n-i}\dfrac{n!}{(n-i)!}\)（組合數分母上的 \(i!\) 與 \(f_i\) 抵消了）。
2.3.2. 染色問題 I（一維）
\(n\) 個格子，\(k\) 種顏色，相鄰格子顏色不同，每個顏色至少出現一次，求方案數。
設 \(g_k\) 表示恰好出現 \(k\) 個顏色的方案數，\(f_k\) 表示至多出現 \(k\) 個顏色的方案數，即 \(k\times (k-1)^{n-1}\)。顯然有 \(f_k=\sum_{i=0}^k \dbinom{k}{i}g_i\)，從而反推出 \(g_k=\sum_{i=0}^k(-1)^{k-i}\dbinom{k}{i}f_i\)。顯然可以 \(\mathcal{O}(k\log n)\) 計算。
2.3.3. 染色問題 II（二維）
\(n\times m\) 個格子，\(1\) 種顏色，每行和每列至少有一個格子被塗色，求方案數。
設 \(g_{i,j}\) 表示恰好有 \(i\) 行 \(j\) 列被塗色，\(f_{i,j}\) 表示至多有 \(i\) 行 \(j\) 列被塗色的方案數，即 \(2^{ij}\)。那么對 \(i\) 和 \(j\) 這一維分別二項式反演，有 \(f_{n,m}=\sum_{i=0}^n\dbinom{n}{i}\sum_{j=0}^m \dbinom{m}{j}g_{i,j}\)，從而反推出 \(g_{n,m}=\sum_{i=0}^n\sum_{j=0}^m (-1)^{(n-i)+(m-j)}\dbinom{n}{i}\dbinom{m}{j}f_{i,j}\)，時間復雜度 \(\mathcal{O}(n^2)\)。
2.3.4. 染色問題 III（三維）
\(n\times m\) 個格子，\(k\) 種顏色，每行和每列至少一個格子被染色，每個顏色至少出現一次，格子可不被塗色，求方案數。
請讀者結合染色問題 I. 與 II. 自行思考。
2.4. 參考文章
2.5. 例題
*I. P6478 [NOI Online #2 提高組] 游戲
很不錯的題目，可以用來入門二項式反演。
設 \(f_x\) 為欽定 \(x\) 對祖先關系的情況總數（即祖先關系至少有 \(x\) 對），\(g_x\) 為答案，那么根據二項式反演有 \(g_x=\sum_{i=x}^{n}(-1)^{i-x}\dbinom{i}{x}f_i\)。
\(f_x\) 可以樹形背包求，別忘了最后乘上 \((m-x)!\)，因為剩下 \(m-x\) 對點對順序可以任意。時間復雜度 \(\mathcal{O}(n^2)\)。
const int N=5e3+5;
ll n,sz[N],a[N],f[N][N],C[N][N];
char s[N];
vector <int> e[N];
void dfs(int id,int fa){
f[id][0]=1;
for(int it:e[id]){
if(it==fa)continue;
dfs(it,id),mem(f[0],0,N);
for(int i=0;i<=sz[id]/2;i++)
for(int j=0;j<=sz[it]/2;j++)
f[0][i+j]=(f[0][i+j]+f[id][i]*f[it][j])%mod;
sz[id]+=sz[it],a[id]+=a[it],cpy(f[id],f[0],N);
}
for(int i=min(a[id],sz[id]-a[id])+1;i;i--)
f[id][i]=(f[id][i]+f[id][i-1]*((s[id]=='1'?a[id]:sz[id]-a[id])-i+1))%mod;
sz[id]++,a[id]+=s[id]=='0';
}
int main(){
scanf("%d%s",&n,s+1);
for(int i=1;i<n;i++){
int u,v; cin>>u>>v;
e[u].pb(v),e[v].pb(u);
} dfs(1,0);
for(int i=0;i<=n/2;i++)
for(int j=2;j<=n/2-i;j++)
f[1][i]=f[1][i]*j%mod;
for(int i=0;i<=n;i++)
for(int j=0;j<=i;j++)
C[i][j]=(j==0||j==i?1:C[i-1][j-1]+C[i-1][j])%mod;
for(int i=0;i<=n/2;i++){
ll ans=0;
for(int j=i;j<=n;j++)
ans=(ans+((j-i)&1?mod-1:1)*C[j][i]%mod*f[1][j])%mod;
cout<<ans<<endl;
}
return 0;
}
*II. P4859 已經沒有什么好害怕的了
和上一題極度類似。DP 方法幾乎一樣，最后也都要乘上階乘。
實際上，上一題就是把這一題搬到了樹上。
ll ksm(ll a,ll b){
ll s=1;
while(b){
if(b&1)s=s*a%mod;
a=a*a%mod,b>>=1;
}
return s;
}
ll inv(ll x){return ksm(x%mod,mod-2);}
ll fc[N],ifc[N];
ll binom(int n,int m){return fc[n]*ifc[m]%mod*ifc[n-m]%mod;}
ll n,k,c,ans;
ll a[N],b[N],f[N];
int main(){
cin>>n>>k,f[0]=1,c=n+k>>1;
if((n&1)^(k&1))puts("0"),exit(0);
for(int i=0;i<=n;i++)fc[i]=i?fc[i-1]*i%mod:1,ifc[i]=inv(fc[i]);
for(int i=1;i<=n;i++)cin>>a[i]; sort(a+1,a+n+1);
for(int i=1;i<=n;i++)cin>>b[i];
for(int i=1;i<=n;i++){
int cnt=0;
for(int j=1;j<=n;j++)cnt+=b[j]<a[i];
for(int j=cnt;j;j--)f[j]=(f[j]+f[j-1]*(cnt-j+1))%mod;
}
for(int i=0;i<=n;i++)f[i]=f[i]*fc[n-i]%mod;
for(int i=c;i<=n;i++)ans=(ans+(i-c&1?mod-1:1)*binom(i,c)%mod*f[i])%mod;
cout<<ans<<endl;
return 0;
}
*III. CF1228E Another Filling the Grid
以前不會做，現在也不會。
設 \(g_{i,j}\) 表示恰好有 \(i\) 行 \(j\) 列的最小值為 \(1\)，\(f_{i,j}\) 表示最多有 \(i\) 行 \(j\) 列的最小值為 \(1\)，即 \((k-1)^{(n-i)(n-j)}\times k^{n^2-(n-i)(n-j)}\)，那么有 \(f_{n,n}=\sum_{i=0}^n\sum_{j=0}^n\dbinom{n}{i}\dbinom{n}{j}g_{i,j}\)，二項式反演得到 \(g_{n,n}=\sum_{i=0}^n\sum_{j=0}^n(-1)^{n-i+n-j}\dbinom{n}{i}\dbinom{n}{j}f_{i,j}\)。可以 \(\mathcal{O}(n^2\log k)\) 計算，時間復雜度足夠優秀。當然，細節實現可以做到平方甚至線性對數。
l ksm(ll a,ll b){
ll s=1;
while(b){
if(b&1)s=s*a%mod;
a=a*a%mod,b>>=1;
}
return s;
}
ll inv(ll x){return ksm(x%mod,mod-2);}
ll fc[N],ifc[N];
ll binom(int n,int m){return fc[n]*ifc[m]%mod*ifc[n-m]%mod;}
ll n,k,ans;
int main(){
cin>>n>>k;
for(int i=0;i<=n;i++)fc[i]=i?fc[i-1]*i%mod:1,ifc[i]=inv(fc[i]);
for(int i=0;i<=n;i++)
for(int j=0;j<=n;j++){
ll c=(n-i)*(n-j);
ll b=binom(n,i)*binom(n,j)%mod;
ans=(ans+(n-i+n-j&1?mod-1:1)*b%mod*ksm(k,c)%mod*ksm(k-1,n*n-c))%mod;
}
cout<<ans<<endl;
return 0;
}
3. 單位根反演
3.1. 公式
常規證法是用等比數列求和公式，我不喜歡，因為這樣不能使初學單位根反演的同學感受到它的優美與自然，從而讓他們感覺單位根反演是一拍腦袋想出來的，就很 artificial。
一個比較直觀的，我自己的理解就是 \(\omega_n^x\) 在 \(n\mid x\) 時等於 \(1\)，\(n\) 個求和后再除以 \(n\) 就是 \(1\)；而在 \(n\nmid x\) 時不妨設 \(d=\gcd(x,n)\)。把單位根放到平面直角坐標系里感受一下，\(\sum_{i=0}^{n-1}(\omega_n^x)^i\) 就是把 \(\omega_n^x\) 旋轉 \(n\) 次，每次旋轉相同的 \(\dfrac{2\pi x}{n}\) 角度。那么最終形成了相同的 \(d\) 組 “向量”，每組內部有 \(\dfrac{n}{d}\) 個在單位圓上均勻分布的 “向量”，根據高中知識可知每一組向量的和都是 \(0\)。
它還有另外一種形式，在做題時經常會遇到：
3.2. 應用
我們為什么要用單位根反演？鑽研這個問題，不僅可以幫助我們深入理解單位根反演，在做題時我們也能更容易看出它的蛛絲馬跡。
3.2.1. 將 mod 轉化為求和
首先看一道經典的單位根反演入門題：給出 \(n,s,a_0,a_1,a_2,a_3\leq 10^{18}\)，求 \(\sum_{i=0}^n\dbinom{n}{i}s^ia_{i\bmod 4}\) 對 \(998244353\) 取模。
看到組合數和冪結合，自然想到二項式定理，可是 \(a_{i\bmod 4}\) 的部分斷絕了我們的出路 …… 嗎？Nope.
如果我們將 \(a_{i\bmod 4}\) 寫作 \(\sum_{j=0}^3a_j[i\equiv j\pmod 4]\)，再使用單位根反演得到 \(\sum_{j=0}^3a_j\times \dfrac{1}{4}\sum_{k=0}^3\omega_4^{(i-j)k}\)，揉進原來的柿子：
交換求和符號並整理，得
逆用二項式定理，得
如果你學過 NTT 就知道在模 \(998244353\) 意義下，原根和單位根可以互換，因此這里的四次單位根就等於 \(g^{\frac{p-1}4}\)，一個快速冪解決。時間復雜度 \(\mathcal{O}(\log n)\)，有 \(16\) 倍常數。
代碼見例題 I.
3.3. 例題
*I. LOJ#6485. LJJ 學二項式定理
ll ksm(ll a,ll b){
ll s=1;
while(b){
if(b&1)s=s*a%mod;
a=a*a%mod,b>>=1;
}
return s;
}
ll T,n,s,a[4],w[4];
int main(){
cin>>T,w[0]=1,w[1]=ksm(3,mod-1>>2),w[2]=w[1]*w[1]%mod,w[3]=w[1]*w[2]%mod;
while(T--){
cin>>n>>s;
for(int i=0;i<4;i++)cin>>a[i];
ll ans=0;
for(int i=0;i<4;i++)
for(int j=0;j<4;j++){
ll base=a[i]*w[(16-i*j)%4]%mod;
ans=(ans+base*ksm((s*w[j])%mod+1,n))%mod;
}
cout<<ans*ksm(4,mod-2)%mod<<endl;
}
return 0;
}
*II. P5591 小豬佩奇學數學
推推柿子：根據 \(\left\lfloor\dfrac{i}k\right\rfloor=\dfrac{i-i\bmod k}k\)，得
略去 \(\dfrac 1 k\)，括號拆開分別考慮。
-
前半部分：\(\sum_{i=0}^n\dbinom n i p^i\times i\)。
根據 \(\dbinom{n}{m}=\dfrac n m\dbinom {n-1}{m-1}\)，得 \(np\sum_{i=0}^n\dbinom {n-1}{i-1} p^ {i-1}\)。
即 \(np(p+1)^{n-1}\)。
-
后半部分：略去負號，得 \(\sum_{i=0}^{n}\dbinom{n}{i}p^i\sum_{j=0}^{k-1}j[i\equiv j\pmod k]\)。
套入單位根反演，得 \(\sum_{i=0}^n\dbinom n ip^i\sum_{j=0}^{k-1}\dfrac j k\sum_{x=0}^{k-1}\omega_k^{(i-j)x}\)。
略去 \(\dfrac 1 k\)，交換求和符號，得 \(\sum_{j=0}^{k-1}j\sum_{x=0}^{k-1}\omega_k^{-jx}\sum_{i=0}^n\dbinom ni(p\omega_k^x)^i\)。
套入二項式定理，得 \(\sum_{j=0}^{k-1}j\sum_{x=0}^{k-1}\omega_k^{-jx}(p\omega_k^x+1)^n\)。
然后發現不會了。交換求和順序，得 \(\sum_{x=0}^{k-1}(p\omega_k^x+1)^n\sum_{j=0}^{k-1}j(\omega_k^x)^{-j}\)。
略去前面的柿子，將 \(\omega_k^{-x}\) 記做 \(c\)，僅關注 \(\sum_{j=0}^{k-1}jc^{j}\)。由於 \(k\) 是常數，因此將其記做 \(f(c)\)。
看到這里，你發現了什么？對！就是我們小學二年級就學過的錯位相減法！
-
如果 \(c\neq 1\)，求 \(cf(c)-f(c)\)：比對每一項的系數，有 \(f(c)-cf(c)=\sum_{j=1}^{k-1}c^j-(k-1)c^k\)。
先特判掉 \(k=1\) 的情況，然后用等比數列求和公式得到 \(\sum_{j=1}^{k-1}c^j=-1+\sum_{j=0}^{k-1}c^j=-1+\dfrac{c^k-1}{c-1}\)，注意到 \(c^k\) 恆為 \(1\)，所以 \(f(c)=\dfrac{-1-(k-1)}{1-c}=\dfrac{k}{c-1}\)。
-
如果 \(c=1\)，那么原式等於 \(\sum_{j=0}^{k-1}j=\dfrac{k(k-1)}2\)。
因此直接計算即可！時間復雜度 \(\mathcal{O}(k\log n)\)。
-
ll ksm(ll a,ll b){
ll s=1; a%=mod;
while(b){
if(b&1)s=s*a%mod;
a=a*a%mod,b>>=1;
}
return s;
}
ll inv(ll x){return ksm(x,mod-2);}
ll n,p,k,w,ans,sub;
int main(){
cin>>n>>p>>k,w=ksm(3,(mod-1)/k);
ans=n*p%mod*ksm(p+1,n-1)%mod;
for(int x=0;x<k;x++){
ll om=ksm(w,x),pw=ksm(p*om+1,n),c=inv(om);
if(c==1)sub=(sub+pw*(k*(k-1)/2%mod))%mod;
else sub=(sub+pw*k%mod*inv(mod+c-1))%mod;
}
cout<<(ans+mod-sub*inv(k)%mod)*inv(k)%mod;
return 0;
}
4. 莫比烏斯反演
極其重要的一類反演，對於推式子的幫助很大。
4.1. 前置知識
4.1.1. 數論函數與積性函數
數論函數，就是定義域是正整數的函數。
積性函數，就是對於任意 \(x,y\in \N^{+}\)，若 \(\gcd(x,y)=1\)，則 \(f(xy)=f(x)f(y)\) 的函數 \(f\)。若完全積性則不需要 \(\gcd(x,y)=1\) 的條件。
一些積性函數舉例（接下來的文章中會出現）：
- 單位函數 \(\epsilon(n)=[n=1]\)。
- 常數函數 \(1(n)=1\)。
- 恆等函數 \(\mathrm{id}_k(n)=n^k\)，\(\mathrm{id}_1(n)\) 通常記作 \(\mathrm {id}(n)\)。
- 除數函數 \(\sigma_k(n)=\sum_{d|n}d^k\)。\(\sigma_0(n)\) 就是約數個數，記作 \(\tau(n)\) 或 \(\mathrm{d}(n)\)。\(\sigma_1(n)\) 就是約數和，記作 \(\sigma(n)\)。
- 大名鼎鼎的歐拉函數 \(\varphi(n)=\sum_{i=1}^n[\gcd(i,n)=1]\)。關於歐拉函數可以看 基礎數論學習筆記。
- 本章節的核心莫比烏斯函數 \(\mu(n)=\begin{cases}0&\exist\ d>1,d^2\mid n \\(-1)^{\omega(n)}&\mathrm{otherwise}\end{cases}\)，其中 \(\omega(n)\) 表示 \(n\) 本質不同質因子個數。它的積性是顯然的。
4.1.2. 狄利克雷卷積（Dirichlet）
由於學習莫比烏斯反演之前需要學會狄利克雷卷積，所以就放在這了。大名鼎鼎的狄利克雷卷積竟淪落到給莫比烏斯反演作鋪墊。
定義：設 \(f,g\) 是數論函數，定義它們的狄利克雷卷積為 \(h(x)=\sum_{ab=x}f(a)g(b)\)。顯然的，狄利克雷卷積具有交換律，結合律和乘法分配律。它的另一種表達形式為 \(h(x)=\sum_{d|x}f(d)g\left(\dfrac x d\right)\)，非常有用。
FMT 和 FWT 叫位運算卷積，FFT 和 NTT 叫加法卷積，不難發現狄利克雷卷積實際上是乘法卷積。
狄利克雷卷積有如下性質：兩個積性函數的狄利克雷卷積也是積性函數，積性函數的逆元也是積性函數。證明略。
4.1.3. 莫比烏斯函數
因為莫比烏斯函數是積性函數，所以它可以由線性篩求得。
莫比烏斯函數有以下幾條性質：
首先 \(n=1\) 時顯然，否則我們令 \(n\) 的本質不同質因數分別為 \(p_1,p_2,\cdots,p_k\)，從中選取 \(i\) 個質數的方案數為 \(\dbinom ki\)。由選奇數個方案數之和等於選偶數個方案數之和可知和為 \(0\)。即 \(k>0\) 時 \(\sum_{i=0}^k1^{k-i}\times\dbinom k i\times (-1)^{i}=(1+(-1))^k=0\)（二項式定定理）。
因此我們證明了 \(\mu*1=\epsilon\)。
該式子在推式子時比較有用。可以先證明 \(\varphi*1=\mathrm{id}\)，兩邊同時卷積 \(\mu\) 有 \(\varphi*(1*\mu)=\varphi*\epsilon=\varphi=\mathrm{id}*\mu\)。得證。
這是極其重要的一個柿子，也是莫比烏斯反演的核心。將 \(\gcd(i,j)\) 帶入式 \(4.1\) 即得。
線性篩莫比烏斯函數。
mu[1]=1;
for(int i=2;i<N;i++){
if(!vis[i])pr[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&i*pr[j]<N;j++){
vis[i*pr[j]]=1;
if(i%pr[j]==0)break;
mu[i*pr[j]]=-mu[i];
}
}
4.1.4. 整除分塊
哪里有莫反，哪里就有整除分塊。
一個經典的例子是求 \(\sum_{i=1}^n\left\lfloor\dfrac n i\right\rfloor\)，\(n\leq 10^{12}\)。注意到不同的 \(\left\lfloor\dfrac n i\right\rfloor\) 的個數實際上是 \(\sqrt n\) 級別的，故使用整除分塊：對於一個 \(i\)，我們要找到最大的 \(j\) 使得 \(\left\lfloor\dfrac n i\right\rfloor=\left\lfloor\dfrac n j\right\rfloor\)，可以證明 \(j=\left\lfloor\dfrac n {\left\lfloor\dfrac n i\right\rfloor}\right\rfloor\)。其實很好理解：找到使 \(\left\lfloor\dfrac n j\right\rfloor=c\) 的最大的 \(j\)，只需要用 \(n\) 除以 \(c\) 下取整。
那么，\(l\) 從 \(1\) 開始，每次令 \(r\gets \min\left(n,\left\lfloor\dfrac n {\left\lfloor\dfrac n l\right\rfloor}\right\rfloor\right)\)，然后將答案加上 \((r-l+1)\times \left\lfloor\dfrac n l\right\rfloor\) 即可。最后 \(l\gets r+1\)，若大於 \(n\) 則退出。
注意當被除數 \(k\) 和枚舉上界 \(n\) 不同且 \(k<n\) 時，需要特判 \(\left\lfloor\dfrac k l \right\rfloor=0\)。
4.2. 公式
若 \(f,g\) 是兩個數論函數，那么有
證明：已知 \(g*1=f\)，則 \(g*1*\mu=g=f*\mu\)。
證明（來自 OI - Wiki）：
4.3. 直觀理解
只會套皮是不行的，要想熟練運用，不僅要多刷題，也要理解莫反的本質：容斥原理。其中莫比烏斯函數就是容斥系數。也許你會問：容斥系數不應該是 \((-1)^k\) 嗎？其實在某種程度上沒錯，畢竟對比 \((-1)^k\) 與莫比烏斯函數的定義，它們的形式非常相近：只有當 \(x\) 是完全平方數的倍數時容斥系數才為 \(0\)（這其實也暗示了這樣的 \(x\) 不同於別的數）。但我們是對於全體自然數做容斥，會有一些奇怪的事情發生：一個自然數被另一個自然數完全包含（這也說明它一定是完全平方數的倍數，接下來會解釋），而在容斥原理中，一個集合被另一個集合包含這樣的事情是不會發生的。
不妨設接下來有 \(f_n=\sum_{d\mid n}g_d\)。
4.3.1. \(n=2\)
當 \(n=2\) 時，為了求出 \(g_2\)，我們只需要用 \(f_2-f_1=(g_2+g_1)-g_1\)。這里已經可以看見容斥的影子了。不妨結合下圖以更好理解。
4.3.2. \(n=6\) 及 \(\prod p_i\)
當 \(n=6\) 時，為了求出 \(g_6\)，我們首先需要 \(f_6\)，但是這樣會多加 \(g_1+g_2+g_3\)，所以減去 \(f_2+f_3\)，但是這樣又多減了 \(g_1\)，所以我們還應加上 \(f_1\)。結合下圖以更好理解。
如果 \(n\) 是若干個不同質數 \(p_1,p_2,\cdots,p_k\) 的乘積，那么畫出能夠表示 \(n\) 的所有約數的韋恩圖（\(k>3\) 時畫不出來，不過可以理解一下）。\(k\) 個集合相交形成的所有子集就是 \(n\) 的所有約數。比如說，如果有一塊區域落在集合 \(p_1\) 和集合 \(p_3\) 中，那么它所表示的約數就是 \(p_1p_3\)。顯然，如果 \(x\mid y\)，那么區域 \(y\) 包含於區域 \(x\)。
圖源 https://blog.csdn.net/Summer__show_/article/details/76269088 。
不妨將 \(f_d\) 看做整塊區域 \(d\)（比如說上圖中區域 \(30\) 也屬於區域 \(6\)），而 \(g_d\) 看做僅屬於區域 \(d\) 的部分，然后將所有區域表示的數 \(i\) 替換成 \(\dfrac n i\)，那么 \(f_n\) 就表示最外層的大區域，即所有 \(g_d\) 的和，而 \(f_1\) 就表示最里層的小區域，即 \(g_1\)，剛好符合 \(f_n=\sum_{d\mid n} g_d\) 的條件。因此，由於 \(f\) 已知，為了得到僅屬於區域 \(n\) 的部分所表示的數 \(g_n\)，我們只需要進行容斥：加上 \(f_n\)，減去 \(f_{p_1},\cdots, f_{p_k}\) …… 不難發現每個數前的容斥系數就是 \(-1\) 的 \(\dfrac n i\) 所含質因子個數次冪，即 \(\mu\left(\dfrac n i\right)\)（\(n\) 由互不相同的質因子相乘得到，那么 \(\dfrac n i\) 也如此）。這不就是 \(g_n=\sum_{d\mid n}f_d\mu\left(\dfrac n d\right)\) 嗎！
4.3.3. \(n=4\) 及任意自然數
有了上面的經驗，你應該很快知道表示約數 \(4\) 的區域被表示約數 \(2\) 的區域 “完全包含”，這里的 \(y\) 被 \(x\) “完全包含” 其實指的是：\(y\) 包含於的所有集合等於 \(x\) 包含於的所有集合且 \(y\mid x\)。畫出韋恩圖就是一大一小兩個圓，其中一個圓完全落在另一個里面。因此，如果計算了 \(2\) 的貢獻，那么完全沒有必要再算 \(4\) 了，因為 \(4\) 被 \(2\) “完全包含”，所以它們的 “多加/多減狀態” 應該是相同的。也就是說算完 \(2\) 的貢獻，\(4\) 的貢獻也恰好算完，既不多加，也不多減，因而 \(\mu(4)=0\)。不難發現，這種情況僅在一個數是非 \(1\) 的完全平方數的倍數時出現，從而推廣至任意自然數。
綜上，莫比烏斯函數可以看作對 \(\N^{+}\) 做容斥時的容斥系數。
4.4. 應用
一般用來化簡帶有 \(\gcd\) 的柿子，幾乎所有莫反題都需要用到整除分塊。
具體應用可以看 OI-Wiki，抄在博客里其實意義不大。
4.5. 例題
下面所有分式都是下取整。
I. P2522 [HAOI2011]Problem b
二維差分將下界化為 \(1\)，然后推式子
將 \(k\) 搞掉：
莫反：
枚舉約數 \(d\)，記 \(c=\min\left(\dfrac n k,\dfrac m k\right)\)：
由於 \(1\sim x\) 中 \(y\) 的倍數有 \(\dfrac x y\) 個，故原式化為：
整除分塊即可，時間復雜度 \(\mathcal{O}(n+T\sqrt n)\)。
long long 改成 int，11s 變成 5s，離譜。
const int N=5e4+5;
int cnt,pr[N],vis[N],mu[N];
int a,b,c,d,k;
int calc(int x,int y){
x/=k,y/=k;
int c=min(x,y),ans=0;
for(int l=1,r;l<=c;l=r+1){
r=min(x/(x/l),y/(y/l));
ans+=(x/l)*(y/l)*(mu[r]-mu[l-1]);
} return ans;
}
void solve(){
a=read(),b=read(),c=read(),d=read(),k=read();
cout<<calc(b,d)-calc(b,c-1)-calc(a-1,d)+calc(a-1,c-1)<<"\n";
}
int main(){
mu[1]=1;
for(int i=2;i<N;i++){
if(!vis[i])pr[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&i*pr[j]<N;j++){
vis[i*pr[j]]=1;
if(i%pr[j]==0)break;
mu[i*pr[j]]=-mu[i];
} mu[i]+=mu[i-1];
}
int T=read();
while(T--)solve();
return 0;
}
II. SP5971 LCMSUM - LCM Sum
來推式子。
*III. P4318 完全平方數
莫比烏斯函數作為容斥系數的經典例子。
顯然答案滿足可二分性，於是考慮如何求 \(x\) 以內不含有平方因子的數的個數。仍然是容斥：我們只需要加上所有由偶數個質數相乘的平方的倍數，減去由奇數個質數相乘的平方的倍數。發現這就是 \(\sum_{i^2\leq n}\mu(i)\dfrac{n}{i^2}\)，直接暴力求就行了。時間復雜度 \(\mathcal{O}(T\log k\sqrt k)\)。
const int N=1e5+5;
int T,n,cnt,pr[N],vis[N],mu[N];
ll calc(ll x){
ll res=0;
for(int i=1;i*i<=x;i++)res+=mu[i]*(x/i/i);
return res;
}
int main(){
cin>>T,mu[1]=1;
for(int i=2;i<N;i++){
if(!vis[i])pr[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&i*pr[j]<N;j++){
vis[i*pr[j]]=1;
if(i%pr[j]==0){
mu[i*pr[j]]=0;
break;
}
mu[i*pr[j]]=-mu[i];
}
}
while(T--){
cin>>n; ll l=n,r=n<<1;
while(l<r){
ll m=l+r>>1;
if(calc(m)<n)l=m+1;
else r=m;
} cout<<l<<endl;
}
return 0;
}
IV. [BZOJ2986]Non-Squarefree Numbers
雙倍經驗 P4318 完全平方數。代碼從 \(\omega(i)=1\) 時開始容斥，此時容斥系數為 \(1\)（從 \(\omega (i)=0\) 即 \(i=1\) 開始容斥也可以），故篩出的莫比烏斯函數是實際定義的相反數。
template <class T> void cmin(T &a, T b){a = a < b ? a : b;}
template <class T> void cmax(T &a, T b){a = a > b ? a : b;}
bool Mbe;
const int N = 2e5 + 5;
int mu[N], pr[N], vis[N], cnt;
ll n, l, r;
void sieve() {
for(int i = 2; i < N; i++) {
if(!vis[i]) pr[++cnt] = i, mu[i] = 1;
for(int j = 1; j <= cnt && i * pr[j] < N; j++) {
vis[i * pr[j]] = 1;
if(i % pr[j] == 0) break;
mu[i * pr[j]] = mu[i] * -1;
}
}
}
ll check(ll lim) {
ll res = 0;
for(ll i = 2; i * i <= lim; i++) res += lim / i / i * mu[i];
return res;
}
bool Med;
int main(){
fprintf(stderr, "%.2lf\n", (&Mbe - &Med) / 1048576.0);
sieve(), cin >> n, l = 1, r = n << 2;
while(l < r) {
ll m = l + r >> 1;
if(check(m) < n) l = m + 1;
else r = m;
}
cout << l << endl;
return 0;
}