对质数以及一些筛法的研究

数论是数学对正整数进行研究的分支。筛法最初起源于找出质数的过程中。以下将浅谈数论中各种各样的筛法以及它们的应用。

质数，是指除了 \(1\) 和本身之外不能被任何正整数整除的数，即恰好只有两个因数。于是不难写出以下判断代码：

bool isPrime(int x){
    for(int i=2;i<=x;i++)
        if(x%i==0)
            return false;
    return true;
}

观察上面这段代码，会发现它有很大的冗余。由于因子总是成对出现的，故如果不能被较小的因子整除，那么也一定不能被较大的因子整除。于是只需要枚举不超过 \(\sqrt{x}\) 的数即可。

bool isPrime(int x){
    for(int i=2;i*i<=x;i++)
        if(x%i==0)
            return false;
    return true;
}

但如果判断很多个数是否为质数，或者求一段范围内的质数，这种方法就显得有点无力了。于是就有了筛法。

按照质数的定义，含有除 \(1\) 和自己之外的因子的数不是质数，那么我们不妨用每个数筛掉它的倍数，最后剩下的就是质数。不难写出如下代码，时间复杂度 \(O(n\log n)\)：

for(int i=2;i<=n;i++)
    for(int j=i*2;j<=n;j+=i)
        isPrime[j]=false;

这样做不够快。在这份代码中，一个合数不仅会被质因子筛掉，还会被质因子的倍数筛掉，是没有必要的。于是我们只需枚举质数的倍数即可。得到的是 \(O(n\log \log n)\) 的埃氏筛：

for(int i=2;i<=n;i++)
    if(isPrime[i])
        for(int j=i*2;j<=n;j+=i)
            isPrime[j]=false;

还可以更快吗？一个合数可能有多个不同的质因子，会被多次筛掉，如果它只会被最小的质因子筛掉，就可以有效优化复杂度，于是有了以下代码：

for(int i=2;i<=n;i++){
    if(isPrime[i])Prime[++tot]=i;
    for(int j=1;j<=tot&&i*Prime[j]<=n;j++){
        isPrime[i*Prime[j]]=false;
        if(i%Prime[j]==0)break;//此时说明 prime_j 已经不是 i 的最小质因数，直接退出即可
    }
}

这样做能够保证每个合数只被它的最小质因数筛掉，得到的就是时间复杂度 \(O(n)\) 的欧拉筛。

欧拉筛在求积性函数方面有非常重要的作用。以下补充一些积性函数相关的定义：

• 定义域为正整数集，值域为复数集的子集的函数称为数论函数。
• 满足 \(\forall x\perp y \in \mathbb{N_+},f(xy)=f(x)f(y)\) 的数论函数称为积性函数。
• 满足 \(\forall x,y\in \mathbb{N_+},f(xy)=f(x)f(y)\) 的数论函数称为完全积性函数。

以下以求 \(\sum\limits_{i=1}^n \varphi(i)\) 为例介绍欧拉筛在求积性函数中的应用。\(\varphi(i)\) 定义为 \(1 \sim i\) 的正整数中与 \(i\) 互质的数的个数，即 \(\varphi(i)=\sum\limits_{j=1}^i [\gcd(j,i)=1]\)。显然欧拉函数是积性函数，证明略。

首先考虑单个数 \(n\) 的欧拉函数怎么求，设 \(n\) 的质因数分解的结果为 \(p_1^{k_1}p_2^{k_2}\cdots p_m^{k_m}\)，那么与 \(n\) 互质的条件是不含有任何一个相同的质因子，所以 \(\varphi(n)=n(1-\frac{1}{p_1})(1-\frac{1}{p_m})\cdots (1-\frac{1}{p_m})\)，就可以用一开始的算法 \(O(\sqrt{n})\) 求欧拉函数。

int Phi(int n){
    int phi=n;
    for(int i=2;i*i<=n;i++)
        if(n%i==0){
            phi=phi/i*(i-1);
            while(n%i==0)
                n/=i;
        }
    return phi;
}

那么如果要求 \(1\sim n\) 的欧拉函数呢？一个一个求总复杂度将会达到 \(O(n\sqrt{n})\)，不能接受。考虑使用欧拉筛求解，分类讨论：

• \(i\perp \text{Prime}_j\)，此时根据积性函数的定义，\(\varphi(i\times \text{Prime}_j)=\varphi(i)\times \varphi( \text{Prime}_j)\)。

• 否则，说明 \(i\) 包含质因子 \(\text{Prime}_j\)，那么括号内的部分不变，只是前面的 \(n\) 的部分多乘上了 \(\text{Prime}_j\)，\(\varphi(i\times\text{Prime}_j)=\varphi(i)\times \text{Prime}_j\)。

