在这篇帖子中,@wangchenyi 提出了一种质数筛法,但其复杂度并没有优于线性,实际上大致在欧拉筛和埃氏筛之间。而 @wkywkywky 给出的这篇论文中提出了一种时间复杂度为 O(loglognn),最低空间复杂度为 O(N0.5+o(1))≈O(60N) 的筛法。然后蒟蒻就实现了一下这个科技,并且其时间的确要优于一般的欧拉筛。
实际上我在实现该算法的时候,因为模板要求输出第 k 个质数,所以需要开一个 int 数组去存储所有质数;并且为了不对这些质数排序而浪费时间,把空间开到了 O(16×N0.5+o(1))≈O(960N) 级别。
目前在未经任何卡常的情况下,线性筛模板单点时间在 350ms 左右,空间在 33M 左右。相比一般的欧拉筛,无论是时间还是空间上都有较大的优化。
代码:
#include<iostream>
#include<cstdio>
#include<ctime>
#include<cstring>
using namespace std;
const int INF=2000000000;
const int N=100000000,B=560000;
int del1[9]={0,1,13,17,29,37,41,49,53};
int del2[5]={0,7,19,31,43};
int del3[5]={0,11,23,47,59};
int cnt,pri[6000001];
int pre[2001];bool ycl[11001];
bool bj1[9][B+5],bj2[5][B+5],bj3[5][B+5];
double be_time,en_time;
template<typename T> void Read(T &x)
{
x=0;bool f=1;
char ch=getchar();
while(ch<'0'||ch>'9') ch=getchar();
while(ch>='0'&&ch<='9')
{
x=x*10+(ch^48);
ch=getchar();
}
x=f==1?x:-x;
}
template<typename T> void Put(T x)
{
if(x<10) putchar(x+'0');
else Put(x/10),putchar(x%10+'0');
}
void be_speed_test(){be_time=clock();}
void en_speed_test(){en_time=clock();}
void show_time(){printf("%.3lfms\n",en_time-be_time);}
void init()
{
int n,tmp=0;Read(n);
for(int i=2;i<=110;i++)
if(!ycl[i])
for(int j=i;i*j<=11000;j++)
ycl[i*j]=1;
for(int i=7;i<=11000;i++)
if(!ycl[i]) pre[++tmp]=i;
pri[1]=2;pri[2]=3;pri[3]=5;cnt=3;
}
void get_pri(int L)
{
for(int i=1;i<=8;i++)
{
int d=del1[i];
for(int f=1;f<=15;f++)
for(int g=1;g<=30;g++)
if((4*f*f+g*g)%60==d)
{
int x=f,y0=g,k0=(4*f*f+g*g-d)/60;
while(k0<L+B) k0+=2*x+15,x+=15;
while(x>15)
{
x-=15;k0-=2*x+15;
while(k0<L) k0+=y0+15,y0+=30;
int k=k0,y=y0;
while(k<L+B) bj1[i][k-L]=1-bj1[i][k-L],k+=y+15,y+=30;
}
}
for(int q=1;pre[q]*pre[q]<60*(L+B);q++)
{
int now=pre[q]*pre[q],be=60*L/now+(60*L%now>0);
for(int j=be*now;j<=60*(L+B);j+=now)
if((j-d)%60==0) bj1[i][(j-d)/60-L]=0;
}
}
for(int i=1;i<=4;i++)
{
int d=del2[i];
for(int f=1;f<=10;f++)
for(int g=1;g<=30;g++)
if((3*f*f+g*g)%60==d)
{
int x=f,y0=g,k0=(3*f*f+g*g-d)/60;
while(k0<L+B) k0+=x+5,x+=10;
while(x>10)
{
x-=10;k0-=x+5;
while(k0<L) k0+=y0+15,y0+=30;
int k=k0,y=y0;
while(k<L+B) bj2[i][k-L]=1-bj2[i][k-L],k+=y+15,y+=30;
}
}
for(int q=1;pre[q]*pre[q]<60*(L+B);q++)
{
int now=pre[q]*pre[q],be=60*L/now+(60*L%now>0);
for(int j=be*now;j<=60*(L+B);j+=now)
if((j-d)%60==0) bj2[i][(j-d)/60-L]=0;
}
}
for(int i=1;i<=4;i++)
{
int d=del3[i];
for(int f=1;f<=10;f++)
for(int g=1;g<=30;g++)
if((3*f*f-g*g+6000)%60==d)
{
int x=f,y0=g,k0=(3*f*f-g*g-d)/60;
while(1)
{
int en=0;
while(k0>=L+B)
{
if(x<=y0){en=1;break;}
k0-=(y0+15),y0+=30;
}
if(en) break;
int k=k0,y=y0;
while(k>=L&&y<x) bj3[i][k-L]=1-bj3[i][k-L],k-=(y+15),y+=30;
k0+=x+5,x+=10;
}
}
for(int q=1;pre[q]*pre[q]<60*(L+B);q++)
{
int now=pre[q]*pre[q],be=60*L/now+(60*L%now>0);
for(int j=be*now;j<=60*(L+B);j+=now)
if((j-d)%60==0) bj3[i][(j-d)/60-L]=0;
}
}
for(int i=0;i<B;i++)
{
int cur1=1,cur2=1,cur3=1;
int cur=60*(L+i);
while(cur1<=8||cur2<=4||cur3<=4)
{
int num1=INF,num2=INF,num3=INF;
if(cur1<=8) num1=del1[cur1];
if(cur2<=4) num2=del2[cur2];
if(cur3<=4) num3=del3[cur3];
if(num1<num2&&num1<num3)
{
if(bj1[cur1][i]) pri[++cnt]=cur+num1,bj1[cur1][i]=0;
cur1++;
}
else if(num2<num1&&num2<num3)
{
if(bj2[cur2][i]) pri[++cnt]=cur+num2,bj2[cur2][i]=0;
cur2++;
}
else
{
if(bj3[cur3][i]) pri[++cnt]=cur+num3,bj3[cur3][i]=0;
cur3++;
}
}
}
}
void cal()
{
for(int i=0;i<=N/(60*B);i++) get_pri(i*B);
}
int main()
{
init();
// be_speed_test();
cal();
// en_speed_test();
// show_time();
int k,x;
Read(k);
while(k--)
{
Read(x);
Put(pri[x]);putchar('\n');
}
return 0;
}