前情提要:某个上午,xtr 教练给我们恶补了一下数论。没错,这整篇文章就是他恶补的内容。
注:这里面大部分内容我都不给出证明,因为实在太多内容要写了/kk。
数论函数
数论函数就是输入一个正整数,输出一个整数的函数。
这里先介绍几个常见的数论函数:
- 值恒为 1 的 1(n)
- 除了 n=1 为 1 外其他位置为 0 的 ϵ 函数
- 欧拉函数 φ(n)
- 数的因子个数 σ(n)
- id(n)=n
Dirichlet 卷积
直接扔定义:
对于两个数论函数
f,g,以以下方式定义卷积:
(f∗g)(n)=i∣n∑f(i)g(in)
也可以容易地验证卷积这个运算满足交换律和结合律。
这个时候,如果我们知道
f,g 的某一个前缀,根据调和级数,直接计算
f∗g 的前缀是
O(nlogn) 的。
不难观察到
ϵ∗f=f,其中
f 是任意数论函数。于是发现
ϵ 是 Dirichlet 卷积的单位元。
最终还有一个非常重要的式子:
φ∗1=id,这个式子非常重要,在我们后面求
φ 会发挥很大的作用。
Dirichlet 逆
卷积的运算还可以有逆元,就像乘法有乘法逆元一样。
定义
g 是
f 的 Dirichlet 逆,当且仅当
g∗f=f∗g=ϵ,可以记
g=f−1。
根据 Dirichlet 卷积公式可以得到:
i∣n∑f−1(i)f(in)=f−1(n)f(1)+i∣n,i<n∑f−1(i)f(in)=[n=1]
移项得:
f−1(n)=f(1)1([n=1]−i∣n,i<n∑f−1(i)f(in))
不难发现,数论函数
f 存在 Dirichlet 逆的充要条件是
f(1)=0。
积性函数
积性函数是数论函数中的重要概念。一般都是讨论积性函数,可能是因为其他函数没有这样的性质,不好研究。
如果对于所有
(a,b)=1,都有
f(a)f(b)=f(ab) 的话,则称一个数论函数
f 是积性函数。
上面举得数论函数的例子中,除了
ϵ 函数以外均是线性的。
在积性函数的基础上,如果对于数论函数
f,还能够满足对于
(a,b)=1,
f(a)f(b)=f(ab) 也成立,则称
f 为
完全积性函数(
id 就是一个例子)。否则称
f 为不完全积性函数。
可以证明如果
f,g 是积性函数,则
f∗g 也是积性函数。
如果
f 是积性函数,那么
f−1 也是积性函数。
以上被称为积性函数的传递性。
注意:传递性只在不完全积性函数上面成立,在完全积性函数上面不成立。
一般来说,积性函数的 f(1) 都是 1,这个事情在以后的文章里面是默认的,如果有变动会特殊注明。
莫比乌斯函数
莫比乌斯这个人小时候就不陌生,莫比乌斯带也和这个人有关。但是莫比乌斯函数和莫比乌斯带没有关联。
通过 Dirichlet 逆,我们可以引出一个函数,被称为莫比乌斯函数
μ。它被定义为数论函数
1(n) 的 Dirichlet 逆,也就是
μ=1−1。
莫比乌斯函数的计算方法:
- 如果 n 含有 >1 的平方因子,则 μ(n)=0。
- 否则,设 k 为 n 含有的质因子个数,则 μ(n)=(−1)k。
另外的,就是还要牢记一个关于莫比乌斯函数的性质:
[gcd(i,j)=1]=d∣gcd(i,j)∑μ(d)
还有一个,可以通过前面的式子
φ∗1=id 得到。将这个式子的左右两边同时卷上
μ,可以得到
φ=id∗μ。
在做题的时候还发现了其他的一些式子:
d∣n∑dμ(d)=nφ(n)
i=1∑n[(i,n)=1]i=2nφ(n)+[n=1]
这个式子的证明见 https://www.luogu.com.cn/article/5k1yg69y 。
实在不想打了,求轻喷
莫比乌斯反演
莫比乌斯反演是指这两个神秘的公式:
-
若
gn=∑i∣nfi,则有
fn=∑i∣nμ(in)gi。
-
若
gn=∑n∣ifi,则有
fn=∑n∣iμ(ni)gi。
乍一看很复杂,实际上感觉挺合理的:你从用
f 表达
g 变成了用
g 来表达
f,很容易联想到逆元,而
1 和
μ 函数恰好就是逆元!
这两个公式被我们称为莫比乌斯反演公式。和二项式反演等科技类似,一般在求解
f 函数的值时,发现和
f 相关联的
g 函数的值更加容易计算,这个时候就可以考虑使用莫比乌斯反演,使用
g 来倒推
f 的值。
所以有时候把简单问题复杂化也是一个好方法。
整除
来讨论整除
⌊in⌋ 的一些基本性质:
-
⌊j⌊in⌋⌋=⌊ijn⌋,也就是先除以
i 下取整再除以
j 下取整和直接除以
ij 下取整的结果是一样的。这也被称为是整除点的传递性。
-
对于
i≤n,
⌊⌊in⌋n⌋ 是最大的使得
⌊xn⌋=⌊in⌋ 的
x,也就是说,会使得
⌊xn⌋ 发生变动(这里定义为
x→x+1 之后会改变)的位置
x 均形如
⌊⌊in⌋n⌋。
不妨记
Dn:={⌊in⌋∣i=1,2,⋯,n},表示
n 的所有整除点的集合。
- ∣Dn∣=O(n)
这可以通过根号分治得到,对于
i≤n,自然只有
O(n) 种数。否则
⌊in⌋≤n,仍然只有
O(n) 种数。
整除分块 / 数论分块
因为整除点集合的大小是根号级别的,因此在涉及到对整除点相关的式子进行求和的时候,可以使用整除分块对计算进行加速,可以将
O(n) 甚至更高的复杂度压到
O(n)。
欧拉定理和扩展欧拉定理
这里笔者想到了这两个东西,顺便就写一些。
欧拉定理:
- 设 a,p∈N+,(a,p)=1,则:
- aφ(p)≡1(mod p)。
这个定理非常常用,但就是有些局限性:
a,p 要是互质的。
扩展欧拉定理:
- 对于任意 b≥φ(p),有:ab≡abmodφ(p)+φ(p)(mod p)。
这个时候,
a,p 可以不互质,但是指数
b 必须要比
φ(p) 大。
P4139 上帝与集合的正确用法
这道题中,我们需要求这个东西:
222⋯modp
前面的被除数是一个无限层乘方塔,怎么办?
看到这种幂次还要模上一个数的式子,很容易想到使用费马小定理或者是欧拉定理。 显然费马小定理这个时候也无能为力,所以考虑使用欧拉定理。
但是这里
p 是会变化的,不能保证
a=2 和
p 一定互质。所以欧拉定理这里是有局限性的。
那怎么办呢?还有扩展欧拉定理!
这个时候虽然底数不太管用,那么就使用指数
b=222⋯。因为这个乘方塔有无限层,所以这个
b 也一定是一个无线层的乘方塔,所以
b 可以到达无穷大的级别,不管怎样一定会满足
b>φ(p)。
于是就变成了:
222⋯≡2222⋯modφ(p)+φ(p)(modp)
发现这个时候指数又变成了一个乘方塔模上一个数的形式,于是考虑递归。容易发现如果一直
p→φ(p) 的话,最终一定会变成
1,所以这个递归最终一定会结束。
线性筛求解积性函数
对于许多常见的积性函数,只需要在线性筛的基础上略施小计,改改就可以在线性复杂度内求积性函数的前几项。
假设对于质数
p,积性函数的值
f(p) 是很容易求出的(例如
μ(p)=−1,φ(p)=p−1)。
对于合数,发现一部分合数的值可以在线性筛将
i×pj 筛去的时候顺便求解。根据积性函数的性质,如果
pj∣i,则
f(i×pj)=f(i)f(pj),可以直接求解。
对于
pj∣i 的情况,这时候的
i×pj 是另一部分合数,这个时候的函数值需要具体问题具体分析。
例如要求解
μ 函数,如果
pj∣i,就说明
i×pj 的质因数分解里面至少有了
2 个
pj,这个时候
μ(i×pj) 显然为
0。
或者是要求解
φ 函数,如果
pj∣i,就说明
i×pj 的质因数分解里面至少有了
2 个
pj,联系
φ 函数的计算方法,这个时候
φ(i×pj)=φ(i)×pj。
顺便给出
μ 和
φ 函数的求解方式,使用线性筛。
CPPvoid get() {
ph[1] = mu[1] = phs[1] = mus[1] = 1;
for (int i = 2; i < N; i++) {
if (!f[i])
f[++f[0]] = i, ph[i] = i - 1, mu[i] = -1;
for (int j = 1; j <= f[0] && i * f[j] < N; j++) {
f[i * f[j]] = 1;
if (i % f[j])
mu[i * f[j]] = -mu[i], ph[i * f[j]] = (f[j] - 1) * ph[i];
else {
mu[i * f[j]] = 0, ph[i * f[j]] = ph[i] * f[j];
break;
}
}
phs[i] = phs[i - 1] + ph[i], mus[i] = mus[i - 1] + mu[i];
}
}
例题 1
我们来看一道例题。
求有多少个数对
(i,j) 满足
i∈[1,n],j∈[1,m],且
i,j 互质。有
T 组数据。
T≤1000,n,m≤105。
互质就是最大公约数
gcd=1,这一点好像不太方便批量计算。当然你如果跟我说
i,j∈[1,n] 的话我还可以使用欧拉函数来算,但是这道题不是这样的。
考虑把问题复杂化,考虑从求
i,j 最大公约数恰好为
1 变成求
gcd=k,将方案数记为
f(k)。答案就是
f(1),但是这里的
f 值好像都不好求。
这个时候会联想到:如果把问题弱化,将
恰好为 k 变成
为 k 的倍数 呢?这个时候显然很容易,因为这样只需要满足
k∣i,k∣j,答案就是
⌊kn⌋⌊km⌋。不妨就设
g(x)=⌊kn⌋⌊km⌋。
考虑
f,g 的联系,不难发现
g(k)=∑k∣df(d)。
注意到这就是
莫比乌斯反演中的一条公式,于是大力反演,可以发现
f(1)=∑kμ(k)⌊kn⌋⌊km⌋。
自然可以利用线性筛,
O(n) 预处理出所有的
μ。随后,根据整除分块的性质可以发现随着
k 的增大,
⌊kn⌋ 和
⌊km⌋ 分别变化不超过
O(W) 次,其乘积显然也是这样的。对于
⌊kn⌋⌊km⌋ 不变的一段区间,利用预处理计算
μ 的前缀和,区间和便也可以求出。
复杂度是
O(n+Tn) 的,可以通过。
发现
发现刚刚的题目中,通过
O(n) 提前处理
μ 函数的值,我们可以在根号的时间内回答一个询问。这种特性使得我们可以解决多组询问的问题。
但是如果我们将问题改为很少组询问,但
n 的范围放大很多,如
n≤1010,那么此时
O(n) 回答一个询问的代价仍然可以接受,但是线性复杂度的预处理就不太快了。于是我们的目标就是在于优化前面的预处理,或者是
减少前面预处理的范围,后面的一段现场计算。
可能我前面说的有些混乱,反正我们的目标就是:希望在低于
O(n) 的时间复杂度内,求出
n 的所有整除点在
Dn 处的
f 的前缀和。
可以先提前预测一波:因为至少有
O(n) 个整除点,所以我们的预期算法的时间显然应该在
O(n) 到
O(n) 之间。
这个预期算法就是:杜教筛。
杜教筛
杜教筛就是用来快速计算一些特殊的积性函数的前缀和的算法,由国家队爷 dyh 提出。
这种算法比线性筛预处理更快一些,不过其只能处理一些特殊的积性函数,适用范围有点小。
首先,我们先考虑 Dirichlet 卷积,设
h=f∗g,则
h 的前缀和根据定义可以表示为:
i=1∑nh(i)=i=1∑nj=1∑⌊in⌋f(i)g(j)=i=1∑nG(⌊in⌋)
其中
G 是
g 的前缀和函数,我们可以通过一些代数推导,改写这个式子:
G(n)=i=1∑nh(i)−i=2∑nf(i)G(⌊in⌋)
这个改写之后的式子恰好告诉了我们如何计算
g 的前缀和。
从右边的一半部分,同样可以由整除分块知识知道,我们需要递归求解所有
g 在
⌊in⌋ 处的前缀和,这恰好就是
n 的所有整除点的位置,所以和问题十分有相同点。
当然,我们必须要使用递归来逐层求解这个前缀和。而在递归求解时,我们还要求解
⌊in⌋ 的整除点的前缀和,不过根据整除点的传递性,这个点显然也是
n 的整除点。所以这个依赖关系对于
n 的全体整除点恰好是封闭的,也就是说要求解
n 的前缀和只需要求解
n 的所有整除点的前缀和即可。
同时为了利用整除分块,自然还需要
f,h 的前缀和。所以为了求解
f,h 的前缀和,
f,h 还得比较简单,这也是为什么只能处理一些特殊的积性函数。
因此这启发我们:假设数论函数
g 的前缀和不方便求解,但是我们可以找到方便快速求解前缀和的数论函数
f,h,满足
h=f∗g,则可以通过上面的式子计算
g 的前缀和。
我们来通过
μ 和
φ 函数练练手,你可以针对
g=μ 和
g=φ 找到合适的
f 和
h 吗?
回顾之前的数论函数的性质,我们知道:
{1∗μ=ϵ1∗φ=id
μ,ϵ,id 的前缀和都是非常显然的,因此
μ,φ 都可以通过上面的式子求出前缀和。
杜教筛复杂度
当然,为了确定杜教筛是不是我们最终的预期算法,我们还需要计算杜教筛的复杂度。
考虑前缀和的计算式子:
G(n)=i=1∑nh(i)−i=2∑nf(i)G(⌊in⌋)
可以知道,求解
n 的前缀和的整除分块计算次数是
O(n) 次。这个操作需要对
n 的所有整除点各进行一次,因此总的复杂度应该是
n 的所有整除点开方之后的和。
根据整除点
⌊in⌋ 的性质,我们可以将这个结果根据
i≤n 和
i>n 分别估计,最终就是这个式子:
i=1∑nin+i=1∑ni
这个式子乍一看很复杂,但是我们可以继续推导:
O(∑i=1nin+∑i=1ni)=O(∑i=1nin)=O(n∑i=1ni1)
这里,我们引入一个使用微积分可以证明的奇妙结论,这里就略过证明了:
O(∑i=1mik)=O(mk+1),其中
k 是
>−1 的实数。
这样,我们代入
m=n,k=−21,就可以得到:
O(ni=1∑ni1)=O(n×n)=O(n43)
最终,我们得到了杜教筛的初步版本的时间复杂度。
杜教筛优化
O(n43) 确实已经比
O(n) 快了不少,但是,我们还可以让他更快吗?
让我们重新考虑一下线性筛。对于在一个小范围
m 内的所有整除点的前缀和,
n 的整除点可能很密集,这个时候采用所有整除点的开方之和可能甚至还比线性的复杂度更劣。
可能你不太听得懂,但是没关系,举个例子就懂了。
不妨设
m=n32,此时,不超过
n32 的整除点的开方之和是:
O(n43)
而如果我们这个时候如果乖乖使用线性筛,最终时间复杂度只有
O(n32),这是更优的。
所以我们有一个更加聪明的做法:让我们钦定一个比较小的整数
m,在
m 以下的整数使用线性筛,而在
m 以上的整数使用整除分块。
恰好也是
m=O(n32),这个时候两部分的复杂度平衡,总复杂度也是
O(n32)。这就实现了完整的杜教筛算法。
而且可以通过不等式证明,这样的
m 是最优的,也就是我们最终的算法就止步于
O(n32),但是也很快了,有些接近根号。
Powerful Number 筛
虽然杜教筛的复杂度已经达到了很优秀的
O(n32),但是我们发现杜教筛有一个比较明显的短板:杜教筛对于可以求解的积性函数
g 的要求过高,需要满足
f∗g=h,其中
f,h 的前缀和还是需要可以直接快速求解这样的形式。
而在现实中,除了
μ,φ 这种特殊的积性函数,几乎不太可能对想要求解的积性函数
g 找到满足条件的
f,h。
那么,是否有什么方法能够在利用杜教筛简单的求解思想的同时,放宽对积性函数的要求呢?答案是有的,这就是 powerful number 筛,简称 PN 筛。
powerful number
PN 筛和 powerful number 是由很大关联的,那么 powerful number 是什么呢?先给定义:
- 如果一个数 n 的每个质因子的次数都至少为 2,则称 n 是一个 powerful number。
我们首先需要证明一个结论:
- n 以内的 powerful number 的个数是 Θ(n) 的。
我们考虑同时证明上界和下界。
因为所有的 powerful number 都可以表示为
a2b3 的这种形式,显然
n 以内,
a≤n,给定
a,b 的个数是
O((a2n)31) 的,因此总
a,b 个数是:
a=1∑nO((a2n)31)=O(a=1∑na32n31)=O(n31×a=1∑na321)=O(n)
(其中,倒数第二步到最后一步是使用的那个要用微积分证明的奇妙结论)
所以,powerful number 的个数至多是
O(n) 的。另一方面,因为平方数肯定都是 powerful number,因此 powerful number 的个数至少是
Ω(n) 的,因此 powerful number 的个数就是
Θ(n) 的。
PN 筛
接下来,我们考虑如何利用 powerful number 来扩展杜教筛的使用范围。
假设我们想要求解积性函数
f 的前缀和,但是
f 不满足杜教筛的适用条件。此时,如果我们能找到一个满足杜教筛的使用条件的积性函数
g,并且对于任意质数
p,满足:
f(p)=g(p)
则我们可以考虑
h=f∗g−1,尝试求解
h(p) 的值,可以发现由于
f(p)=g(p),我们有
h(1)g(p)+h(p)g(1)=g(p)+h(p)=f(p)+h(p)=f(p),也就是
h(p)=0。这样根据积性函数的性质可以发现,当且仅当
n 是 powerful number 时,
h(n) 才有可能是非
0 的。
因为
f=h∗g,此时
f 的前缀和可以表示为:
i=1∑nf(i)=i=1∑nj=1∑⌊in⌋h(i)g(j)=i=1∑nh(i)j=1∑⌊in⌋g(j)
由于
g 可以使用杜教筛的条件,可以在
O(n32) 内求解所有
n 的整除点处的前缀和。而
h 仅在 powerful number 处可能非
0,因此只需要枚举所有 powerful number 的位置并求解
h,就可以
O(n) 求解
f 的前缀和。