于是就可以在 \(O(n)\) 的时间复杂度求出 \(\sum\limits_{i=1}^n \varphi(i)\)。

int Euler_Sumphi(int n){
    phi[1]=1;
    for(int i=2;i<=n;i++){
        if(isPrime[i])Prime[++tot]=i,phi[i]=i-1;
        for(int j=1;j<=tot&&i*Prime[j]<=n;j++){
            isPrime[i*Prime[j]]=false;
            if(i%Prime[j]==0){
                phi[i*Prime[j]]=phi[i]*Prime[j];
                break;
            }
            phi[i*Prime[j]]=phi[i]*(Prime[j]-1);
        }
    }
    int sum=0;
    for(int i=1;i<=n;i++)
        sum+=phi[i];
    return sum;
}

但是当 \(n\) 逐渐变大时，这样做也不够快了。有没有更快的办法？有的，这就是杜教筛。以下先普及狄利克雷卷积和莫比乌斯反演。

狄利克雷卷积：对于数论函数 \(f(x),g(x)\)，定义 \(h(x)=\sum\limits_{d|x}f(d)g(\frac{x}{d})\) 为 \(f,g\) 的狄利克雷卷积，记为 \(h=f*g\)。满足交换律，结合律。

\[I(n)=1,\varepsilon(n)=[n=1],Id(n)=n \]

\[\mu(n)=\left\{ \begin{array}{l} 1,n=1,\\ (-1)^k,n=p_1p_2\cdots p_k,p 为不同质数,\\ 0,\text{Otherwise} \end{array} \right. \]

莫比乌斯反演：

\[g(n)=\sum_{d|n}f(d)\Leftrightarrow f(n)=\sum_{d|n}\mu(d)g(\frac{n}{d}) \]

证明：由 \(\mu*I=\varepsilon\)，\(g=f*I\)，得 \(g*\mu=f*I*\mu=f*\varepsilon=f\)，得证。

欧拉函数补充性质：\(n=\sum\limits_{d|n}\varphi(d)\)。证明：设 \(\gcd(n,k)=d,1\le k<n\)，则有 \(\gcd(\frac{n}{d},\frac{k}{d})=1\)，设 \(f_i\) 表示 \(d=i\) 的 \(d\) 的个数，那么显然 \(n=\sum\limits_{i=1}^n f(i)\)，注意到 \(f(x)=\varphi(\frac{n}{x})\)，因此 \(n=\sum\limits_{d|n}\varphi(\frac{n}{x})\)，由对称性得 \(n=\sum\limits_{d|n}\varphi(d)\)，证毕。

接下来推导杜教筛公式。假设要求 \(S(n)=\sum\limits_{i=1}^n f(i)\)。

考虑求狄利克雷卷积前缀和：

\[\begin{align*} \sum_{i=1}^n (f*g)(i)&=\sum_{i=1}^n \sum_{d|i}f(d)g(\frac{n}{d})\\ &=\sum_{d=1}^ng(d)\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} f(i)\\ &=\sum_{d=1}^n g(d)S(\lfloor \frac{n}{d} \rfloor) \end{align*} \]

移项得 \(g(1)S(n)=\sum\limits_{i=1}^n (f*g)(i)-\sum\limits_{d=2}^n g(d)S(\lfloor \frac{n}{d} \rfloor)\)

有了上面这些，就可以开始推式子了。由 \(\varphi*I=Id\)，令 \(f=\varphi,g=I\)，则 \(\sum\limits_{i=1}^n (f*g)(i)=\sum\limits_{i=1}^ni=\frac{n(n+1)}{2}\)，第二部分只需要对 \(S\) 整除分块即可。此时时间复杂度 \(O(n^{\frac{3}{4}})\)。

int Dujiao_Sumphi(int n){
    int sumphi=n*(n+1)/2;
    for(int l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        sumphi-=(r-l+1)*Dujiao_Sumphi(n/l);
    }
    return sumphi;
}

如果预处理 \(1\sim n^{\frac{2}{3}}\) 的 \(S\)，那么时间复杂度可以优化到 \(O(n^{\frac{2}{3}})\)。优化后的代码如下：

int Dujiao_Sumphi(int n){
    int sumphi=n*(n+1)/2;
    if(n<=N)return Pre[N];//N=n^{2/3}
    for(int l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        sumphi-=(r-l+1)*Dujiao_Sumphi(n/l);
    }
    return sumphi;
}

对于其它积性函数的求和，也可以通过构造适当的 \(g\) 从而简化公式右边的计算。当然还有其它不错的解法比如 min25 筛，洲阁筛等，将会在后面继续学习探讨。