浅谈素数筛优化

前置知识：
质数筛、质数判定

之前我没有找到亚线性筛时没有在这篇文章刊登在洛谷日报之前修改，非常抱歉。

现在发现了原来已有的亚线性筛，分别为Sieve of Atkin与Wheel factorization。而这篇文章原来讲的就是Wheel factorization。

下面是这些筛法在wiki上的介绍： https://en.wikipedia.org/wiki/Sieve_of_Atkin
https://en.wikipedia.org/wiki/Wheel_factorization

再下面是一个什么都不会的人自己YY出Wheel factorization的过程，可以不用看了。

大家应该都记得这样一种质数判定的方法：

bool is_prime(int a)
{
    if(a==1)return 0;
    if(a==2||a==3)return 1;
    if(a%6!=1&&a%6!=5)return 0;
    for(int i=5;i<=sqrt(a);i+=6)
        if(a%i==0||a%(i+2)==0)return 0;
    return 1;
}

虽然这种方法貌似没有名字，但这种方法常数极小，带常数的规模大概是\frac{\sqrt n}{3}，是一般素数判断方法效率的3倍，其原理就在于这个方法应用了一个原理：

除2、3外的任何质数，除以6的余数一定是1或5。

证明也很简单：因为除以6余2、4、0的数可以被2整除，除以6余3、0的数可以被3整除。
仅这一简单的技巧，程序效率瞬间快了2倍，质数密度为原来3倍。
而我们既然能用这个技巧判定质数，能不能把它用到筛质数里面呢？
我们知道，平常我们用到的最快的质数筛为欧拉筛，时间复杂度为O(n)，相比之下，慢一点但简单一点的埃氏筛时间复杂度为O(nloglogn)。想要突破O(n)这个时间复杂度，貌似只有在你可以预先知道所有质数的位置时，时间复杂度约\Theta(\frac{n}{log(n)})（即1~n内大概的质数个数）。任何素数筛法都达不到这个

可能大家不是太懂素数公式的有关内容，一篇百度百科介绍了素数及伪素数的通项公式，这里简单介绍一下：

素数及伪素数的通项公式为p(n)=3n+1+sin^2(\frac{n \pi}{2})。看起来十分简单，但其中会夹有伪素数。百科里面提到了根据n的值判断p(n)是否为伪素数，但程序实现复杂度较高。也可以用O(\sqrt n)的判断方法或Miller Rabin算法判断，但达不到我们想要的O(1)复杂度。

另一篇百度百科和一篇知乎都有关于素数公式的介绍，但这些公式计算起来都十分复杂，并没有降低计算素数的时间复杂度。

那么，在无法优化时间复杂度的情况下，我们自然地想起优化常数。

我们又想起了上面的判定质数的方法，有了一个筛质数的想法：能不能只筛掉除以6余1或5的数，减少冗余判断呢？
先从埃氏筛开始。埃氏筛的核心思想是：每遇到一个质数n，就筛掉比n\times n大的能整除n的所有合数。

#include<cstdio>
#include<bitset>
#include<cassert>
using namespace std;
const int N=1e8;
bitset<N+10>a;
bool ok=1;
int main(){
    //freopen("prime1.out","w",stdout);
    a[0]=a[1]=1;
    for(int i=2;i<=N;i++)
        if(!a[i]){
            //printf("%d\n",i);
            if(ok&&i*i>N)ok=0;//ok防止i乘法上溢
            if(ok)for(int j=i*i;j<=N;j+=i)a[j]=1;
        }
    return 0;
}

不难看出，每次筛掉的质数构成一个等差数列。它具备两个因素：首项与公差。我们只筛掉除以6余1或5的数时，也应该考虑这两个条件。
下面详细分析一下素数筛常数优化的方法：首先，每次被筛掉的数必定是当前考虑的质数n的倍数，而这个被筛掉的数模6应该为1或5，所以，我们把可能要被筛掉的数分为两类：模6余1与模6余5。这样，我们可以分别从两类数中找到n的最小倍数，然后因为6与任何质数互质，所以我们把公差设为6\times n，然后就可以筛了。
与n模6余数相同的那边总是很好找，因为n的最小倍数就是n（不过要从n\times 7开始筛）。但另一边呢？由同余我们想到了扩展欧几里得算法，在logn时间内可以解出同余方程6x+1\equiv0\pmod{n}（以模6余1那边为例），即6x\equiv n-1\pmod{n}。
有没有更简单更快的方法？有。
我们可以从不从x入手，而从n入手。因为x<n，所以我们可以考虑枚举n的满足条件的最小倍数kn中k的系数：枚举k为2、3、4、5即可。
我们又想到了最开始的定理：模6不为1或5的数为合数。我们可以联系到k的系数，k的系数明显不可能是2、3、4。这样，我们就只有一种可能了：k为5。
检验也很简单：设n=6i+1或n=6i-1，则第一种情况5n\%6=5，第二种情况5n\%6=1。
写代码就很简单了。注释部分供大家分析和思考。

