筛法
前言
本文介绍了下述筛法：
- 素数筛法。
- 线性筛求积性函数。
- 杜教筛。
- Powerful Number 筛。
素数筛法
在很多时候，我们想知道对于 \(\forall 1\le x\le n\)（\(n\) 为给定的正整数），有哪些 \(x\in\mathbb P\)。
一个自然的想法是对于每个 \(x\)，都进行一次素性检验。
显然，这一暴力的做法无法通过较大的数据范围。
于是，素数筛法就出现了。
埃拉托斯特尼筛（埃筛）
我们发现，对于一个正整数 \(x>1\)，它除本身之外的倍数必定是合数。
利用这一结论，我们可以避免很多不必要的检测。
若我们从小到大考虑每个数 \(x\)，同时将它除本身之外的倍数 \(y\) 都标记为合数，那么最后没被标记的就是质数了。
显然，我们只要考虑所有的 \(2\le x\le\lfloor\sqrt{n}\rfloor\) 即可。
实现
int P[N],cnt;
bool np[N];
void Eratosthenes(int n){
cnt=0;
np[0]=np[1]=true;
for(int x=2;x*x<=n;++x)if(!np[x])for(int y=x*x;y<=n;y+=x)np[y]=true;
for(int x=2;x<=n;++x)if(!np[x])P[++cnt]=x;
}
可以证明，埃筛的时间复杂度为 \(O(n\log\log n)\)。
一些优化
只筛奇数
由于除了 \(2\) 的偶数都是合数，因此我们只需考虑奇数即可。
这不仅能让我们的内存需求减半，还能让所需的操作减半。
压位
我们使用的 np 数组是 bool 类型的。一般一个 bool 类型变量占 \(1\) 字节（即 \(8\) 比特），但储存一个 bool 值只需 \(1\) 比特。
利用位运算相关知识，我们可以将每个 bool 值压到一个比特位里，这样我们仅需 \(n\) 比特（即 \(\displaystyle\frac{n}{8}\) 字节）而并非 \(n\) 字节，可以显著减少内存占用。这种方式被称为压位。
在 STL 中，vector<bool> 和 bitset<> 都会进行压位，并且有常数优化（后者稍微好一点）。
若觉得常数还是过大，可以手写 bitset。不过这就不是本篇文章要讨论的东西了。
线性筛（欧拉筛）
注意到在埃筛中会存在一个合数被多次标记的情况。
若能使每个合数只被标记一遍，就可以把时间复杂度改进为 \(O(n)\)。
int P[N],cnt;
bool np[N];
void Euler(int n){
cnt=0;
np[0]=np[1]=true;
for(int x=2;x<=n;++x){
if(!np[x])P[++cnt]=x;
for(int i=1;i<=cnt&&P[i]*x<=n;++i){
np[P[i]*x]=true;
if(x%P[i]==0)break;
}
}
}
解释一下为什么每个合数只会被标记一次：
当 \(x\bmod P_i=0\) 时，换言之，\(x\) 之前被 \(P_i\) 筛过了。
由于 \(P\) 保存的质数是从小到大的，因此 \(x\) 乘上后续质数的结果在后面一定会被 \(P_i\) 的倍数筛掉。
因此，这里不需要先筛一次，直接break即可。
线性筛求积性函数
在一类问题中，我们要预处理对于 \(\forall1\le x\le n\)（\(n\) 为给定的正整数），一些积性函数 \(f(x)\) 的值。
如：欧拉函数 \(\varphi(x)\)，莫比乌斯函数 \(\mu(x)\)，除数函数 \(\sigma_k(x)\)。
接下来将一一介绍如何用线性筛求出它们。
在以下的讲述中，皆有 \(\displaystyle x=\prod_{i=1}^qp_i^{\alpha_i}\)，\(p_1\) 为 \(x\) 最小的素因子，\(x'=\dfrac{x}{p_1}\)。
当提到 \(x\not\in\mathbb P\) 时，默认 \(x\neq1\)。
欧拉函数 \(\varphi(x)\)
首先，有 \(\varphi(1)=1\)。
注意到在线性筛中，每一个合数 \(x\) 都是被 \(p_1\) 筛掉的。
显然，在线性筛的过程中 \(x\) 是通过 \(p_1\times x'\) 筛掉的。
下面对 \(x'\bmod p_1\) 进行分类讨论：
- 若 \(x'\bmod p_1=0\)：
那么 \(x'\) 包含了 \(x\) 的所有的质因子。
此时有
- 若 \(x'\bmod p_1\neq0\)：
那么 \(x'\) 和 \(p_1\) 是互质的，所以有
当 \(x\in\mathbb P\) 时，则有 \(\varphi(x)=x-1\)。
实现
int P[N],cnt;
bool np[N];
int phi[N];
void CalcPhi(int n){
cnt=0;
np[0]=np[1]=true;
phi[1]=1;
for(int x_=2;x_<=n;++x_){
if(!np[x_]){
P[++cnt]=x_;
phi[x_]=x_-1;
}
for(int i=1;i<=cnt&&P[i]*x_<=n;++i){
int x=P[i]*x_;
np[x]=true;
if(x_%P[i]==0){
phi[x]=P[i]*phi[x_];
break;
}
phi[x]=phi[P[i]]*phi[x_];
}
}
}
莫比乌斯函数 \(\mu(x)\)
当 \(x=1\)：\(\mu(x)=1\)。
当 \(x\in\mathbb P\) 时：\(\mu(x)=-1\)。
当 \(x\not\in\mathbb P\) 时：
- \(x'\bmod p_1=0\)：
此时有 \(x\bmod p_1^2=0\)，所以 \(\mu(x)=0\)。 - \(x'\bmod p_1\neq0\)：
此时 \(\mu(x)=\mu(p_1)\times\mu(x')=-\mu(x')\)。
实现
int P[N],cnt;
bool np[N];
int mu[N];
void CalcMu(int n){
cnt=0;
np[0]=np[1]=true;
mu[1]=1;
for(int x_=2;x_<=n;++x_){
if(!np[x_]){
P[++cnt]=x_;
mu[x_]=-1;
}
for(int i=1;i<=cnt&&P[i]*x_<=n;++i){
int x=P[i]*x_;
np[x]=true;
if(x_%P[i]==0)break;
mu[x]=-mu[x_];
}
}
}
除数函数 \(\sigma_k(x)\)
令 \(f(x)=p_1^{\alpha_1}\)。
当 \(x=1\)：\(\sigma_k(x)=1\)。
当 \(x\in\mathbb P\)：\(f(x)=x,\sigma_k(x)=x^k+1\)。
当 \(x\not\in\mathbb P\)：
- \(x'\bmod p_1=0\)：
\(f(x)=p_1\times f(x'),\sigma_k(x)=\begin{cases}\sigma_k(f(x))\times\sigma_k(\dfrac{x}{f(x)}),&f(x)\neq x\\ p_1^k\times\sigma_k(x') +1,&f(x)=x\end{cases}\)。 - \(x'\bmod p_1\neq0\)：
\(f(x)=p_1,\sigma_k(x)=\sigma_k(x')\times\sigma_k(p_1)\)。
实现
int QuickPow(int x,int y){
int ans=1;
while(y){
if(y&1)ans=ans*x;
x=x*x;
y>>=1;
}
return ans;
}
int P[N],cnt;
bool np[N];
int f[N],Pow[N],sigma[N];
void CalcSigma(int n,int k){
cnt=0;
np[0]=np[1]=true;
sigma[1]=1;
for(int x_=2;x_<=n;++x_){
if(!np[x_]){
++cnt;
P[cnt]=x_;
Pow[cnt]=QuickPow(x_,k);
f[x_]=x_;
sigma[x_]=Pow[cnt]+1;
}
for(int i=1;i<=cnt&&P[i]*x_<=n;++i){
int x=P[i]*x_;
np[x]=true;
if(x_%P[i]==0){
f[x]=P[i]*f[x_];
if(f[x]!=x)sigma[x]=sigma[f[x]]*sigma[x/f[x]];
else sigma[x]=Pow[i]*sigma[x_]+1;
break;
}
f[x]=P[i];
sigma[x]=sigma[P[i]]*sigma[x_];
}
}
}
一般的积性函数 \(f(x)\)
和求 \(\sigma_k(x)\) 类似，现在将过程概括如下。
令 \(g(x)=p_1^{\alpha_1}\)。
当 \(x=1\)：显然有 \(f(x)=1\)。
当 \(x\in\mathbb P\)：\(g(x)=x\)，\(f(x)\) 根据函数的性质 \(O(1)\) 计算。
当 \(x\not\in\mathbb P\)：
- \(x'\bmod p_1=0\)：
\(g(x)=g(x')\times p_1\)。
若 \(g(x)=x\)，则说明 \(x=p_1^{\alpha_1}\)，可以 \(O(1)\) 计算 \(f(x)\)；否则 \(f(x)=f(g(x))\times f\left(\dfrac{x}{g(x)}\right)\)。 - \(x'\bmod p_1\neq0\)：
\(g(x)=p_1\)，\(f(x)=f(p_1)\times f(x')\)。
杜教筛
杜教筛是一种利用狄利克雷卷积和数论分块，在低于线性时间的复杂度内，计算数论函数 \(f(n)\) 前缀和 \(S(n)\) 的方法。
前置知识
欧拉函数 \(\varphi(n)\)
考虑到 \(\operatorname{id}=\varphi*1\)，则有：
因此 \(\displaystyle S(n)=\dfrac{n(n+1)}{2}-\sum_{d=2}^xS(\lfloor\dfrac{n}{d}\rfloor)\)
对于较小的 \(n\)，我们可以使用线性筛预处理出来；
若 \(n\) 的最大值为 \(N\)，则可以筛出 \(N^{\frac{2}{3}}\) 内的值，此时时间复杂度为 \(O(n^{\frac{2}{3}})\)。
为了进一步加快计算，可以使用记忆化：将计算过的值用 map/unordered_map 存下来。
实现
int P[N],cnt;
bool np[N];
int phi[N],SumOfPhi[N];
void EulerSieveSumOfPhi(int n){
cnt=0;
np[0]=np[1]=true;
phi[1]=1;
for(int x_=2;x_<=n;++x_){
if(!np[x_]){
P[++cnt]=x_;
phi[x_]=x_-1;
}
for(int i=1;i<=cnt&&P[i]*x_<=n;++i){
int x=P[i]*x_;
np[x]=true;
if(x_%P[i]==0){
phi[x]=P[i]*phi[x_];
break;
}
phi[x]=phi[P[i]]*phi[x_];
}
}
for(int x=1;x<=n;++x)SumOfPhi[x]=SumOfPhi[x-1]+phi[x];
}
unordered_map<int,int> MapSumOfPhi;
int CalcSumOfPhi(int n){
if(n<=MAXN)return SumOfPhi[n];
if(MapSumOfPhi[n])return MapSumOfPhi[n];
int ans=n*(n+1)/2;
for(int l=2,r;l<=n;l=r+1){
r=n/(n/l);
ans-=(r-l+1)*CalcSumOfPhi(n/l);
}
return MapSumOfPhi[n]=ans;
}
莫比乌斯函数 \(\mu(n)\)
原理与欧拉函数 \(\varphi(n)\) 类似，只不过利用的是 \(\epsilon=\mu*1\)：
因此 \(S(n)=1-\sum_{d=2}^nS(\lfloor\dfrac{n}{d}\rfloor)\)。
实现
int P[N],cnt;
bool np[N];
int mu[N],SumOfMu[N];
void EulerSieveSumOfMu(int n){
cnt=0;
np[0]=np[1]=true;
mu[1]=1;
for(int x_=2;x_<=n;++x_){
if(!np[x_]){
P[++cnt]=x_;
mu[x_]=-1;
}
for(int i=1;i<=cnt&&P[i]*x_<=n;++i){
int x=P[i]*x_;
np[x]=true;
if(x_%P[i]==0)break;
mu[x]=-mu[x_];
}
}
for(int x=1;x<=n;++x)SumOfMu[x]=SumOfMu[x-1]+mu[x];
}
unordered_map<int,int> MapSumOfMu;
int CalcSumOfMu(int n){
if(n<=MAXN)return SumOfMu[n];
if(MapSumOfMu[n])return MapSumOfMu[n];
int ans=1;
for(int l=2,r;l<=n;l=r+1){
r=n/(n/l);
ans-=(r-l+1)*CalcSumOfMu(n/l);
}
return MapSumOfMu[n]=ans;
}
一般的数论函数 \(f(n)\)
找到一个合适的数论函数 \(g(n)\)，令 \(h=f*g\)，则有：
因此 \(\displaystyle S(n)=\dfrac{\sum_{i=1}^nh(i)-\sum_{d=2}^ng(d)S(\lfloor\dfrac{n}{d}\rfloor)}{g(1)}\)。
这就要求：
- 能快速计算 \(h(n)\) 的前缀和；
- 能快速计算 \(g(n)\) 的前缀和，以用数论分块求 \(\displaystyle\sum_{d=2}^ng(d)S(\lfloor\dfrac{n}{d}\rfloor)\)。
Powerful Number 筛
Powerful Number 筛是一种利用 Powerful Number（以下简称 PN）的性质，在低于线性的时间复杂度下，计算一些积性函数 \(f(n)\) 的前缀和 \(F(n)=\displaystyle\sum_{i=1}^nf(i)\) 的筛法。
PN 筛对积性函数的性质要求较为严苛，要求存在一个积性函数 \(g(n)\)，使得 \(g(n)\) 易求前缀和 \(G(n)=\displaystyle\sum_{i=1}^ng(i)\)，并且对于质数 \(p\)，都有 \(g(p)=f(p)\)。
什么是 PN
对于一正整数 \(\displaystyle n=\prod_{i=1}^kp_i^{\alpha_i}\)，若对于 \(\forall 1\le i\le k\)，都有 \(\alpha_i>1\)，则称 \(n\) 是 PN。（为了方便表述，下文中我们令 \(1\) 也是 PN）
下面是 PN 两个重要的性质：
- 性质 \(1\)：所有 PN 都可以表示成 \(a^2b^3\) 的形式。
证明：
若 \(2\mid \alpha_i\)，则 \(p_i^{\alpha_i}=\left(p_i^{\frac{\alpha_i}{2}}\right)^2\times1^3\)；否则 \(p_i^{\alpha_i}=\left(p_i^{\frac{\alpha_i-3}{2}}\right)^2\times p_i^3\)。
故 PN 可以表示成 \(a^2b^3\) 的形式。
- 性质 \(2\)：\(n\) 以内的 PN 个数至多有 \(O(\sqrt{n})\) 个。
证明：
考虑枚举 \(a\)，再计算符合条件的 \(b\) 的个数，则 \(n\) 以内的 PN 个数约为\[\int_1^{\sqrt{n}}\sqrt[3]{\dfrac{n}{x^2}}\mathrm{d}x=3\sqrt{n}-3\sqrt[3]{n}=O(\sqrt{n}) \]
其中，性质 \(2\) 是 PN 筛复杂度的保障。
PN 筛的过程
首先找到一个积性函数 \(g(n)\)，使得 \(\forall p\in \mathbb{P},g(p)=f(p)\)，并令 \(h(n)\) 为满足 \(g*h=f\) 的函数。
根据积性函数的性质，\(h(n)\) 也为积性函数；
由 \(f(p)=g(p)h(1)+g(1)h(p)\) 可得，\(h(p)=0\)。
因此，\(h(n)\) 只有当 \(n\) 为 PN 时取到非 \(0\) 值。
然后便有以下式子：
所以，我们可以用 DFS 在 \(O(\sqrt{n})\) 内找出所有 PN，并计算出 \(h\) 的所有有效值；
对于 \(h\) 有效值的计算，只需计算出所有 \(h(p^c)\) 的值，即可通过积性函数的性质计算出 \(h\) 的所有有效值；
对于每个有效值 \(d\)，计算 \(h(d)G\left(\left\lfloor\dfrac{n}{d}\right\rfloor\right)\) 并累加即可得到 \(F(n)\)。
至于如何计算 \(h(p^c)\)，也有两种方法：
第一种是直接推出 \(h(p^c)\) 与 \(p\) 和 \(c\) 有关的公式，然后根据公式计算 \(h(p^c)\)；
第二种是根据 \(f=g*h\)，即 \(\displaystyle f(p^c)=\sum_{i=0}^cg(p^i)h(p^{c-i})\)，可得 \(\displaystyle h(p^c)=f(p^c)-\sum_{i=1}^cg(p^i)h(p^{c-i})\)，现在就可以通过枚举 \(p\) 和 \(c\) 来计算 \(h(p^c)\)。
现在整理下 PN 筛的过程：
- 构造合适的 \(g\)；
- 构造快速计算 \(G\) 的方法；
- 计算 \(h(p^c)\)；
- 搜索 PN，过程中累加答案；
对于计算 \(h(p^c)\) 可以直接根据公式计算，可以预处理，也可以搜到了再临时计算。
PN 筛的复杂度分析
以上文第二种计算 \(h(p^c)\) 的方式为例进行分析，可将过程分为两部分：计算 \(h(p^c)\) 和搜索 PN。
对于第一部分，根据 \(O(\sqrt{n})\) 的素数个数为 \(O\left(\dfrac{\sqrt{n}}{\log n}\right)\)，而指数 \(c\) 至多为 \(\log n\)，而计算 \(h(p^c)\) 需要循环 \(O(c)\) 次。
所以第一部分的时间复杂度为 \(O\left(\dfrac{\sqrt{n}}{\log n}\cdot\log n\cdot\log n\right)=O(\sqrt{n}\log n)\)，这是一个较为宽松的上界。
根据题目的不同还可以加相应的剪枝，进一步减小第一部分的时间复杂度。
对于第二部分，由于 \(n\) 以内的 PN 至多有 \(O(\sqrt{n})\) 个，所以至多搜索 \(O(\sqrt{n})\) 次。
若计算 \(G\left(\left\lfloor\dfrac{n}{d}\right\rfloor\right)\) 的时间复杂度为 \(O(T)\)，则第二部分的时间复杂度为 \(O(T\sqrt{n})\)。
特别地，若用杜教筛计算 \(G\left(\left\lfloor\dfrac{n}{d}\right\rfloor\right)\)，则第二部分的时间复杂度为杜教筛的时间复杂度 \(O(n^{\frac{2}{3}})\)；
其中，杜教筛需要线性筛预处理并记忆化，这样就能保证过程中计算到的 \(G\) 值在计算 \(G(n)\) 时已经计算过。
对于空间复杂度，其瓶颈在于存储 \(h(p^c)\)；
若使用二维数组 \(h\) 记录，\(h_{i,j}\) 表示 \(h(p_i^j)\) 的值，则空间复杂度为 \(O\left(\dfrac{\sqrt{n}}{\log n}\cdot\log n\right)=O(\sqrt{n})\)。
PN 筛的实现
int n,F,h[N][M];bool vish[N][M];
void init(){
F=0;
sievePrimeAndG();
for(int i=1;i<=cnt;++i){
h[i][0]=1,vish[i][0]=true;
h[i][1]=0,vish[i][1]=true;
}
}
void dfs(int pid,int num,int H){
F+=H*G(n/num);
for(int id=pid;id<=cnt;++id){
int p=P[id];
if(num>n/p/p)break;
for(int tnum=num*p*p,pc=p*p,c=2;tnum<=n;tnum*=p,pc*=c,++c){
if(!vish[id][c]){
int res=f(pc);
for(int i=1;i<=c;++i)res-=g(p,c)*h[id][c-i];
h[id][c]=res;
vish[id][c]=true;
}
dfs(id+1,tnum,H*h[id][c]);
}
}
}
int calcF(){
init();
dfs(1,1,1);
return F;
}
//sievePrimeAndG() 是线性筛，顾名思义，是筛质数和 g 的前缀和的。
//P[] 和 cnt 的含义，详见上文介绍线性筛的部分。
//在具体的问题中，f(pc) 可以直接计算，g(p,c) 可以递推计算。
浙公网安备 33010602011771号