再探欧式筛——一种泛用性更强的欧拉筛法/线性筛法实现

JustinRochester · 2023-11-27 20:36:18 · 个人记录

一、引言

欧式筛/欧拉筛法/线性筛法（Euler Sieve）是一种能够在 O(n) 时间复杂度内，处理 [1,n] 内质数的方法。

其相比埃氏筛/埃拉托斯特尼筛法（Eratosthenes Sieve）的 O(n\log\log n) 时间复杂度，主要的优化在于欧式筛保证了所有正整数 n 均只被其最小质因数 {minp}_n 筛到，从而保证了线性复杂度。

此外，欧式筛同时还具有线性时间内预处理欧拉函数 \boldsymbol \varphi(n)、莫比乌斯函数 \boldsymbol \mu(n) 等简单积性函数的能力。但对于除数个数函数 \boldsymbol d(n) （或记为 \boldsymbol {\sigma_0}(n)）、除数和函数 \boldsymbol \sigma(n) 等比较复杂的积性函数，则实现起来比较复杂。

本文提出了一个能处理上述复杂函数，泛用性更强的一种欧拉筛实现；并通过实验证明，这个实现的效率也非常高。

二、相关基础内容

2.1 筛法

我们称，一个只有 1 和其本身是自己因数的正整数 n 为质数。一个比较简单的求解 [1,n] 范围内质数的方式为：

对已知的质数，划去所有其倍数。未被划去的数即为质数。参考代码如下：

bool nprime[MAXN];
for(int i=2; i<=n; ++i) if(!nprime[i])
    for(int j=i+i; j<=n; j+=i)
        nprime[j]=1;

此即为埃氏筛法，时间复杂度可证明为 O(n\log \log n)。

为保证时间复杂度为线性，我们需要保证每个数 n 仅被自身的最小质因数 {minp}_n 筛到。当然，这个描述的实现方法就和动态规划中的填表法一样，不容易理解；我们也可以参考动态规划的思想，转化为刷表法实现。参考代码如下：

vector<int> pr;
bool nprime[MAXN];
for(int i=2; i<=n; ++i) {
    if(!nprime[i])
        pr.push_back(i);
    for(int p : pr)
        if(p*i>n)
            break;
        else {
            nprime[p*i]=1;
            if(i%p==0) //再大的就不是枚举最小质因数了
                break;
        }
}

于是，我们就实现了 O(n) 的欧式筛法。

2.2 积性函数和狄利克雷卷积

我们称，对任意满足 \gcd(a,b)=1 的 a,b，都满足 \boldsymbol f(ab)=\boldsymbol f(a)\boldsymbol f(b) 的函数为积性函数，用粗体表示。

很显然，根据积性函数的性质可以有：

其中 p_i 是互不相同的质数。因此，一个比较简单的描述积性函数的方法就是，直接描述它在质数幂次的结果

通过上式可以很简单地推出，积性函数 \boldsymbol f(n) 和积性函数 \boldsymbol g(n) 的点积 \boldsymbol f(n)\cdot \boldsymbol g(n) 也是积性函数。为此，我们在后文将直接使用 (\boldsymbol f\cdot \boldsymbol g)(n) 来表示这种点积后的积性函数。

此外，我们还能定义狄利克雷卷积运算，称积性函数 \boldsymbol f(n) 和积性函数 \boldsymbol g(n) 的狄利克雷卷积为：

通过定义式也可以证明，狄利克雷卷积也是积性函数。

2.3 简单积性函数的线性筛法

我们定义积性函数：欧拉函数 \boldsymbol \varphi(n) 表示 [1,n] 中与 n 互质的正整数个数。

此篇文章 中证明了，欧拉函数在质数幂次 p^e(e>0) 处的结果为：\boldsymbol \varphi(p^e)=p^{e-1}(p-1)。

此外，另一个比较常见的积性函数是莫比乌斯函数 \boldsymbol \mu(n)，其在质数幂次 p^e(e>0) 处的结果为：

这两个函数的一个共性是：其在质数幂次处的结果可以很容易地整体乘法递推。我们假设质数 p<{minp}_n,e>0，则从已知 \boldsymbol f(p^e\cdot n) 递推出 \boldsymbol f(p^{e+1}\cdot n) 是很容易的：

