浅谈积性函数的线性筛法
前置知识
线性筛
线性筛可以在严格$O(n)$的时间内筛出积性函数的值，
它有常见的套路
假设$n = p_1^{a_1} p_2^{a_2} \dots p_k^{a_k}$
如果我们能快速得出$f(p_i)$和$f(p_i^{k+1})$的取值，那么直接套板子就行了
在下文中如无特殊说明，默认$p_i$表示$n$质因数分解之后第$i$个质数，$a_i$表示$p_i$的指数
常见的有以下几种
线性筛素数
比较简单，这也是筛其他积性函数的基础
#include<cstdio>
const int MAXN = 1e4 + 10;
int N, prime[MAXN], vis[MAXN], tot;
void GetPrime(int N) {
vis[1] = 1;
for(int i = 2; i <= N; i++) {
if(!vis[i]) prime[++tot] = i;
for(int j = 1; j <= tot && i * prime[j] <= N; j++) {
vis[i * prime[j]] = 1;
if(!(i % prime[j])) break;
}
}
}
int main() {
GetPrime(1e3);
for(int i = 1; i <= tot; i++)
printf("%d ", prime[i]);
return 0;
}线性筛莫比乌斯函数
这个也是比较常见的
根据莫比乌斯函数的定义
$$\mu =\begin{cases}\left( -1\right) ^{k}\left( n=p_{1}p_{2}\ldots p_{k}\right) \ 0\left( \exists P^{2}|n\right) \ 1\left( n=1\right) \end{cases}$$
直接筛就可以了
#include<cstdio>
const int MAXN = 1e4 + 10;
int N, prime[MAXN], vis[MAXN], mu[MAXN], tot;
void GetMu(int N) {
vis[1] = mu[1] = 1;
for(int i = 2; i <= N; i++) {
if(!vis[i]) prime[++tot] = i, mu[i] = -1;
for(int j = 1; j <= tot && i * prime[j] <= N; j++) {
vis[i * prime[j]] = 1;
if(!(i % prime[j])) {
mu[i * prime[j]] = 0;
break;
}
mu[i * prime[j]] = mu[i] * mu[prime[j]];
//根据莫比乌斯函数的定义，这里也可以写为
//mu[i * prime[j]] = -mu[i];
}
}
}
int main() {
GetMu(1e3);
for(int i = 1; i <= tot; i++)
printf("%d ", mu[i]);
return 0;
}线性筛欧拉函数
这个貌似更常用一点qwq。
我在以前的文章中也详细的讲过，这里只放一下代码
#include<cstdio>
const int MAXN = 1e4 + 10;
int N, prime[MAXN], vis[MAXN], phi[MAXN], tot;
void GetPhi(int N) {
vis[1] = phi[1] = 1;
for(int i = 2; i <= N; i++) {
if(!vis[i]) prime[++tot] = i, phi[i] = i - 1;
for(int j = 1; j <= tot && i * prime[j] <= N; j++) {
vis[i * prime[j]] = 1;
if(!(i % prime[j])) {
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
phi[i * prime[j]] = phi[i] * phi[prime[j]];
}
}
}
int main() {
GetPhi(1e3);
for(int i = 1; i <= tot; i++)
printf("%d ", phi[i]);
return 0;
}线性筛约数个数
这个就稍微高端一点了，按照上面的套路，我们只需要考虑最小的质因子对每个数的贡献
$n = p_1^{a_1} p_2^{a_2} \dots p_k^{a_k}$
设$d(i)$表示$i$的约数个数
那么根据约数定理
$d(i) = \prod_{i = 1}^k (a_i + 1) $
记$a(i)$表示$n$的最小的质因子($a_1$)的指数。
$d(p_i) = 2$
当$i % p_j = 0$时，考虑$i * p_j$，实际上也就是$a_i$的指数多了$1$
我们先除去原来的，再加上新的就OK了
#include<cstdio>
const int MAXN = 1e4 + 10;
int N, prime[MAXN], vis[MAXN], D[MAXN], a[MAXN], tot;
void GetD(int N) {
vis[1] = D[1] = a[1] = 1;
for(int i = 2; i <= N; i++) {
if(!vis[i]) prime[++tot] = i, D[i] = 2, a[i] = 1;
for(int j = 1; j <= tot && i * prime[j] <= N; j++) {
vis[i * prime[j]] = 1;
if(!(i % prime[j])) {
D[i * prime[j]] = D[i] / (a[i] + 1) * (a[i] + 2);
a[i * prime[j]] = a[i] + 1;
break;
}
D[i * prime[j]] = D[i] * D[prime[j]];
a[i * prime[j]] = 1;
}
}
}
int main() {
GetD(1e3);
for(int i = 1; i <= tot; i++)
printf("%d ", D[i]);
return 0;
}线性筛约数和
同样，按照上面的套路考虑
$n = p_1^{a_1} p_2^{a_2} \dots p_k^{a_k}$
设$SD(i)$表示$i$的约数和
$SD(n) = \prod_{i = 1}^k (\sum_{j = 1}^{a_i} p_i^j)$
$sum(i)$表示$i$最小的质因子的贡献，即$sum(i) = \sum_{i = 1}^{a_1}p_1^j$
$low(i)$表示$i$最小质因子的指数，$low(i) = a_1$
有了这三个我们就可以转移了
同样是考虑$i$的最小的因子对答案的贡献，应该比较好推，看代码吧
#include<cstdio>
const int MAXN = 1e4 + 10;
int N, prime[MAXN], vis[MAXN], SD[MAXN], sum[MAXN], low[MAXN], tot;
void GetSumD(int N) {
vis[1] = SD[1] = low[1] = sum[1] = 1;
for(int i = 2; i <= N; i++) {
if(!vis[i]) prime[++tot] = i, sum[i] = SD[i] = i + 1, low[i] = i;
for(int j = 1; j <= tot && i * prime[j] <= N; j++) {
vis[i * prime[j]] = 1;
if(!(i % prime[j])) {
low[i * prime[j]] = low[i] * prime[j];
sum[i * prime[j]] = sum[i] + low[i * prime[j]];
SD[i * prime[j]] = SD[i] / sum[i] * sum[i * prime[j]];
break;
}
low[i * prime[j]] = prime[j];
sum[i * prime[j]] = prime[j] + 1;
//这里low和sum不是积性函数
SD[i * prime[j]] = SD[i] * SD[prime[j]];
}
}
}
int main() {
GetSumD(1e3);
for(int i = 1; i <= tot; i++)
printf("%d ", SD[i]);
return 0;
}非常见积性函数的筛法
很多情况下我们会遇到求两个积性函数狄利克雷卷积的情况
很显然，这个函数也是积性函数，我们考虑如何求得
为了方便筛，我们需要把问题无限简化，
设$low(i)$表示$p_1^{a_1}$
考虑筛法中最关键的地方
$i \% p_j = 0$，
有了$low(i)$，此时我们需要分两种情况讨论
- $low(i) = i$，此时$i$一定是某个素数的幂的形式（否则就会break掉）
这里就用到了我最开始说的那个套路
如果我们能快速的利用$f(p_i^{k})$更新出$f(p_i^{k + 1})$，那这个素数就很容易筛了
- $low(i) \not = i$，那么$i / low(i)$一定与$low(i) * p_j$是互质的，我们可以直接利用积性函数的性质去更新
C++版的伪代码
vis[1] = low[1] = 1; H[1] = 初始化
for(int i = 2; i <= N; i++) {
if(!vis[i]) prime[++tot] = i, mu[i] = -1, H[i] = 质数的情况, low[i] = i;
for(int j = 1; j <= tot && i * prime[j] <= N; j++) {
vis[i * prime[j]] = 1;
if(!(i % prime[j])) {
low[i * prime[j]] = (low[i] * prime[j]);
if(low[i] == i) H[i * prime[j]] = 特殊判断;
else H[i * prime[j]] = H[i / low[i]] * H[prime[j] * low[i]];
break;
}
H[i * prime[j]] = H[i] * H[prime[j]];
low[i * prime[j]] = prime[j];
}
}