Loading [MathJax]/jax/output/HTML-CSS/jax.js

「从零开始的莫比乌斯反演 & 杜教筛」

无聊开个坑……为了学弟的数论提前做好准备。
本文的写作目的是在你只接触过一点数论的情况下教会你莫比乌斯反演和杜教筛……
顺便加深一下我的理解。

下文规定,[P] 表示,若 P 成立,[P]=1,否则 [P]=0

数论分块

之所以要放到这里讲,是因为……不知道。

数论分块可以在 O(n) 的复杂度内计算 ni=1f(i)g(ni) 的值,但前提是可以求出 f(i) 的前缀和并算出所有 g(ni)(1in)
简单地说,就是,打表可以发现,ni 有很多连续的相同段,故可以用乘法分配律把相同的段放在一起算,那么复杂度应该就是段数。

然后有一个结论,ni(1in) 的取值只有 O(n) 种(2n 种左右)。
证明:
分类讨论。
in,则此时的 i 只有 O(n) 种,ni 自然也只有 O(n) 种。
i>n,则易得 ni<n,此时 ni 也只有 O(n) 种。
故结论成立。

于是,如果要实现程序,需要知道每一段的最右端,即在给定 i 的情况下求得满足 ni=nx 的最大的正整数 x
又一个结论:x=nni
证明:
由于 ni=nx,故 ninx<ni+1
故有 nni+1<xnni
x 最大,即为 nni
证毕。

参考代码

1
2
3
4
5
for(int l = 1,r;l <= n;l = r + 1)
{
r = n / (n / l);
//do sth with [l,r]...
}

筛法

埃氏筛

埃氏筛的起源是,首先假设 [2,n] 以内的整数都是质数,然后对于每个数,枚举其除了自身以外的倍数,并把其标记为合数。
最后得到的就是质数了。

然而这样子非常的慢……
根据调和级数,复杂度是 O(ni=1ni)=O(nlogn)

于是,埃氏筛就出现了:只枚举素数的倍数。
这个比较显然,对于合数 i,其倍数 d,也是 i 的因数的倍数,故合数的倍数均已经被筛过。

以上就是埃氏筛的主要思想,但是实际实现的时候,我们一般枚举 i 的倍数的时候,会从 i2 开始筛,因为其以下的倍数都已经被筛过了。

时间复杂度

埃氏筛的复杂度为 O(nloglogn)

证明如下(可以跳过):
i 个素数可以近似地表示为 ilni
π(n)n 以内的素数个数)也可以近似地表示为 nlnn
故算法的复杂度为 O(nlnni=2nilni)=O(nlnn2nxlnxdx)=O(nlnlnn)=O(nloglogn) 实际上这个复杂度并没有考虑到从平方开始枚举倍数的优化。

欧拉筛 / 线性筛

从名字就能看出来这个算法是线性的。

注意到埃氏筛优化之后仍然会有重复筛到的数。
于是欧拉筛的思想就是:令每个数被其最小质因子筛到。
这样就有了代码实现的问题。

于是,我们原来是枚举质数 i,然后枚举其倍数 ij(1jni)
现在我们换个顺序,先枚举这个 j,然后再枚举一个质数 i(1ini)

那么,令 j 的最小质因子为 p,则对于 ipi 显然都是 ij 的最小质因子。
而对于 i>p,则 ij 的最小质因子为 p,故此时可以退出。

参考代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#include <cstdio>
const int N = 1e7;
int n;
int vis[N + 5],cnt,prime[N + 5];
int main()
{
scanf("%d",&n);
for(int i = 2;i <= n;++i)
{
if(!vis[i])
prime[++cnt] = i;
for(int j = 1;j <= cnt && i * prime[j] <= n;++j)
{
vis[i * prime[j]] = 1;
if(i % prime[j] == 0)
break;
}
}
for(int i = 1;i <= cnt;++i)
printf("%d ",prime[i]);
}

上述代码输出了给定的正整数 n 以内的所有质数,其中 n107

时间复杂度

从代码似乎不能明显地看出复杂度,然而根据其思想,容易得到其复杂度为 O(n)