因此，我们可以拓展一下欧式筛，在刷表的时候，确认一下当前枚举的质数是不是自己的最小质因数，即可根据积性性质或上述递推性质分别转移了。

vector<int> pr;
int phi[MAXN], mu[MAXN];
bool nprime[MAXN];
for(int i=2; i<=n; ++i) {
    phi[1]=1;
    mu[1]=1;
    if(!nprime[i]) {
        pr.push_back(i);
        phi[i]=i-1;
        mu[i]=-1;
    }
    for(int p : pr)
        if(p*i>n)
            break;
        else {
            nprime[p*i]=1;
            if(i%p==0) {//递推转移
                phi[p*i]=p*phi[i];
                mu[p*i]=0;
                break;
            }
            else {//积性转移
                phi[p*i]=phi[p]*phi[i];
                mu[p*i]=mu[p]*mu[i];
            }
        }
}

这个方法的使用很广泛，然而当积性函数不满足可乘法递推性质时，这个实现会使得递推转移部分的代码写得非常复杂；甚至若维护不好，可能还会在维护过程中失去线性时间复杂度。

例如，对于除数个数函数，其在质数幂次 p^e(e>0) 处的结果为 \boldsymbol {\sigma_0}(p^e)=e+1。若强行维护，则 \boldsymbol {\sigma_0}(n\cdot p^{e+1})={e+2\over e+1}\cdot \boldsymbol {\sigma_0}(n\cdot p^e)。若想到了维护最小质因数的幂次数组，还尚且能 O(n) 转移，否则还要在每次花费 O(e) 的复杂度处理 {e+2\over e+1}，时间复杂度为：

其中 \text{Prime} 表示 [1,n] 中的质因数集合。

很显然，这个实现已经不够理想了，若是更复杂的积性函数 \displaystyle \boldsymbol {\sigma_k}(p^e)=\sum_{i=0}^e p^{ki}，维护起来只会更加复杂。

三、再探欧拉筛

3.1 优化转移

上述欧拉筛造成积性函数只方便递推转移的关键在于，我们在递推过程中只能方便地表示每个数的最小质因数。而在大部分的积性函数中，我们要知道 \boldsymbol f(p^e\cdot n)(p<{minp}_n,e>0)，是需要知道最小质因数的幂次 p^e。

这就启示我们能不能额外维护最小质因数的幂次 {pe}_n，这样我们即可通过查表快速计算积性函数 \boldsymbol f(n)=\boldsymbol f({pe}_n)\cdot \boldsymbol f({n\over {pe_n}})。

事实上，{pe}_n 是很好处理的。我们在枚举不超过 {minp}_n 的质数 p 时，若 p<{minp}_n，则 {pe}_{n\cdot p}=p，否则 {pe}_{n\cdot p}={pe}_n\cdot p；这在欧式筛中是很容易实现的。

并且，由于维护了 {pe}_n，我们可以不必再维护 {nprime}_n 来判断一个数是否不是质数。参考代码如下：

vector<int> pr;
//...要筛的的积性函数
int pe[MAXN];
//f[1]=1;
for(int i=2, x; i<=n; ++i) {
    if(!pe[i]) {
        pr.push_back(i);
        pe[i]=i;
        //...要筛的积性函数
    }
    for(int p : pr)
        if((x=p*i)>n)//记忆 x=p*i，避免后续重复乘法
            break;
        else {
            if(i%p==0) {
                pe[x]=pe[i]*p;
                //f[x]=f[x]*f[x/pe[x]]
                break;
            }
            else {//积性转移
                pe[x]=p;
                //f[x]=f[i]*f[p];
            }
        }
}

3.2 质数部分的处理

若仔细思考上述代码，很快就会发现一个问题：

虽然 \boldsymbol f(n)=\boldsymbol f({pe}_n)\cdot \boldsymbol f({n\over {pe_n}}) 在大部分情况下都是对的，但按照该写法，n=p^e 就会有：

但是，这个转移式并不能求出在质数幂次时的积性函数值！

