莫比烏斯反演入門
莫比烏斯反演
莫比烏斯反演是數論中的一個重要內容，可以化簡很多數論函式的計算。
本文形式化過程偏多，一定要耐心看完並試著自己推導。
前置芝士
>莫比烏斯函式<
定義
對於定義在自然數域的兩個函式 \(F(x)\) 與 \(f(x)\) ，若兩函式滿足
\[F(n)=\sum_{d|n}f(d) \]則有
\[f(n)=\sum_{d|n}\mu(d)F(\frac{n}{d}) \]莫比烏斯反演還有另外一種常用的形式:
\[F(n)=\sum_{n|d}f(d)\\ \Rightarrow f(n)=\sum_{n|d}\mu(\frac{d}{n})F(d) \]證明：
\[\begin{aligned} \sum_{d|n}\mu(d)F(\frac{n}{d}) & =\sum_{d|n}\mu(d)\sum_{d^{\prime}|\frac{n}{d}}f(d^{\prime})\\ & =\sum_{d^{\prime}|n}\mu(d)\sum_{d|\frac{n}{d}}f(d)\\ & =f(n) \end{aligned} \]使用例：
求
\[\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=1] \]不妨假設這裡的\(n<m\)。
對於gcd這個東西很套路，先設
\[f(x)=\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=x]\\ g(x)=\sum_{x|d}f(d) \]其中 \(d\leq n\)。
代入莫比烏斯反演公式2，得到：
\[f(x)=\sum_{x|d}\mu(\frac{d}{x})g(d)\\ \Rightarrow f(1)=\sum_{1|d}\mu(d)g(d)\\ \]修改列舉項，可以得到
\[f(1)=\sum_{i=1}^{n}\mu(i)g(i) \]\(\mu(i)\) 再好求不過了，關鍵在於 \(g(i)\) 是個什麼東西。
從含義出發，易得到：
\[g(x)=\sum_{i=1}^{n}\sum_{j=1}^{m}[x|gcd(i,j)]\\ g(x)=\sum_{i=1}^{n/x}\sum_{j=1}^{m/x}[1|\gcd(i,j)]\\ \]\(1|gcd(i,j)\) 顯然恆成立，則
\[g(x)=\frac{n}{x}\cdot \frac{m}{x} \]\(O(1)\) 算就可以了。
最後我們求的是
\[Ans=f(1)=\sum_{i=1}^{n}\mu(i)\lfloor\frac{n}{i}\rfloor\lfloor\frac{m}{i}\rfloor \]所以只需要 \(O(n)\) 算 \(\mu(i)\) 字首和即可。
最後程式碼和>莫比烏斯函式<裡的沒有太大區別
例題
P2398 GCD SUM
>題目連結<
題目描述
求
\[\sum_{i=1}^{n}\sum_{j=1}^{n}gcd(i,j) \]輸入格式
一行一個正整數 \(n\)
輸出格式
一行一個整數表示答案
資料範圍
\(n\leq 10^5\)
解析
這裡直接講 \(\sum_{i=1}^{n}\sum_{j=1}^{m}gcd(i,j)\) 的做法
直接拿給我們是沒有辦法做的。
不妨設\(n < m\)，實際程式碼中 \(n=min(n,m)\) 即可。這類有限項列舉求和交換列舉項順序一般是沒有問題的。
考慮列舉gcd，能化得如下式子
\[\sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m}d\cdot[gcd(i,j)=d]\\ =\sum_{d=1}^{n}d\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=d]\\ =\sum_{d=1}^{n}d\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)=1] \]後面這一團不就是之前推的式子？
直接套結論
\[\sum_{d=1}^{n}d\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)=1]\\ =\sum_{d=1}^{n}d\sum_{i=1}^{n/d}\mu(i)\lfloor\frac{n}{id}\rfloor\lfloor\frac{m}{id}\rfloor \]可以知道，對於特定區間內的 \(d\) ，\(n/d\) 是無變化的。
所以可以對前面的 \(d\) 進行數論分塊。
後面的求值參考上面，也可以數論分塊。
總時間複雜度 \(O(\sqrt{n}\times\sqrt{n})=O(n)\) 。
落到本題上，本題的 \(n=m\) , 數論分塊的取 \(min(n/(n/i),m/(m/i))\) 操作都免了。
參考code:
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e6;
int primes[N],tot=0;
bool mp[N];
int mu[N],sum_[N];
void init(int n)
{
mu[1]=1;
for(int i=2; i<=n; i++)
{
if(!mp[i])
{
primes[++tot]=i;
mu[i]=-1;
}
for(int j=1; primes[j]*i<=n; j++)
{
int tmp=primes[j]*i;
mp[tmp]=1;
if(i%primes[j]==0)
{
mu[tmp]=0;
break;
}
mu[tmp]=mu[i]*-1;
}
}
for(int i=1; i<=n; i++)
{
sum_[i]=sum_[i-1]+mu[i];
}
}
ll solve(ll a,ll b)
{
ll res=0;
for(int i=1,j; i<=a; i=j+1)
{
j=min(a/(a/i),b/(b/i));
res+=(ll)(sum_[j]-sum_[i-1])*(a/i)*(b/i);
}
return res;
}
int main()
{
init(1e5+10);
int n;
scanf("%d",&n);
ll ans=0;
for(int i=1,j; i<=n; i=j+1) //第一遍整除分塊前面的d
{
j=n/(n/i);
ll tmp=(ll)(i+j)*(j-i+1)/2;
ans+=(ll)tmp*solve(n/i,n/i);//第二遍整除分塊
}
printf("%lld",ans);
return 0;
}
P2568 GCD
>題目連結<
題目描述
給定正整數 \(n\)，求 \(1\le x,y\le n\) 且 \(\gcd(x,y)\) 為素數的數對 \((x,y)\) 有多少對。
輸入格式
只有一行一個整數，代表 \(n\)。
輸出格式
一行一個整數表示答案。
資料範圍
\(1\le n \le 10^7\)
解析
題目求
\[\sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)\text{是素數}] \]令 \(gcd(i,j)=d\)，則
\[\text{原式}=\sum_{d=1}^{n}d[d\text{是素數}]\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}[gcd(i,j)=1] \]後面那個東西怎麼做不用重複說了吧。
至於前面，對 \(d\) 再進行一次數論分塊 (有向下取整除法的地方就可以考慮一下數論分塊)
再處理一個素數個數字首和就行了。
code:(僅供思路參考，下面這份程式碼以 2.82s/247.67MB 的好成績驚險卡過)
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e7+10;
ll primes[N],psum[N],tot=0;
bool mp[N];
ll mu[N],sum[N];
void init(ll n)
{
mu[1]=1,mp[1]=1;
for(int i=2; i<=n; i++)
{
if(!mp[i])
{
primes[++tot]=i;
mu[i]=-1;
}
for(int j=1; (ll)primes[j]*i<=n; j++)
{
ll tmp=(ll)primes[j]*i;
mp[tmp]=1;
if(i%primes[j]==0)
{
mu[tmp]=0;
break;
}
mu[tmp]=mu[i]*-1;
}
}
for(int i=1; i<=n; i++)
{
sum[i]=sum[i-1]+mu[i];
if(mp[i]) psum[i]=psum[i-1];
else psum[i]=psum[i-1]+1;
}
}
ll solve(ll a,ll b)
{
ll res=0;
for(int i=1,j; i<=a; i=j+1)
{
j=min(a/(a/i),b/(b/i));
res+=(ll)(sum[j]-sum[i-1])*(a/i)*(b/i);
}
return res;
}
int main()
{
init(1e7+10);
int n;
scanf("%d",&n);
ll ans=0;
for(int i=1,j; i<=n; i=j+1)
{
j=n/(n/i);
ans+=(ll)(psum[j]-psum[i-1])*solve(n/i,n/i);
}
printf("%lld",ans);
return 0;
}
P2257 YY的GCD
>題目連結<
題目描述
神犇 YY 虐完數論後給傻× kAc 出了一題
給定 \(N, M\)，求 \(1 \leq x \leq N\) ，\(1 \leq y \leq M\) 且 \(\gcd(x, y)\) 為質數的 \((x, y)\) 有多少對。
輸入格式
第一行一個整數 \(T\) 表述資料組數。
接下來 \(T\) 行，每行兩個正整數，\(N, M\)。
輸出格式
共 \(T\) 行，每行一個整數表示第 \(i\) 組資料的結果。
資料範圍
\(T=10^4\ , \ N,M\le 10^7\)
解析
乍看跟上面那個題幾乎完全相同
但是是多組詢問，按照上面那個級別的優秀時空複雜度是肯定過不了的。
再往下推一下式子(這裡已經把 \(n,m\) 換上去了，\(n,m\) 相當於題中的 \(N,M\) ，且 \(n\le m\) )
令 \(T=id\)
\[\begin{aligned} Ans & =\sum_{d=1}^{n}d[d\text{是素數}]\sum_{i=1}^{n/d}\mu(i)\lfloor\frac{n}{id}\rfloor\lfloor\frac{m}{id}\rfloor\\ & =\sum_{T=1}^{n}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{d|T}[d\text{是素數}]\mu(\frac{T}{d}) \end{aligned} \]按照慣例，我們要想一下後面那團東西的字首和怎麼求。
可以從線性篩出發考慮。
令函式 \(\lambda(x)=\sum_{d|x}[d\text{是素數}]\mu(\frac{x}{d})\) ，求這個函式的字首和。
首先，\(\lambda(1)=0,\lambda(\text{質數})=1\) 。
設 \(x=i\cdot y\) ， \(y\) 是 \(x\) 的最小質因子。
1. 當 \(y|i\) 時，即 \(x\) 有多個最小質因子：
-
若 \(i\) 質因子互不相同時，當且僅當列舉的素數 \(d=y\) 時 \(\mu(\frac{x}{d})\ne 0\)。
此時：\(\lambda(x)=\mu(\frac{x}{y})\)
-
若 \(i\) 有相同質因子，則對於任意列舉的素數 \(d\) ，\(\frac{x}{d}\) 內都有相同質因子， 即 \(\mu(\frac{x}{d})=0\)。
此時仍有 \(\lambda(x)=\mu(\frac{x}{y})\)
2. 當 \(y\nmid i\) 時，即 \(x\) 只有一個最小質因子：
\[\lambda(x)=\sum_{d|x}[d\text{是素數}]\mu(\frac{x}{d})\\ \Rightarrow \lambda(i\cdot y)=\sum_{d|(i\cdot y)}[d\text{是素數}]\mu(\frac{i\cdot y}{d}) \]我們已經知道，\(\mu(i\cdot y)=-\mu(i)\) (如果看不懂可以再看看莫比烏斯函式定義)
所以對於 \(\lambda(i)\) 其中的每一項的相反數都被包括在了 \(\lambda(x)\) 中，且 \(\lambda(x)\) 只是多了一個 \(\mu(\frac{i\cdot y}{y})=\mu(i)\)。
所以此時的 \(\lambda(x)=-\lambda(i)+\mu(i)\)
綜上：
\[\lambda(x)= \begin{cases} \mu(1) , & x\text{是質數}\\ \mu(i) , & y|i\\ -\lambda(i)+\mu(i) , & y\nmid i \end{cases} \]於是 \(\lambda(x)\) 就可以用線篩求了。
code: 7.79s/90.86MB
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e7+10;
int primes[N],tot=0;
int mu[N],lam[N];
bool mp[N];
void init(int n)
{
mu[1]=mp[1]=1;
for(int i=2; i<=n; i++)
{
if(!mp[i]) primes[++tot]=i,mu[i]=-1,lam[i]=1;
for(int j=1; primes[j]*i<=n; j++)
{
int x=i*primes[j];
mp[x]=1;
if(i%primes[j]==0)
{
lam[x]=mu[i];
mu[i]=0;
break;
}
lam[x]=-lam[i]+mu[i];
mu[x]=-1*mu[i];
}
}
for(int i=1; i<=n; i++)
lam[i]+=lam[i-1];//字首和
}
int main()
{
init(1e7);
int t;
scanf("%d",&t);
while(t--)
{
int n,m;
scanf("%d%d",&n,&m);
ll ans=0;
if(n>m) swap(n,m);
for(int i=1,j; i<=n; i=j+1)
{
j=min(n/(n/i),m/(m/i));
ans+=(ll)(lam[j]-lam[i-1])*(n/i)*(m/i);
}
printf("%lld\n",ans);
}
return 0;
}
P1829 [國家集訓隊]Crash的數字表格/JZPTAB
P5518 [MtOI2019]幽靈樂團 / 莫比烏斯反演基礎練習題
(在做了在做了)