【瞎口胡】数学 5 莫比乌斯反演
莫比乌斯反演用来解决一类数论问题。
在开始之前,让我们先了解一些基础知识。
数论函数
在某一整数集合上(一般为全体正整数集)定义的函数称为数论函数。
还记得积性函数的概念吗?对于函数 \(f\) 和互质的整数 \(a,b\),如果 \(f(ab)=f(a)f(b)\),那么 \(f\) 就是积性函数。如果 \(a,b\) 不互质 \(f(ab)=f(a)f(b)\) 仍然成立,那么 \(f\) 就是完全积性函数。
举几个常见的数论函数例子,它们也是积性函数:
-
常数函数 \(I\)(完全积性函数)
\(I(n)=1\)
有时也写做 \(1(n)\)。
-
单位元函数 \(\epsilon\)(完全积性函数)
\(\epsilon(n)=[n=1]\)
有时也写作 \(e(n)\)。
-
恒等函数 \(\operatorname{id}\)(完全积性函数)
\(\operatorname{id}(n)=n\)
有时也写作 \(Id(n)\)。
-
约数个数函数 \(\tau\)
\(\tau(n)=\sum \limits_{d|n} 1\)
-
约数和函数 \(\sigma\)
\(\sigma(n)=\sum\limits_{d|n} d\)
-
欧拉函数 \(\varphi\)
\(\varphi(n)=\sum \limits_{i=1}^{n} [\gcd(n,i)=1]\)
-
莫比乌斯函数 \(\mu\)
\(\mu(n) = \begin{cases} 0 & \exists i(p_i>1)\\ 1&n=1 \or \forall i(p_i=1) \and 2 \mid m \\ -1 & \text{otherwise.}\end{cases}\)
其中 \(n=\prod\limits_{i=1}^{m} p_i^{k_i}\),\(p_i\) 是两两不同的质数。
狄利克雷卷积
定义
对于数论函数 \(f,g\),定义二者的狄利克雷卷积为 \(f*g\),其中
性质
容易验证,狄利克雷卷积具有以下性质:
-
交换律
\[f*g=g*f \] -
结合律
\[(f*g)*h=f*(g*h) \]-
证明
注意到
\[((f*g)*h)(n) \]\[=\sum \limits_{abc=n} f(a)g(b)h(c) \]\[= \sum \limits_{a|n} f(a) \sum \limits_{b| \dfrac na} g(b)h(\dfrac {n}{ab}) \]\[= \sum \limits_{a|n} f(a) (g*h)(\dfrac na) \]\[=(f*(g*h))(n) \]
-
-
分配律
\((f\pm g)*h=f*h \pm g*h\)
-
等式性质
\[f*h=g*h \Leftrightarrow f=g (h(1) \neq 0) \]-
证明(参考 OI-Wiki)
右推左显然。
证明左推右,我们先假设 \(f\neq g\),那么存在最小的 \(x\) 使得 \(f(x) \neq g(x)\)。设 \(r=(f-g)*h\),那么有
\[r(x)=\sum \limits_{d|x} (f-g)(d)\times h(\dfrac xd) \]\[=\sum \limits_{d|x} (f(d)-g(d))\times h(\dfrac xd) \]\[=(f(x)-g(x))\times h(1) \]\[\neq 0 \]因为 \(f*h=g*h\),所以 \(r=(f-g)*h=0\),矛盾,于是一定有 \(f=g\),证毕。
-
-
一个结论
两个积性函数 \(f,g\) 的狄利克雷卷积仍然是积性函数。
证明:设 \(a,b\) 互质,则
\[(f*g)(a) \times (f*g)(b) \]\[=\sum \limits_{d_1|a} f(a)g(\dfrac {a}{d_1}) \sum \limits_{d_2|b} f(b)g(\dfrac {b}{d_1}) \]因为 \(a\) 和 \(b\) 互质,因此原式可以写作
\[\sum \limits_{d|ab} f(ab)g(\dfrac{ab}{d}) \]\[=(f*g)(ab) \]证毕。
单位元与狄利克雷逆
对于任意数论函数,容易证明,\(\epsilon*f=f\),所以 \(\epsilon\) 是狄利克雷卷积中的单位元。
如果 \(f*g=\epsilon\),则称 \(f,g\) 互为狄利克雷逆。由等式性质可知,一个数论函数的狄利克雷逆唯一。
-
另一个结论
积性函数的狄利克雷逆也是积性函数。
证明可见 OI-Wiki,这里不再证明。
常见数论函数的狄利克雷卷积
-
\(\mu * I = \epsilon\)
证明:\(n=1\) 时显然成立。当 \(n>1\) 时,设 \(n\) 的唯一分解式为 \(\prod \limits_{n=1}^{m} p_i^{k_i}\),则
\[(\mu * I)(n)= \sum \limits_{d|n} \mu(d) \]\[=\mu(1)-\mu(p_1)-\mu(p_2)\cdots-\mu(p_m)+\mu(p_1p_2)+\mu(p_1p_3)... \]\[=1-\dbinom m1 + \dbinom m2-\dbinom m3 + ... \]\[=(1-1)^m \] -
\(I * I = \tau\)
-
\(I * \operatorname{id} = \sigma\)
-
\(\varphi * I = \operatorname {id}\)
即证 \(\sum \limits_{d|n} \varphi(d)=n\)。
观察到 \(1\) 到 \(n\) 中,与 \(n\) 的最大公约数恰好为 \(d\) 的数有 \(\varphi(\dfrac nd)\) 个。那么上面的式子就是枚举 \(d\) 求和。
-
\(\mu * \operatorname {id}=\varphi\)
在上一个式子两侧卷上 \(\mu\),得 \(\varphi * I * \mu = \operatorname{id} * \mu\),化简即为上式。
这同时也推导出了 \(\varphi\) 是积性函数。
莫比乌斯反演
对于数论函数 \(f,g\),
证明:由左侧可知 \(f=g*I\),两边同时卷上 \(\mu\),得 \(g=f * \mu\)。
接下来看一些例题,如无说明均默认对某一大质数取模。
例题 1 Luogu P2398 GCD SUM
题意
求
\(1 \leq n \leq 10^5\)
题解
枚举 \(d\),\(i,j\) 必须是 \(d\) 的倍数才有贡献,于是有 \(\left \lfloor \dfrac nd \right \rfloor\) 个 \(i\) 和 \(j\) 符合条件,即
线性筛 \(\varphi\) 的前缀和,然后整除分块即可。时间复杂度 \(O(n)\)。
例题 2
题意
求
\(10^4\) 次询问,每次询问 \(1 \leq n \leq m \leq 10^7\)
题解
做法和例题 1 类似。
例题 3 Luogu P3455 ZAP-Queries
题意
求
\(5 \times 10^4\) 组询问,\(1 \leq n,m \leq 5 \times 10^4\)
题解
做法和例题 2 类似。
例题 4 [HAOI2011]Problem b
题意
求
\(5 \times 10^4\) 组询问,\(1 \leq a \leq b \leq 5 \times 10^4,1 \leq c \leq d \leq 5 \times 10^4,1 \leq k \leq 5 \times 10^4\)
题解
令 \(S(n,m)\) 表示 \(\sum \limits_{i=1}^{n} \sum \limits_{j=1}^{m} [\gcd(i,j)=k]\) 的答案,这就是例题 3。
可以容斥一下,最终的答案就是 \(S(b,d)-S(a-1,d)-S(b,c-1)+S(a-1,a-1)\)。
例题 5 [国家集训队]Crash的数字表格
题意
求
\(1 \leq n,m \leq 10^7\)
题解
设 \(S(n,m)= \sum \limits_{i=1}^{n} i \sum \limits_{j=1}^{m}j \cdot [\gcd(i,j)=1]\)。考虑化简 \(S\),
可以整除分块求解 \(S(n,m)\)。
此时原式
继续整除分块。
时间复杂度 \(O(n)\)。
例题 6 Luogu P2257 YY 的 GCD
题意
求
\(10^4\) 组询问,\(1 \leq n,m \leq 10^7\)
题解
发现好像化不动了。我们发现瓶颈似乎在于 \(\sum \limits_{p \in \mathbb P}\) 这个求和下标不连续,于是我们把求和顺序交换一下,原式
设 \(k=pd\),则原式
于是这个求和下标就连续起来了。设 \(f(x)= \sum \limits_{p \in \mathbb P \and p \mid x} \mu(\dfrac xp)\),只要求出 \(f\) 的前缀和就可以整除分块解决。
来尝试线性筛一下 \(f(i)\):
-
\(i\) 是质数
只有 \(p=i\) 有贡献,\(f(i)=\mu(1)=1\)。
-
枚举到 \(x=i \times y\),其中 \(y\) 是 \(x\) 的最小质因子
-
\(i\) 中不包含 \(y\) 这个质因子
此时 \(f(x)\) 里面多了一项 \(\mu(\dfrac xy)= \mu(i)\),并且原来的每一个 \(\mu(\dfrac {i}{p})\) 都变成了 \(\mu(\dfrac {iy}{p})\),多了一个质因子,于是 \(\mu\) 值会变为相反数,则有 \(f(x)=-f(i)+\mu(i)\)。
-
\(i\) 中包含 \(y\) 这个质因子
如果 \(i\) 中没有多个质因子,那么当 \(p=y\) 时 \(\mu(\dfrac xp)=\mu(i)\) 不为 \(0\),其余时候都为 \(0\),则有 \(f(x)=\mu(i)\)。
如果 \(i\) 中也有多个质因子,那么无论什么时候 \(\mu(\dfrac xp)\) 都为 \(0\),则有 \(f(x)=0\)。观察到这个时候有 \(\mu(i)=0\),于是 \(f(x)=\mu(i)\)。可以通过这个手段来简化代码。
-
综上,我们得到了 \(f\) 的计算式
此时原式
可以整除分块解决。预处理时间复杂度 \(O(n)\),单词询问 \(O(\sqrt n)\)。
例题 7 HDU5382 GCD?LCM!
题意
设
求 \(\sum \limits_{i=1}^n f(i)\)。
\(1 \leq n \leq 10^6\)
题解
我们发现不等式很令人头疼,于是我们找到了另外一个函数
观察到 \(f(n)=f(n-1)+2n-1-g(n-1)\)。其中 \(2n-1\) 就是 \((1,n),(2,n),\cdots,(n,n),(n,1),(n,2),...\) 这些有贡献的 \(i,j\)。观察到 \((n,n)\) 算重了,于是要减掉。
我们来化简 \(g\):
令 \(t(x)= \sum \limits_{i} \sum \limits_{j} [ij=d][\gcd(i,j)=1]\),那么 \(t(x)=2^m\),其中 \(m\) 是 \(x\) 不同的质因子个数,因为同一种质数不能同时分到 \(i,j\)。
我们线性筛预处理 \(t\),并枚举 \(d\) 和 \(k\),将 \(t(\dfrac {kd}{d}-1)=t(k-1)\) 的贡献加到 \(g(kd)\) 中,就得到了 \(g\)。按照递推式便可求出 \(f\)。时间复杂度 \(O(n \log n)\)。
例题 8 [SDOI2014]数表
题意
求
其中 \(\sigma(x) = \sum \limits_{d\mid x} d\)。
\(2 \times10^4\) 组询问,\(1 \leq n,m \leq 10^5, |a| \leq 10^9\)
题解
式子里出现了 \(x\) 和 \(d\) 的乘积。很套路地把它换成 \(T=xd\),然后枚举 \(T\),原式
整理一下,原式
也就是说只有 \(\sigma(d) \leq a\) 的 \(d\) 才对答案有贡献。
令
,则原式
观察到随着 \(a\) 的增大,有贡献的 \(d\) 越来越多,我们可以离线处理,将 \(a\) 排序,然后动态维护 \(S\)。
每加入一个 \(a\),我们都更新 \(a' < d \leq \min\{a,n\}\) 的 \(d\) 的贡献,其中 \(a'\) 是上一次加入的 \(a\)。具体来讲,加入 \(d\) 的贡献时将所有 \(S(kd)\) 加上 \(\sigma(d)\mu(\dfrac {kd}d) = \sigma(d)\mu(k)\),\(k\) 满足 \(kd \leq n\)。
此时我们完成了 \(a\) 限制的处理。因为询问需要整除分块处理,所以我们需要让 \(S\) 支持单点修改区间求和。树状数组就可以维护。
最多加入 \(n\) 次 \(d\) 的贡献,每次需要在树状数组上修改,这一部分的时间复杂度为 \((n + \dfrac n2 + \dfrac n3 + \cdots + 1 ) \log n = O(n \log^2 n)\)。
每次查询需要整除分块并树状数组查询,复杂度 \(O(\sqrt n \log n)\)。
总时间复杂度 \(O(n \log^2 n + q \sqrt n \log n)\),其中 \(q\) 是询问次数。
# include <bits/stdc++.h>
const int N=100010,INF=0x3f3f3f3f;
typedef long long ll;
const ll MOD=2147483648ll;
struct Node{
int id;
ll val;
bool operator < (const Node &rhs) const{
return val<rhs.val;
}
}sigs[N];
struct Queries{
int n,m,a,id;
bool operator < (const Queries &rhs) const{
return a<rhs.a;
}
}qu[N];
ll sigma[N],mu[N];
bool vis[N];
int prime[N],psum;
int T;
ll ans[N];
int pt;
ll sum[N];
inline int read(void){
int res,f=1;
char c;
while((c=getchar())<'0'||c>'9')
if(c=='-')f=-1;
res=c-48;
while((c=getchar())>='0'&&c<='9')
res=res*10+c-48;
return res*f;
}
inline void init(void){
mu[1]=1;
sigma[1]=1;
for(int i=2;i<=N-10;++i){
if(!vis[i]){
prime[++psum]=i;
mu[i]=-1,sigma[i]=i+1;
}
for(int j=1;i*prime[j]<=N-10&&j<=psum;++j){
vis[i*prime[j]]=true;
if(i%prime[j]==0){
sigma[i*prime[j]]=sigma[i]+prime[j]*(sigma[i]-sigma[i/prime[j]]);
break;
}else{
sigma[i*prime[j]]=sigma[i]*(prime[j]+1);
mu[i*prime[j]]=-mu[i];
}
}
}
for(int i=1;i<=N-10;++i){
sigs[i].id=i,sigs[i].val=sigma[i];
}
std::sort(sigs+1,sigs+1+N-10);
}
inline int lowbit(int x){
return x&(-x);
}
inline void bitchange(int x,int v){
for(;x<=N-10;x+=lowbit(x)){
sum[x]=(sum[x]+v)%MOD;
}
return;
}
inline ll bitquery(int x){
ll res=0;
for(;x;x-=lowbit(x)){
res=(res+sum[x])%MOD;
}
return res;
}
inline void change(int x){
for(int p=1;x*p<=N-10;++p){
bitchange(x*p,sigma[x]*mu[p]);
}
return;
}
inline ll query(ll n,ll m){
ll res=0;
if(n>m){
std::swap(n,m);
}
for(ll l=1,r;l<=n;l=r+1){
r=std::min(n/(n/l),m/(m/l));
res=(res+(n/l)*(m/l)%MOD*(bitquery(r)-bitquery(l-1))%MOD+MOD)%MOD;
}
return res;
}
int main(void){
// freopen("in.txt","r",stdin);
// freopen("out.txt","w",stdout);
init();
T=read();
for(int i=1;i<=T;++i){
qu[i].n=read(),qu[i].m=read(),qu[i].a=read(),qu[i].id=i;
}
std::sort(qu+1,qu+1+T);
pt=1;
for(int i=1;i<=T;++i){
while(pt<=N-10&&sigs[pt].val<=qu[i].a){
change(sigs[pt].id),++pt;
}
ans[qu[i].id]=query(1ll*qu[i].n,1ll*qu[i].m);
}
for(int i=1;i<=T;++i){
printf("%lld\n",ans[i]);
}
return 0;
}
例题 9 [CQOI2015]选数
题意
在区间 \([L,H]\) 中选取 \(N\) 个整数(一个整数可以被选多次),求它们的最大公约数为 \(K\) 的方案数。
两个方案是不同的当且仅当至少存在一个满足 \(1 \leq i \leq N\) 的 \(i\),使得两个方案选取的第 \(i\) 个整数不同。
\(1 \leq N,K \leq 10^9,1 \leq L \leq H \leq 10^9,H-L \leq 10^5\)
题解
我们观察到原式等价于
注意到 \(\dfrac{L}{K}\) 是上取整,因为可能会出现 \(a \times \left \lfloor \dfrac ba\right \rfloor \leq b\) 的情况,比如 \(a=2,b=3\),这显然不是我们所需要的。
继续化简:
设 \(l=\left \lceil \dfrac{L}{K} \right \rceil,r=\left \lfloor \dfrac{H}{K} \right \rfloor\),原式
\(=\sum \limits_{d = 1}^{r} \mu(d)(\left \lfloor \dfrac{r}{d} \right \rfloor - \left \lfloor \dfrac{l-1}{d} \right \rfloor)^N\)
因为 \(H \leq 10^9\),所以需要快速求出 \(\mu\) 的前缀和,然后整除分块。杜教筛即可。关于这种算法,详见本系列的第六篇《数学 6 杜教筛》。
# include <bits/stdc++.h>
typedef long long ll;
const ll MOD=1e9+7;
const int N=1000010,MAXN=1000000;
int mu[N];
int prime[N],psum;
bool vis[N];
std::unordered_map <int,int> camu;
int n,k,L,R;
inline int read(void){
int res,f=1;
char c;
while((c=getchar())<'0'||c>'9')
if(c=='-')f=-1;
res=c-48;
while((c=getchar())>='0'&&c<='9')
res=res*10+c-48;
return res*f;
}
inline ll qpow(ll d,ll p){
ll ans=1;
while(p){
if(p&1ll){
ans=ans*d%MOD;
}
d=d*d%MOD,p>>=1ll;
}
return ans;
}
inline void init(void){
mu[1]=1;
for(int i=2;i<=MAXN;++i){
if(!vis[i]){
prime[++psum]=i;
mu[i]=-1;
}
for(int j=1;j<=psum&&i*prime[j]<=MAXN;++j){
vis[i*prime[j]]=true;
if(i%prime[j]==0){
break;
}else{
mu[i*prime[j]]=-mu[i];
}
}
}
for(int i=1;i<=MAXN;++i){
mu[i]+=mu[i-1];
}
return;
}
inline int getmu(int x){
if(x<=MAXN){
return mu[x];
}
if(camu[x]){
return camu[x];
}
ll res=1ll;
for(int l=2,r;l<=x;l=r+1){
r=x/(x/l);
res-=getmu(x/l)*(r-l+1)%MOD;
}
return camu[x]=(int)(res);
}
int main(void){
init();
n=read(),k=read(),L=read(),R=read();
L=L/k+((L%k)>0);
R/=k;
ll ans=0;
for(int l=1,r;l<=R;l=r+1){
if((L-1)/l){
r=std::min((L-1)/((L-1)/l),R/(R/l));
}else{
r=R/(R/l);
}
ans=(ans+1ll*(getmu(r)-getmu(l-1))*qpow(R/l-(L-1)/l,n))%MOD;
}
// printf("%d ",getmu(19260817));
printf("%lld",(ans+MOD)%MOD);
return 0;
}
例题 10 [SDOI2015]约数个数和
题意
给定 \(n,m\),求
其中 \(d(x)\) 代表 \(x\) 的约数个数,和 \(\tau(x)\) 的定义相同。
多组数据,数据组数 \(T \leq 5 \times 10^4\),\(1 \leq n,m \leq 5 \times 10^4\)
题解
对于 \(ij\) 形式的式子我们不是很好处理。考虑转化:
证明:记 \(P(x,y)\) 表示命题「\(x \mid i\) 且 \(y \mid j\) 且 \(\gcd(x,y)=1\)」。只需要证明,满足 \(P(x,y)\) 的有序对 \((x,y)\) 和 \(ij\) 的约数一一对应。
设 \(i = \prod \limits_{l=1}^{r} p_l^{k_l}\),\(ij\) 的某个约数 \(t\) 的唯一分解式为 \(\prod \limits_{l=1}^{r} {p}_l^{{k'}_l}\)。我们可以构造一组 \((x,y)\) 与之对应。最初,\(x,y\) 均为 \(1\)。对于每一个 \(p_l\),如果 \({k'}_l \leq k_l\),\(x\) 乘上 \(p_l^{{k'}_l}\)。否则 \(y\) 乘上 \(p_l^{k'_l - k_l}\)。
容易证明,\(x,y\) 互质。观察到 \(t\) 是 \(ij\) 的约数,因此不存在某个 \(l\) 满足 \({k'}_l-k_l\) 大于 \(j\) 中 \(p_l\) 的指数。这样,\(y\) 总是 \(j\) 的因子。注意到 \(x\) 也是 \(i\) 的因子,这样,构造出的有序对 \((x,y)\) 就满足 \(P(x,y)\)。按照这种构造方法,对于不同的 \(t\),\(P(x,y)\) 一定不同。这样就证明了 \(t\) 到满足 \(P(x,y)\) 的 \((x,y)\) 的映射是一对一的。
接下来证明,对于任意满足 \(P(x,y)\) 的有序对 \((x,y)\),存在一个 \(t\),使得 \(t\) 可以构造出 \((x,y)\)。这样的 \(t\) 总是存在的。设 \(x=\prod \limits_{l=1}^{r} p_l^{kx_l}\),\(y = \prod \limits_{l=1}^{r} p_l^{ky_i}\)。对于每一个 \(l\),\(kx_l\) 和 \(ky_i\) 至多有一个大于 \(0\)。如果 \(kx_l\) 大于 \(0\),将 \(t\) 的 \(p_l\) 指数设为 \(kx_l\);否则,将 \(t\) 的指数设置为 \(k_l + ky_i\)。容易证明,\(t\) 是 \(ij\) 的约数。这样就证明了 \(t\) 到满足 \(P(x,y)\) 的 \((x,y)\) 的映射是映上的。
至此,\(t\) 到 \((x,y)\) 的映射既是一对一的,又是映上的,于是 \(t\) 和满足 \(P(x,y)\) 的 \((x,y)\) 一一对应。证毕。
那么,原式
设
是标准的整除分块形式。可以用 \(O(n\sqrt n)\) 的时间处理出所有 \(S(x)\),此时原式
这仍然可以整除分块处理。
总时间复杂度 \(O(\max \{n\} \sqrt{\max\{n\}} + T \sqrt n)\)。