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

· · 个人记录

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

质数，是指除了 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) 的欧拉筛。

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

以下以求 \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})，不能接受。考虑使用欧拉筛求解，分类讨论：

于是就可以在 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 \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)。

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

\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 筛，洲阁筛等，将会在后面继续学习探讨。