数论函数

定义

数论函数指定义域为正整数,值域为复数的函数。
在一些文献里也叫做算术函数。

积性函数

本来还有个加性函数,这里略过。

对于一个数论函数 f,若对于所有 gcd(x,y)=1 都有 f(xy)=f(x)f(y) 成立,则称 f 是一个积性函数。
特别地,若 f 对于所有 x,y 都满足上式,则称 f 为完全积性函数。

易知积性函数 f 一定满足 f(1)=1

常见的积性函数有:

  • 元函数:ϵ(n)=[n=1]
  • 单位函数:1(n)=1
  • 标号函数:ID(n)=n
  • 约数函数:σk(n)=dndk
  • 欧拉函数:φ(n)=ni=1[gcd(i,n)=1]
  • 莫比乌斯函数:μ(n)={1,n=10,d>1 d2n(1)ω(n),d>1 d2n,其中 ω(n) 表示 n 的不同质因子个数。
  • 刘维尔函数:λ(n)=mi=1(1)ci,其中 n=mi=1pcii

(比较显然的)性质

f,g 为积性函数,则:

  • h(x)=f(xC)
  • h(x)=fC(x)
  • h(x)=f(x)g(x)
  • h(x)=dxf(x)g(xd)

都是积性函数。

欧拉函数

定义:φ(n)=ni=1[gcd(i,n)=1]

通式

对于 n>1,令 n=mi=1pcii,其中 pi 为互不相同的质数(即分解质因数)。
φ(n)=nmi=1(11pi)

证明:考虑容斥。
首先假设前 n 个正整数都与 n 互质,然后从中减去 p1 的倍数,p2 的倍数……pm 的倍数。
但是这样会重复减 p1,p2 的公倍数,p1,p3 的公倍数……pm1,pm 的公倍数,于是再加回来。
然而又会重复加 p1,p2,p3 的公倍数,p1,p2,p4 的公倍数……pm2,pm1,pm 的公倍数,于是再减去。
……

注意到前 n 个正整数中 i 的倍数的个数为 ni,而以上的 i 都为 n 的因子,取整符号可去掉。
φ(n)=nnp1np2npm+np1p2+np1p3++npm1pmnp1p2p3np1p2p4npm2pm1pm

然后,就是欧拉函数通式推导的神奇的一步。
注意到每一项的符号取决于分母是多少个质因子的乘积,故可以视作是若干个 1pi 的乘积乘 n 的结果。
φ(n)=nmi=1(11pi)

由上式易得欧拉函数为积性函数。

性质及证明

欧拉定理

gcd(a,n)=1,则有 aφ(n)1(modn)

严谨的证明要牵扯到各种剩余系之类的名词,这里用一种比较通俗的说法代替。
证明:
x1,x2,,xφ(n) 为前 n 个正整数中与 n 互质的 φ(n) 个数。
再令 yi=axi

这里要先证明两个性质:

  1. yiyj(modn)(ij)
    反证法。假设存在 yiyj(modn)
    则根据同余式的性质与 gcd(a,n)=1,有 xixj0(modn)
    但由于 1xi<n,易知假设不成立。

  2. gcd(yimodn,n)=1
    由于 yi=axi,又因为 gcd(a,n)=gcd(xi,n)=1,又根据欧几里得算法,原式显然成立。

根据上述两个性质可得 φ(n)i=1yiφ(n)i=1xi(modn)φ(n)i=1axiφ(n)i=1xi(modn)aφ(n)φ(n)i=1xiφ(n)i=1xi(modn)aφ(n)1(modn)

证毕。

n 为质数时,φ(n)=n1,就等价于费马小定理。
由欧拉定理可以导出 ababmodφ(n)(modn)(gcd(a,n)=1)

扩展欧拉定理

对于任意的 a,nbφ(n),有 ababmodφ(n)+φ(n)(modn)

证明:
n=pcs,其中 p 为一个质数,gcd(pc,s)=1(即提取 n 的一个质因子)。
根据欧拉定理,有 pφ(s)1(mods)
而由于欧拉函数具有积性,可得 φ(s)φ(n),故 pφ(s)pφ(n)(modn)
pφ(s)=xs+1,则 pφ(s)+c=xn+pc
pφ(s)+cpc(modn)
根据以上过程同样可证 pφ(n)+cpc(modn)
bc 时,pcpbcpφ(n)+cpb+φ(n)(modn)
而由于 cφ(n),当 b2φ(n) 时,pbpbφ(n)(modn)
pbpbmodφ(n)+φ(n)(modn)

对于 a 的一个质因子 p,当 pn,根据以上过程可得 pbpbmodφ(n)+φ(n)(modn) 成立;当 pn,根据欧拉定理同样可得上式成立。
故原式成立。
证毕。

其他(比较显然的)性质

  • φ(2n)=φ(n)(n1(mod2))
  • φ(pc)=pcpc1,其中 p 为质数。
  • φ(np)={φ(n)p,pnφ(n)(p1),pn,其中 p 为质数。
    由通式显然易得。
    由这条性质可以得出线性筛欧拉函数的方法。
  • 对于 n>2φ(n)0(modn)
    由更相减损术(gcd(n,i)=gcd(n,ni)(n>i))可得。
  • ni=1[gcd(i,n)=1]i=φ(n)n+[n=1]2
    同样由更相减损术可得,对于 n>1,与 n 互质的数是成对出现的。
  • dnφ(d)=n
    证明:
    将前 n 个正整数按与 ngcd 分类,则由于每个数与 ngcd 是唯一的,有 n=dnni=1[gcd(i,n)=d]
    ni=1[gcd(i,n)=d]=di[gcd(i,n)=d]=ndi=1[gcd(i,n)=1]=φ(nd)
    其中用到了 gcd(agcd(a,b),bgcd(a,b))=1 这个性质,用反证法易证。

欧拉反演

虽然是个假名字,但还挺有用的。
即利用 dnφ(d)=n 这条性质。
之所以叫欧拉反演,因为它可以用来做一些莫反的题,但局限性很强。

例:求 ni=1nj=1gcd(i,j)
则以 dnφ(d)=n 代入其中的 gcd(i,j) 项可得 ni=1nj=1dgcd(i,j)φ(d)
接下来是反演题推导的一个常用套路:交换求和号。
一般的策略即在枚举约数时,考虑枚举约数,再考虑其倍数与其产生的贡献。
nd=1φ(d)ni=1nj=1[dgcd(i,j)]
注意到 dgcd(i,j) 等价于 di,dj,并且 i,j 的枚举之间互不影响,由乘法原理可得 ni=1nj=1[dgcd(i,j)]=(nd)2
nd=1φ(d)ni=1nj=1[dgcd(i,j)]=nd=1φ(d)(nd)2
看起来这个式子依然要 O(n) 来做。
但是,如果利用上文中提到的数论分块的技巧(f(n)=φ(n),g(n)=n2),可以把复杂度降到 O(n)(前提是线性筛出了欧拉函数的前缀和)。

莫比乌斯函数

莫比乌斯函数的定义式为 μ(n)={1,n=10, d>1,d2n(1)ω(n),otherwise,其中 ω(n) 表示 n 的不同质因子个数。
有的文献中也把定义式写作:设 n=mi=1pcii,则 μ(n)={1,n=10, ci>1(1)m, ci=1
其实两者是等价的。

一个重要的性质

dnμ(d)=[n=1]
其实类似这种东西的证明可以根据积性,分解质因数后证明对于质数的幂成立,并证明原式。
这样子比较套路,这里有另外一种证明。

对于 n=1,原式显然成立。
对于 n>1,注意到若 d 有完全平方数因子(简称为平方因子),则 μ(d)=0
故有贡献的 d 肯定得在 n 的所有质因子中,每个至多取一个,然后乘起来得到。
此时 μ(d)=(1)ω(d)
dnμ(d)=ω(n)i=1Ciω(n)(1)i
由二项式定理((a+b)k=ki=1Cikaibki)可知 ω(n)i=1Ciω(n)(1)i=ω(n)i=1Ciω(n)1ω(n)i(1)i=(11)ω(n)=0
证毕。

为什么称它为一个重要的性质呢,因为莫比乌斯反演的证明需要用到这个性质,所以它可以看作一个引理。
但这个引理的用途不止这点。

莫比乌斯反演

若有数论函数 f,F 满足 F(n)=dnf(d),则 f(n)=dnμ(nd)F(d)
以上即莫比乌斯反演公式。

证明:
F(n)=dnf(d) 代入上式,可得 dnμ(nd)kdf(k)
考虑将 k 提到最前面,即 knf(k)kd,dnμ(nd)
d=dk,则由 dkn,kn 可得 dnk
knf(k)kd,dnμ(nd)=knf(k)dnkμ(ndk)
注意到 dnkμ(ndk)=dnkμ(d)=[nk=1]=[n=k]
故原式成立。

不过一般会用到的是莫比乌斯反演公式的另一种形式:若 F(n)=ndf(d),则 f(n)=ndμ(dn)F(d)
证明方式类似。
其实在学到下文中的数论函数的狄利克雷卷积之后,莫比乌斯反演公式只需要几步即可证明。

然而实际上,通常情况下推导式子的时候使用莫比乌斯反演公式的话会非常麻烦,所以使用莫比乌斯函数的那个重要的性质足矣。

例:求 ni=1mj=1[gcd(i,j)=k]
首先 i,j 必为 k 的倍数,故 kikj[gcd(i,j)=k]
同样用到 gcd(agcd(a,b),bgcd(a,b))=1 的性质,得 nki=1nkj=1[gcd(i,j)=1]
dnμ(d)=[n=1] 代入 [gcd(i,j)=1],得 nki=1nkj=1dgcd(i,j)μ(d)
d 提到最前面,得 min(nk,mk)d=1μ(d)ndkmdk
做到这一步,就可以使用数论分块的技巧求解了,复杂度是 O(nk+mk)
你可能会问两个参数该如何数论分块,实际上,只需要每次令 r=min(nknkl,mkmkl) 即可。
此题即「『POI2007』Queries」。

如果你有在认真读,可能会注意到我把 nkd 换成了 ndk(另外一个也是)。
证明:
nk=nk+c,则 ndk=nkd+cd
注意到 nkdnkd+d1dcd<1d
ndk=nkd+cd=nkd
证毕。

掌握套路之后,不妨再做点题。
例(「『SDOI2015』约数个数和」):求 ni=1mj=1σ0(ij)
这里需要用到一个性质:σ0(ij)=xiyj[gcd(x,y)=1]
这个性质在我当时做这道题的时候并不会证明,当时写的题解中也没有写。
现在来证一下。

对于一个 dij,令 x=igcd(i,d),y=dgcd(i,d)
则显然易得 gcd(x,y)=1,xi
ygcd(i,d)ij 可得 yxj,又因为 gcd(x,y)=1,故 yj
于是可得 gcd(x,y)=1,xi,yj
而对于任意的 gcd(x,y)=1,xi,yj 都可以得到 iyxij,并且与所有的 dij 一一对应。
故原式得证。

然后接着来推导,这里就直接略去文字注释了: ni=1mj=1σ0(ij)=ni=1mj=1xiyj[gcd(x,y)=1]=ni=1mj=1xiyjdgcd(x,y)μ(d)=min(n,m)d=1μ(d)dxdyxiyi1=min(n,m)d=1μ(d)dxdynxmy=min(n,m)d=1μ(d)ndx=1mdy=1ndxmdy=min(n,m)d=1μ(d)f(nd)f(md)

之所以把后面的东西写成函数,是因为这两个式子只和 nd,md 有关。
易知 f(n)=ni=1ni
这题数据范围比较小,所以说直接暴力预处理 O(nn) 是可以过的(视 n,m 同阶)。

但是如果数据范围开到 108(不考虑多组数据)?
不做预处理,直接在用到的时候暴力跑,因为 f 函数只需要用到形如 ni 处的值。
复杂度是 O(ni=1(i+ni))=O(n1nxdx)=O(n3/4)
(实际上如果使用杜教筛求 μ 前缀和的话,可以做到总共 O(n2/3)

然而实际上 f 函数是可以使用线性筛求出的。
注意到 f(n)=ni=1ni=ni=1id1=nd=1σ0(d) 于是 f 其实就是 σ0 的前缀和。
线性筛出 σ0 即可。
具体方式是根据约数个数定理(σ0(n)=mi=1(ci+1)ci 是前文提到的质因子分解的形式),线性筛时记录最小质因子的次数即可(因为只会被最小质因子筛到)。

在下文中还会再提到此题的 O(n2/3) 做法。

更多的习题还可以在我博客的莫比乌斯反演标签中找到。

狄利克雷卷积

定义两个数论函数 f,g 的狄利克雷卷积 fg 仍然是一个数论函数,且 (fg)(n)=dnf(d)g(nd)
狄利克雷卷积满足交换律 fg=gf,结合律 (fg)h=f(gh)
此处证明略去。

根据上文可得 μ1=ϵ,φ1=ID,fϵ=f

上文中有提到莫比乌斯反演公式的狄利克雷卷积证明,即证明 F=f1f=μF
证明:将 F 代入后式,得 f=μ(f1)
根据交换律与结合律得 f=(μ1)f=ϵf=f,故得证。

数论函数前缀和

我们常常会遇到数论函数的前缀和问题,尽管有时只是核心算法的一部分。
而很多时候,我们需要的是所有 ni 处的前缀和值。
然后呢,有一种常用的算法,可以用来对一些常见的数论函数解决这个问题。

杜教筛

对于积性函数 f,g,h=fg,考虑 ni=1h(i)
ni=1h(i)=ni=1dig(d)f(id)=nd=1g(d)dif(id)=nd=1g(d)ndi=1f(i)

S(n)=ni=1f(i),则原式即 nd=1g(d)S(nd)
从和式中提出 d=1 一项,得 ni=1h(i)=g(1)S(n)+nd=2g(d)S(nd)
注意到 g(1)=1,故 S(n)=ni=1h(i)nd=2g(d)S(nd)

若能较快地求出 g,h 的前缀和,即可递归地求 S(n)
记忆化之后,实际上就是对于所有 ni 花费了 O(ni) 的复杂度。
总复杂度即 O(ni=1i+ni)=O(n1nxdx)=O(n3/4) 若预处理 nc(c>12) 以内的值,则复杂度变为 O(nc+n1c2)
c=23 时两项平均达到最优值 O(n2/3)

另一种利用狄利克雷卷积的技巧

如果我们有 f=gh,而要求 ni=1f(i),则可以做得更简单一点。

ni=1f(i)=ni=1dig(d)h(id)=nd=1g(d)dih(id)=nd=1g(d)ndi=1h(i)

于是只需要求得 g,h 的前缀和即可。
用和杜教筛类似的预处理方式可以做到同样的 O(n2/3)

应用

  • μ 的前缀和:
    注意到 μ1=ϵ,使用杜教筛即可。
  • φ 的前缀和:
    注意到 φ1=ID,使用杜教筛即可。
  • σ0 的前缀和:
    注意到 σ0=11,使用另一种技巧即可。
  • (μID)1 的前缀和:
    首先使用另一种技巧转化为求 μID 的前缀和。
    然后注意到 (μID)ID=ϵ
    使用杜教筛即可。
    需要注意的是预处理可能会比较麻烦,不过可以在质数次幂处分类讨论。
    可参考 EI 的「复杂度分析:积性函数的狄利克雷卷积」
  • φID 的前缀和:
    直接使用另一种技巧即可。

扯淡

某些特殊的函数除了这些基本套路以外,还有别的做法。
比如求 μ2 的前缀和,可以通过 μ 的定义得到 O(n3/5) 的做法(如果只要求 n 处的前缀和,可以 O(n1/3))。
这里就不扯了,读者可以自己思考。

Related Issues not found

Please contact @Alpha1022 to initialize the comment