Optimize Prime Sieve

Posted 2020-12-05 tmzbot

/// A heavily optimized sieve
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
typedef unsigned int u32;
typedef unsigned long long ull;
const char pr60[]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59};
const char masks[][4]={
    {3,7,11,13},
    {3,17,19,23},
    {2,29,31},
    {2,37,41},
    {2,43,47},
    {2,53,59}
};
const u32 segsize=65536;
void Apply_mask(u32*a,u32*b,u32 l1,u32 l2){
    u32 t=0;
    for(u32 q=0,r=l1/l2;q<r;++q)
        for(u32 i=0;i<l2;++i)
            a[t++]|=b[i];
    for(u32 i=0;t<l1;++i)
        a[t++]|=b[i];
}
void Gen_mask_sub(u32*a,u32 l1,u32 b){
    u32 st=b>>1,rt=0;
    while(rt<l1){
        a[rt]|=1<<st;
        st+=b;
        if(st>=30)st-=30,++rt;
        if(st>=30)st-=30,++rt;
    }
}
void PrintMask(u32*a,u32 len){
    printf("Mask of len %u
",len*60);
    for(u32 i=0;i<len;++i){
        for(u32 j=0;j<30;++j)
            if((a[i]&(1<<j)))
                printf("%llu
",i*60ull+j*2ull+1ull);
    }
}
u32 Gen_mask(u32*a,int id){
    int len=masks[id][0];
    u32 ll=1;
    for(int i=1;i<=len;++i)
        ll*=masks[id][i];
    memset(a,0,4*ll);
    for(int i=1;i<=len;++i)
        Gen_mask_sub(a,ll,masks[id][i]);
//  PrintMask(a,ll);
    return ll;
}
const u32 mask=0x1a4b3496;
const u32 pr60_m=0xdb4b3491;
u32 pr[10000][4],prl;
int main(){
    ull ma,tma,tmx;scanf("%llu",&ma);
    tma=(ma-1)/60+1;
    tmx=tma*60;//upper limit
    u32*sieve=new u32[tma];// getting a sieve ready
    u32*maske=new u32[7429];
    std::fill(sieve,sieve+tma,mask);
    for(int i=0;i<6;++i)
        Apply_mask(sieve,maske,tma,Gen_mask(maske,i));

    ull preseg=std::min(tmx,ull(sqrt(ma)/60)+1);
    u32 j=61;
    for(;ull(j)*j<=preseg*60;j+=2){
        u32 v=j/60,u=(j%60)>>1;
        if(!(sieve[v]&(1<<u))){
            v=j/30,u=j%30;
            u32 rt=j*3/60,st=(j*3%60)>>1;
            while(rt<preseg){
                sieve[rt]|=1<<st;
                rt+=v;
                st+=u;
                if(st>=30)st-=30,++rt;
            }
            pr[prl][0]=v;
            pr[prl][1]=u;
            pr[prl][2]=rt;
            pr[prl][3]=st;
            prl++;
        }
    } // Non-segmented sieve core
    if(preseg==tmx)goto end;
    for(u32 segl=preseg;segl<tma;segl+=segsize){
        u32 segr=std::min(segl+segsize,u32(tma));
        for(;ull(j)*j<=segr*60;j+=2){
            u32 v=j/60,u=(j%60)>>1;
            if(!(sieve[v]&(1<<u))){
                v=j/30,u=j%30;
                ull t=j*ull(j);
                u32 rt=t/60,st=t%60>>1;
                pr[prl][0]=v;
                pr[prl][1]=u;
                pr[prl][2]=rt;
                pr[prl][3]=st;
                prl++;
            }
        }
        for(int i=0;i<prl;++i){
            u32 v=pr[i][0],u=pr[i][1],rt=pr[i][2],st=pr[i][3];
            while(rt<segr){
                sieve[rt]|=1<<st;
                rt+=v;
                st+=u;
                if(st>=30)st-=30,++rt;
            }
            pr[i][0]=v;
            pr[i][1]=u;
            pr[i][2]=rt;
            pr[i][3]=st;
        }
    }
    end:
    sieve[0]=pr60_m;
    int count=1;
    for(u32 i=0;i<tma;++i){
        for(u32 j=0;j<30;++j)
            if(!(sieve[i]&(1<<j)))++count;
    }
    for(ull a=tmx-1;a>ma;a-=2){
        u32 i=a/60,j=a%60>>1;
        if(!(sieve[i]&(1<<j)))--count;
    }
    printf("%d
",count);
    return 0;
}

一个Eratosthenes筛。单线程，筛1e9<0.5s，程序运行时间<0.7s。(7700k 4.2GHz)
（注意后面的统计部分效率极其低下，可用查表位运算替代。）

0.只留奇数项
1.压位,30pack int
2.对于<60的素数用很快的筛子一遍
3.分段筛法