浅谈积性函数的线性筛法
前置知识
线性筛
线性筛可以在严格$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)$,此时我们需要分两种情况讨论
1. $low(i) = i$,此时$i$一定是某个素数的幂的形式(否则就会break掉)
这里就用到了我最开始说的那个套路
如果我们能快速的利用$f(p_i^{k})$更新出$f(p_i^{k + 1})$,那这个素数就很容易筛了
2. $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]; } }
参考资料