#include<cstdio>
#include<bitset>
#include<cassert>
using namespace std;
const int N=1e8/6;
bool ok=1;
bitset<2>a[N+10];//a[i][0]表示%6=5的数(可以表示为6i-1)，a[i][1]表示%6=1的数(可以表示为6i+1) 
void work(int x,int d,int id){//筛数的部分，x为首项，d为公差 
    for(int i=x;i<=N;i+=d){
        a[i][id]=1;
    }
}
int main(){
    //#define OUT
    #ifdef OUT
    freopen("prime1.out","w",stdout);
    printf("2\n3\n");
    #endif
    for(int i=6,j=1;j<=N;i+=6,j++){
        if(5*i-1>N)ok=0;//如果不用输出所有质数这里就可以退出了
        if(!a[j][0]){
            int x=i-1;
            #ifdef OUT
            printf("%d\n",x);
            #endif
            assert((j+x)*6-1==x*7);//不会抛出异常
            if(ok){//ok只是为了加速，实际对程序安全性无影响，因为这里的j+x不会超过int
                work(j+x,x,0);
                work((x*5-1)/6,x,1);
                //work(5*i-1,x,1);与上一行等价
            }
        }
        if(!a[j][1]){
            int x=i+1;
            #ifdef OUT
            printf("%d\n",x);
            #endif
            assert((j+x)*6+1==x*7);//不会抛出异常
            if(ok){
                work(j+x,x,1);
                work((x*5+1)/6,x,0);
                //work(5*i+1,x,0);与上一行等价
            }
        }
    }
    return 0;
}

不过优化后埃筛速度提升没有欧拉筛快，于是我又检查了一遍代码，发现我推导首项那部分一开始就弄错了。最快的埃氏筛都是从n\times n开始筛的，我却还从n\times 2开始。
那就只好再按首项为n\times n推导一遍了。

1.(6i \pm 1)^2\%6=1

推导很简单，只看常数项就可以了。1\times1\%=5\times5\%6=1。

2.n\%6=1时，n\times(n+4)\%6=5；

看一下n\%6,(n+2)\%6,(n+4)\%6的值就明白了。

#include<cstdio>
#include<bitset>
#include<cassert>
using namespace std;
const int N=1e8/6,maxn=1e8;
bitset<2>a[N+10];//a[i][0]表示%6=5的数(可以表示为6i-1)，a[i][1]表示%6=1的数(可以表示为6i+1) 
bool ok=1;//又要用到避免溢出的技巧了
void work(int x,int d,int id){//筛数的部分，x为首项，d为公差，id标记x%6的余数 
    for(int i=x;i<=N;i+=d){
        a[i][id]=1;
    }
}
int main(){
    //#define OUT
    #ifdef OUT
    freopen("prime1.out","w",stdout);
    printf("2\n3\n");
    #endif
    for(int i=6,j=1;j<=N;i+=6,j++){
        if(!a[j][0]){
            int x=i-1;
            #ifdef OUT
            printf("%d\n",x);
            #endif
            if(ok&&x*x>maxn)ok=0;
            if(ok){
                work((x*x+1)/6,x,1);
                work((x*(x+2)+1)/6,x,0);
            }
        }
        if(!a[j][1]){
            int x=i+1;
            #ifdef OUT
            printf("%d\n",x);
            #endif
            if(ok){
                work((x*x+1)/6,x,1);
                work((x*(x+4)+1)/6,x,0);
            }
        }
    }
    return 0;
}

类似的，我们得出欧拉筛的优化方法。欧拉筛的核心是对任意正整数都筛一遍，碰到可以整除这个正整数的质数时停止。证明方法和上面类似。
下面是可能不完全严谨的证明：
先附原欧拉筛代码：

