zoukankan      html  css  js  c++  java
  • Optimize Prime Sieve

    /// 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.分段筛法

  • 相关阅读:
    P3387 【模板】缩点 tarjan
    P2831 愤怒的小鸟 状压dp
    交流帖
    P3959 宝藏 模拟退火。。。
    B1060 [ZJOI2007]时态同步 dfs
    P1850 换教室 概率dp
    树链刨分(待修改)
    B3403 [Usaco2009 Open]Cow Line 直线上的牛 deque
    B3402 [Usaco2009 Open]Hide and Seek 捉迷藏 最短路
    B5248 [2018多省省队联测]一双木棋 状压dp
  • 原文地址:https://www.cnblogs.com/tmzbot/p/9326890.html
Copyright © 2011-2022 走看看