前言
莫比乌斯反演是数论中的重要内容
对于一些函数(f(n))，如果很难直接求出它的值
而容易求出其倍数和或约数和(g(n))，那么可以通过莫比乌斯反演简化运算，求得(f(n))的值
前置知识：整除分块
对于类似于 (sum_{i=1}^{n}lfloorfrac{n}{i} floor) 的式子
我们可以很容易地用 (O(n)) 的效率求出它的值
但是在一些题目中，我们需要更优秀的时间效率，这时候就要用到整除分块
主要的思想是对于任意一个 (i(i leq n)) 找到一个最大的 (j)
使得 (i leq j leq n) 并且 (lfloorfrac{n}{i}
floor=lfloorfrac{n}{j}
floor)
此时 (j=leftlfloorfrac{n}{leftlfloorfrac{n}{i} ight floor} ight floor)
因为对于 $lfloorfrac{n}{i} floor $，当 (i leq sqrt{n}) 时只有 (sqrt{n}) 种，当 (i>sqrt{n})时，(lfloorfrac{n}{i} floor< sqrt{n})，也只有 (sqrt{n}) 种取值，共计 (2sqrt{n}) 种
所以该过程的时间复杂度为 (O(sqrt{n}))
关于正确性的证明
首先显然 (j leq n)
又有 (j=leftlfloorfrac{n}{leftlfloorfrac{n}{i} ight floor} ight floor geq lfloor frac{n}{frac{n}{i}} floor =i)
设 (k=lfloor frac{n}{i} floor),那么我们要证明的就是 当$ lfloor frac{n}{j} floor=k(时，)j$ 的最大值为 (lfloor frac{n}{k} floor)
因为 (j) 为整数，所以 (max(j)=lfloor frac{n}{k} floor)
练习题
分析
(ans=sum_{i=1}^n(kmod i)=sum_{i=1}^nk-ileftlfloorfrac{k}{i} ight floor)
代码
#include<cstdio>
#include<cmath>
#include<algorithm>
#define rg register
int n,k,mmax;
long long ans=0;
int main(){
scanf("%d%d",&n,&k);
mmax=std::min(n,k);
for(int l=1,r;l<=mmax;l=r+1){
r=k/(k/l);
r=std::min(r,mmax);
ans+=1LL*(k/l)*(l+r)*(r-l+1)/2;
}
ans=1LL*n*k-ans;
printf("%lld
",ans);
return 0;
}
拓展：二维整除分块
求 $ sum_{i=1}^{min (n,m)}leftlfloorfrac{n}{i} ight floorleftlfloorfrac{m}{i} ight floor $的值
此时我们需要将代码中 (r) 的取值改为 (min(m/(m/l),n/(n/l)))
莫比乌斯函数
(μ(d)) 为莫比乌斯函数的函数名，是一个由容斥系数所构成的函数。
其定义为：
(1)、当 (d=1) 时， (mu(d)=1) 。
(2)、当 (d = prod_{i=1}^k p_i)，并且所有的 (p_i) 为不同的素数的时候， (mu(d)=(-1)^k) 。
(3)、当 (d) 含有的任一质因子的幂大于等于 (2) 次，那么 (mu(d)=0) 。
求法（线性筛）
const int maxn=4e5+5;
bool not_pri[maxn];
int t,mu[maxn],pri[maxn];
void xxs(){
not_pri[0]=not_pri[1]=1;
mu[1]=1;
for(rg int i=2;i<maxn;i++){
if(!not_pri[i]){
pri[++pri[0]]=i;
mu[i]=-1;
}
for(rg int j=1;j<=pri[0] && 1LL*i*pri[j]<maxn;j++){
not_pri[i*pri[j]]=1;
if(i%pri[j]==0) break;
mu[i*pri[j]]=-mu[i];
}
}
}
性质一
对于任意正整数 (n)，(sum_{d|n}mu(d)=[n=1])
其中 ([n=1]) 表示只有当 (n=1) 成立的时候返回值为 (1)，否则返回值为 (1)
在反演的时候，这一个式子经常会将 (gcd(i,j)=1) 这个条件代换掉
证明
当 (n=1)时，返回值显然为 (1)
当 (n eq1)时，根据唯一分解定理，(n) 一定能划分为若干质因数的乘积
我们设这样的质因数一共有 (k) 种
如果任意一种质因数选择了大于一个，那么组成的数的 (mu) 值一定为 (0)
所以每种质因子最多选择 (1) 个
所以选择 (i) 种质因子的方案数为 (C_k^i)，返回值为 ((-1)^i)
贡献的总和为 (sum_{i=1}^k C_k^i(-1)^i)
如果再加上一项 (C_k^0(-1)^0)
根据二项式定理，恰好是 ((1-1)^k=1)
减去加上的一项，结果就是 (0)
性质二
对于任意正整数(n)，(sum_{d|n}frac{mu(d)}{d}=frac{varphi(n)}{n})
这个性质阐述了莫比乌斯函数和欧拉函数之间的联系
证明
(varphi(n)=sum_{i=1}^n[gcd(n,i)=1]=sum_{i=1}^{n} sum_{d|gcd(n,i)} mu(d)=sum_{d|n} mu(d) imes frac{n}{d})
转换一下就变成了(sum_{d|n}frac{mu(d)}{d}=frac{varphi(n)}{n})
莫比乌斯反演
形式一
(F(n)=sum_{d|n}f(d)=>f(n)=sum_{d|n}mu(d)F(frac{n}{d}))
证明
(f(n)=sum_{d|n}mu(d)F(frac{n}{d})=sum_{d|n}mu(d)sum_{k|frac{n}{d}}f(k)=sum_{k|n}f(k)sum_{d|frac{n}{k}}mu(d))
根据莫比乌斯函数性质一，只有当 (n=1) 时，(sum_{d|n} mu(d)=1)
所以只有当 (n=k) 时，求出的值才不为 (0)
所以右边等于 (f(n))
形式二
(F(n)=sum_{n|d}f(d)=>f(n)=sum_{n|d}mu(frac{d}{n})F(d))
证明
令(k=frac{d}{n})，
(f(n)=sum^{+infty}_{k=1}mu(k)F(nk)=sum^{+infty}_{k=1}mu(k)sum_{nk|t}f(t)=sum_{n|t}f(t)sum_{k|frac{t}{n}}mu(k))
当且仅当 (t=n) 时，(sum_{k|frac{t}{n}}mu(k)=1)，其它情况都为 (0)
所以右边等于 (f(n))
例题一、ZAP-Queries
分析
要求的式子是
(sum_{i=1}^a sum_{j=1}^b[gcd(i,j)=k])
默认 (a>b)
把 (a) 和 (b) 的上界同时除以 (k)，可以得到
(sum_{i=1}^{a/k} sum_{j=1}^{b/k}[gcd(i,j)=1])
把后面的部分用莫比乌斯函数带入，可以得到
(sum_{i=1}^{a/k} sum_{j=1}^{b/k}sum_{d|gcd(i,j)}mu(d))
把 (d) 提到前面
(sum_{d=1}^bmu(d)sum_{i=1}^{a/kd} sum_{j=1}^{b/kd})
也就是
(sum_{d=1}^bmu(d) frac{a}{kd} frac{b}{kd})
我们用整除分块处理一下，就可以做到单次询问 (sqrt{b}) 的复杂度
代码
#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=4e5+5;
bool not_pri[maxn];
int t,n,m,d,mu[maxn],pri[maxn],sum[maxn];
void xxs(){
not_pri[0]=not_pri[1]=1;
mu[1]=1;
for(rg int i=2;i<maxn;i++){
if(!not_pri[i]){
pri[++pri[0]]=i;
mu[i]=-1;
}
for(rg int j=1;j<=pri[0] && 1LL*i*pri[j]<maxn;j++){
not_pri[i*pri[j]]=1;
if(i%pri[j]==0) break;
mu[i*pri[j]]=-mu[i];
}
}
for(rg int i=1;i<maxn;i++){
sum[i]=sum[i-1]+mu[i];
}
}
int main(){
xxs();
t=read();
while(t--){
n=read(),m=read(),d=read();
if(n>m) std::swap(n,m);
n/=d,m/=d;
rg int mmax=n;
rg long long ans=0;
for(rg int l=1,r;l<=mmax;l=r+1){
r=std::min(n/(n/l),m/(m/l));
if(r>mmax) r=mmax;
ans+=1LL*(n/l)*(m/l)*(sum[r]-sum[l-1]);
}
printf("%lld
",ans);
}
return 0;
}
例题二、YY的GCD
分析
要求的式子是
(sum_{i=1}^n sum_{j=1}^m[gcd(i,j)=k](k in prime))
默认 (n>m)
和上一道题一样，我们把 (k) 提出来
(sum_{k=1}^msum_{i=1}^{n/k} sum_{j=1}^{m/k}[gcd(i,j)=1](k in prime))
把后面的部分用莫比乌斯函数带入，可以得到
(sum_{k=1}^msum_{i=1}^{n/k} sum_{j=1}^{m/k}sum_{d|gcd(i,j)}mu(d)(k in prime))
把 (d) 扔到前面
(sum_{k=1}^msum_{d=1}^m mu(d)sum_{i=1}^{n/kd}sum_{j=1}^{m/kd}(k in prime))
其实就是
(sum_{k=1}^msum_{d=1}^m mu(d)frac{n}{kd}frac{m}{kd}(k in prime))
设(kd=T)，改为枚举 (T)，则有
(sum_{T=1}^m sum_{k=1}^mmu(T/k)frac{n}{T}frac{m}{T}(k in prime))
把求和符号换一下位置
(sum_{T=1}^m frac{n}{T}frac{m}{T}sum_{k=1}^mmu(T/k)(k in prime))
前半部分可以整除分块计算，后半部分可以用 (Dirichlet) 前缀和预处理
代码
#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
#define rg register
typedef long long ll;
const int maxn=1e7+5;
bool not_pri[maxn];
int n,m,mu[maxn],pri[maxn];
ll sum[maxn];
void xxs(){
not_pri[0]=not_pri[1]=1;
mu[1]=1;
for(rg int i=2;i<maxn;i++){
if(!not_pri[i]){
pri[++pri[0]]=i;
mu[i]=-1;
}
for(rg int j=1;j<=pri[0] && 1LL*i*pri[j]<maxn;j++){
not_pri[i*pri[j]]=1;
if(i%pri[j]==0) break;
mu[i*pri[j]]=-mu[i];
}
}
for(rg int i=1;i<=pri[0];i++){
for(rg int j=1;1LL*j*pri[i]<maxn;j++){
sum[j*pri[i]]+=mu[j];
}
}
for(rg int i=1;i<maxn;i++){
sum[i]+=sum[i-1];
}
}
void solve(int n,int m){
ll ans=0;
for(rg int l=1,r;l<=n;l=r+1){
r=std::min(n/(n/l),m/(m/l));
ans+=1LL*(sum[r]-sum[l-1])*(n/l)*(m/l);
}
printf("%lld
",ans);
}
int main(){
int n,m,t;
scanf("%d",&t);
xxs();
while(t--){
scanf("%d%d",&n,&m);
if(n>m) std::swap(n,m);
solve(n,m);
}
return 0;
}
例题三、gcd
题目描述
有 (n) 个正整数 (x_1 sim x_n)，初始时状态均为未选。
有 (m) 个操作，每个操作给定一个编号 (i)，将 (x_i) 的选取状态取反。
每次操作后，你需要求出选取的数中有多少个互质的无序数对。
输入格式
第一行两个整数 (n,m)。第二行 (n) 个整数 (x_1 sim x_n)。接下来 (m) 行每行一个整数。
输出格式
(m) 行，每行一个整数表示答案。
样例
样例输入
4 5
1 2 3 4
1
2
3
4
1
样例输出
0
1
3
5
2
数据范围与提示
对于 (20\%) 的数据，(n,m<=1000)。
对于另外 (30\%) 的数据，(x_i<=100)。
对于 (100\%) 的数据，(n,m<=200000,x_i<=500000,1<=i<=n)。
分析
我们设 (s[i]) 为当前选取的数中有多少数是 (i) 的倍数
设 (f[i]) 为 (gcd) 是 (i) 的倍数的数对的个数
那么就有 (f[i]=frac{s[i](s[i]+1)}{2})
设 (g[i]) 为 (gcd) 是 (i) 的数对的个数
我们要求的就是 (g[1])
根据莫比乌斯反演的形式二
(f[i]=sum_{i|d} g[d])
则
(g[i]=sum_{i|d} mu(frac{d}{i})f[d])
而 (f[d]) 可以在 (O(sqrt{n})) 的时间内更新
所以我们只要开一个变量一直记录答案即可
代码
#include<cstdio>
#include<iostream>
#include<cmath>
#include<algorithm>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=1e6+5;
bool not_pri[maxn],vis[maxn];
int n,m,mu[maxn],pri[maxn],a[maxn],cnt[maxn];
long long g[maxn],ans;
void xxs(){
not_pri[0]=not_pri[1]=1;
mu[1]=1;
for(rg int i=2;i<maxn;i++){
if(!not_pri[i]){
pri[++pri[0]]=i;
mu[i]=-1;
}
for(rg int j=1;j<=pri[0] && 1LL*i*pri[j]<maxn;j++){
not_pri[i*pri[j]]=1;
if(i%pri[j]==0) break;
mu[i*pri[j]]=-mu[i];
}
}
}
int main(){
xxs();
n=read(),m=read();
for(rg int i=1;i<=n;i++){
a[i]=read();
}
rg int aa,bb,cc,now;
for(rg int i=1;i<=m;i++){
aa=read();
vis[aa]^=1;
cc=a[aa];
bb=sqrt(cc);
if(vis[aa]){
for(rg int j=1;j<=bb;j++){
if(cc%j==0){
now=j;
cnt[now]++;
ans-=1LL*g[now]*mu[now];
g[now]=1LL*(cnt[now]-1)*cnt[now]/2;
ans+=1LL*g[now]*mu[now];
if(j*j!=cc){
now=cc/j;
cnt[now]++;
ans-=1LL*g[now]*mu[now];
g[now]=1LL*(cnt[now]-1)*cnt[now]/2;
ans+=1LL*g[now]*mu[now];
}
}
}
} else {
for(rg int j=1;j<=bb;j++){
if(cc%j==0){
now=j;
cnt[now]--;
ans-=1LL*g[now]*mu[now];
g[now]=1LL*(cnt[now]-1)*cnt[now]/2;
ans+=1LL*g[now]*mu[now];
if(j*j!=cc){
now=cc/j;
cnt[now]--;
ans-=1LL*g[now]*mu[now];
g[now]=1LL*(cnt[now]-1)*cnt[now]/2;
ans+=1LL*g[now]*mu[now];
}
}
}
}
printf("%lld
",ans);
}
return 0;
}