因此，关于质数幂次处的积性函数值，必须在筛到质数的时候同步处理：当我们枚举到质数 p 时，我们枚举 [1,n] 范围内所有 p 的幂次 p^e，根据定义计算其结果。

一个比较朴素的实现方法是：通过不停让一个初始为 p 的变量乘以 p，若超过 n 了则直接退出。

但是这个实现方法是非常危险的，当我们筛到的质数 p 接近 n 时，我们进行一次乘法则会直接达到 n^2 的规模，很容易发生溢出。

除了开 long long/unsigned long long 以外，还有一个比较巧妙但常数较小的方法如下：

我们只要在维护 p^e 的同时再维护 p^{e-1} ，通过 p^{e-1} 是否大于 \lfloor{n\over p}\rfloor 即可；而由于我们又需要在后续判断 i\cdot p>n ，这也可以同样等价于 i>\lfloor{n\over p}\rfloor，这一次除法很好地解决了溢出问题和后续多次乘法的问题。

3.3 处理质数幂次的时间复杂度

很显然，若积性函数满足原来的可乘法递推性时，它在 \boldsymbol f(n) 转移到 \boldsymbol f(n\cdot p) 的时间复杂度是 O(1) 的，因此它从 \boldsymbol f(p^e) 转移到 \boldsymbol f(p^{e+1}) 显然也是 O(1) 的。因此这一部分的时间复杂度和所有质因数的幂次数显然是等价的，而该数显然不超过 n，因此时间复杂度仍然为 O(n)。

我们不妨考虑不满足可乘法递推性时，均摊到质数幂次的时间复杂度容许为多高？不妨假设所有质数幂次的均摊复杂度为 T_2(n)，其中筛到各质数时需要先花费均摊为 T_1(n) 的时间复杂度预处理，则：

根据质数定理，\pi(n)\sim O({n\over \ln n})；为保证总复杂度为 O(n)，故 T_1(n)=O(\log n)；而 T_2(n) 则比较繁琐：

因此，同理可得，T_2(n)=O(\log n) 即可保证线性复杂度

3.4 一些细节优化

当然，我们知道，维护 n 中最小质因数 {pe}_n 和维护 n 中除了最小质因数的部分 {fr}_n 是等价的。那两个到底谁的效率更高呢？

我们分别考虑枚举的 p<{minp}_n 和 p={minp}_n 的情况：

当 p<{minp}_n 时，

\begin{cases} {pe}_{n\cdot p}&={pe}n\cdot p\ {fr}{n\cdot p}&={fr}_n\ \end{cases}

其中 (n/d) 表示 \lfloor{n\over d}\rfloor

按照乘法递推的方法即可筛欧拉函数。一个实现方法是：在递推函数中，维护静态变量 lst，表示上一次返回的结果；通过传入参数 p,e 可以快速判断结果的变化，例如：

static uint lst;//unsigned
if(e==1) return lst=p-1;
else return lst=lst*p;

该实现方法的评测 记录 为 31ms；而除本文外，目前最快的 记录 为 24ms。本文耗时仅为其约 1.375 倍。

4.1.2 除数和函数

LOJ #124. 除数函数求和 1

考虑积性函数在质数幂次位置的结果：

因此，可以同上维护静态变量 lst 表示上次的输出结果，也是容易转移的。

static uint lst;
if(e==1) return lst=fpow(p, k)-1;
else return lst=lst*fpow(p, k);

但实际上，由于每次反复计算 fpow(p, k) 会造成额外的 O(\log k) 开销，这里可以再额外引入静态变量 pk 进行维护。

static uint pk = -1, lst = -1;
if (e == 1)
    pk = fpow(p, k), lst = 1 + pk;
else
    lst = ((ull)lst * pk + 1) % P;
return lst;

该实现方法的评测 记录 为 1819ms；而除本文外，目前最快的 记录 为 1646ms。本文耗时约为其 1.105 倍。

4.2 较为复杂的积性函数

4.2.1 \boldsymbol \varphi*\boldsymbol \mu 的欧拉筛

BZOJ #4804. 欧拉心算 法一

考虑积性函数在质数幂次位置的结果：

故 e=1 时 (\boldsymbol\varphi*\boldsymbol \mu)(p^e)=p-1-1=p-2；