筛法

LXcjh4998 · 2025-04-14 21:07:42 · 算法·理论

本文原载于作者的博客园。

前言

本文介绍了下述筛法：

素数筛法。
线性筛求积性函数。
杜教筛。
素数筛法

在很多时候，我们想知道对于 \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 的所有的质因子。
此时有

\varphi(x)&=x\times\prod_{i=1}^q\frac{p_i-1}{p_i}\\ &=p_1\times x'\times\prod_{i=1}^q\frac{p_i-1}{p_i}\\ &=p_1\times\varphi(x') \end{aligned}

若 x'\bmod p_1\neq0：
那么 x' 和 p_1 是互质的，所以有 \varphi(x)=\varphi(p_1)\times\varphi(x')

当 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^2=0$，所以 $\mu(x)=0$。
此时 $\mu(x)=\mu(p_1)\times\mu(x')=-\mu(x')$。 ### 实现 ```cpp 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)$。
### 实现
```cpp
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：

$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)$。
$g(x)=p_1$，$f(x)=f(p_1)\times f(x')$。 # 杜教筛杜教筛是一种利用狄利克雷卷积和数论分块，在低于线性时间的复杂度内，计算数论函数 $f(n)$ 前缀和 $S(n)$ 的方法。 ## 前置知识
数论分块；
狄利克雷卷积。
欧拉函数 \varphi(n)

考虑到 id=\varphi*1，则有：
\dfrac{n(n+1)}{2} &=\sum_{i=1}^ni\\ &=\sum_{i=1}^n\sum_{d\mid i}\varphi(\dfrac{i}{d})\\ &=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\varphi(i)\\ &=\sum_{d=1}^nS(\lfloor\dfrac{n}{d}\rfloor)\\ &=S(n)+\sum_{d=2}^nS(\lfloor\dfrac{n}{d}\rfloor) \end{aligned}
因此 \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：

1 &=\sum_{i=1}^n\epsilon(i)\\ &=\sum_{i=1}^n\sum_{d\mid i}\mu(\dfrac{n}{d})\\ &=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\mu(i)\\ &=\sum_{d=1}^nS(\lfloor\dfrac{n}{d}\rfloor)\\ &=S(n)+\sum_{d=2}^nS(\lfloor\dfrac{n}{d}\rfloor) \end{aligned}

因此 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，则有：

\sum_{i=1}^nh(i) &=\sum_{i=1}^n\sum_{d\mid i}f(\dfrac{n}{d})g(d)\\ &=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}g(d)f(i)\\ &=\sum_{d=1}^ng(d)S(\lfloor\dfrac{n}{d}\rfloor)\\ &=g(1)S(n)+\sum_{d=2}^ng(d)S(\lfloor\dfrac{n}{d}\rfloor) \end{aligned}

因此 \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 值。

然后便有以下式子：

F(n) &=\sum_{i=1}^nf(i)\\ &=\sum_{i=1}^n(g*h)(i)\\ &=\sum_{i=1}^n\sum_{d\mid i}h(d)g\left(\dfrac{i}{d}\right)\\ &=\sum_{d=1}^nh(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}g(i)\\ &=\sum_{d=1}^nh(d)G\left(\left\lfloor\dfrac{n}{d}\right\rfloor\right)\\ &=\sum_{\substack{d\le n\\d \text{ is PN}}}h(d)G\left(\left\lfloor\dfrac{n}{d}\right\rfloor\right) \end{aligned}

所以，我们可以用 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) 可以递推计算。

筛法

前言

素数筛法

埃拉托斯特尼筛（埃筛）

实现

一些优化

只筛奇数

压位

线性筛（欧拉筛）

线性筛求积性函数

欧拉函数 \varphi(x)

实现