#include<cstdio>
#include<bitset>
using namespace std;
bitset<100000000+10>a;
int p[5761455+5],len,n;
int main(){
    scanf("%d",&n);
    a[0]=a[1]=1;
    for(int i=2;i<=n;i++){
        if(!a[i])p[++len]=i;
        for(int j=1;j<=len&&i*p[j]<=n;j++){
            a[i*p[j]]=1;
            if(i%p[j]==0)break;
        }
    }
    for(int i=1;i<=len;i++)printf("%d\n",p[i]);
    return 0;
}

可以发现，欧拉筛的核心部分是对每个数i（不管是质数还是合数），都在质数序列p[]中检查一遍，遇到i\%p[j]==0时停止。其中i\times p[j]分为两类：因数中有2或3的，和因数中没有2和3的。而且i的分类一定和i\times p[j]的分类相同（因为p[j]为质数）。所以由数学归纳法得：当我们不用因数中有2或3的数i的数跑欧拉筛时，也不会得到因数中有2或3的i\times p[j]，而其它因数中没有2和3的数还会被筛到。
所以我们得出了常数较小、效率较高的欧拉筛。

#include<cstdio>
#include<bitset>
#include<cassert>
using namespace std;
const int N=1e8/6,maxn=1e8;
bitset<2>a[N+10];
int p[5761455+5],len;//p数组存质数，5761455是计算得出的1e8以内的质数个数
void work(int x){
    for(int i=1;i<=len&&p[i]*x<=maxn;i++){
        int t=p[i]*x;
        int d=t%6==1;
        assert(a[(t+1)/6][d]==0);//不会抛出异常，代表所有数只被筛了一次
        a[(t+1)/6][d]=1;//t%6为1为5都有可能，不能确定
        if(x%p[i]==0)return ;
    }
}
int main(){
    //#define OUT
    #ifdef OUT
    freopen("prime1.out","w",stdout);
    printf("2\n3\n");
    //不要把这两个质数数预先加入p数组里，不过可以在p数组最前面预留两个位置，筛质数时不考虑最前面两个位置
    #endif
    for(int i=6,j=1;j<=N;i+=6,j++){
        if(!a[j][0])p[++len]=i-1;
        work(i-1);
        if(!a[j][1])p[++len]=i+1;
        work(i+1);
    }
    #ifdef OUT
    for(int i=1;i<=len;i++)printf("%d\n",p[i]);
    #endif
    return 0;
}

通过本地测试比较上面几种素数筛的效率，可以发现：

输出所有质数大约是筛质数时间的3倍或更多。
（下面不考虑程序输出和O2）
对埃氏筛与欧拉筛，优化后运行时间都为原来的40%。
优化后的埃氏筛比没有优化的欧拉筛快。
优化后的时间不是原来的1/3。我又试着减小优化后程序常数，比如去掉函数。（虽然代码更短，但代码可读性可能更差，不推荐编写）
不过这样的优化不管怎么样都不能达到原来1/3的时间。因为优化前后有的东西时间复杂度还是没有变，比如p数组优化前后都要把所有质数存下来。
但空间复杂度也有很大的优化：优化后埃氏筛直接变为原来1/3的空间，但优化后欧拉筛大于原来的1/3，因为优化前后都要有一个同样大小的数组存所有的质数。
洛谷IDE我也跑了一下：
我在洛谷IDE测开O2时n=10^8的情况，普通欧拉筛只要600ms左右，但优化欧拉筛总是有900ms以上。而不开O2两种方法都会超时。
貌似洛谷O2的优化效果貌似等同于几个指令集。
想暴力碾标算，n方过百万吗？就靠洛谷O2了

另外是几个优化：

1. 如果不需要处理出包含所有质数的数组p[]，只需要得到1~n的所有数的bool质数表，就可以在找到最后一个可以筛掉合数的数之后退出。如以6为基底（见下一条优化）时，此时只需要在5i>n时退出循环就可以了（5是大于2、3的质数）。如果n是确定的，还可以本地预处理出这个数i。虽然work函数都会把整个bool表格边里一遍，但主程序中就可以很快退出循环。
2. 还可以把6换成前k个素数的乘积n，如30、210等，以它们为基底筛素数。但由于要考虑的不仅是1~n以内的质数，还有另外一些非质数，如以210为基底时，121虽然是合数，但210i+121就可能是质数。所以实质上是要预处理出与基底互质的所有数。 所以我们只要控制好基底，就可以最小化素数筛的时间复杂度。这就是Wheel factorization筛的思想。