夺取不了线性筛素数最优解手把手教程 / 浅谈埃氏筛优化

Jason331 · 2026-02-04 19:34:36 · 算法·理论

\gdef\dis{\displaystyle}

前言

你说的对，标题和某文章很像，但是真的推荐夺取多项式全家桶最优解手把手教程。

你说的对，但是我不想用线性筛，我要用埃氏筛。

下面的内容假设读者有至少一个大脑。

埃氏筛

题目：P3383 【模板】线性筛素数。

题目：U658905 【模板】埃氏筛素数（n = 10^9），本篇文章以上一题为例，此题仅提供更大范围的测试，读者可以尝试修改代码后通过此题。

基本方法

我们首先建立一个标记数组，然后从小到大枚举每个质数，将这个数所有比自己大的倍数标记为合数。

时间复杂度（不想证）：O(n\log\log n)。

于是你写出了一个小巧清新的代码：

:::success[code]

#include <bits/stdc++.h>
using namespace std;

constexpr int N = 100000010;
bool is_prime[N];
int prime[N],ptop;

void premake()
{
    memset(is_prime,true,sizeof is_prime);
    is_prime[0] = is_prime[1] = false;
    ptop = 0;
    for(int i = 2;i < N;i++)
    {
        if(!is_prime[i]) continue;
        prime[++ptop] = i;
        for(int j = i * 2;j < N;j += i)
            is_prime[j] = false;
    }
}

int main()
{
    premake();
    int n,q;
    scanf("%d %d",&n,&q);
    while(q--)
    {
        int x;
        scanf("%d",&x);
        printf("%d\n",prime[x]);
    }
    return 0;
}

:::

虽然 AC 了，但是一看耗时：总用时 10.35s！这也太慢了！

先来几个常见优化

我们发现埃氏筛之所以比线性筛复杂度劣，是因为每个数字会被它的所有质因子筛一遍，有没有办法减少不必要的筛选？

一、平方 - 平方根优化（名字瞎起的）

对于每一个质数 p，我们现在要筛掉 kp，如果 k < p，那 kp 肯定在之前被 k（或其质因子）筛过了，所以我们可以从 p^2 开始筛选。

又发现，我们既然从 p^2 开始筛选，那我们只需要枚举小于等于 \sqrt n 的质数即可。

二、bitset

可以使用 bitset 来代替 bool 数组，这样既能节省空间，又可以减小常数。

:::info[为什么埃氏筛使用 bitset 可以优化时间，而线性筛常出现负优化？] 说说我个人的理解：比较线性筛和埃氏筛的算法逻辑，发现对于每个数筛掉的数字数量来说，线性筛是较均匀的，而埃氏筛则明显偏向于小数字。

而这就导致了：埃氏筛有更多较为连续的内存访问，表现为，小质数筛数时隔了几个数就会访问一次；而线性筛的内存访问则非常随机，表现为大部分的筛选都由比较大的数字完成，两次访问的内存间隔较大。

可以用下面的代码测试：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

constexpr int N = 100000010;

struct node
{
    bool data[N];
    int last,cnt;
    ll sum;
    void set(){last = cnt = sum = 0;memset(data,true,sizeof data);}
    bool &operator[](int note)
    {
        cnt++,sum += abs(last - note);
        return data[last = note];
    }
}is_prime;

int prime[N],ptop;

pair <ll,int> Eratosthenes1()
{
    is_prime.set();
    is_prime[0] = is_prime[1] = false;
    ptop = 0;
    for(int i = 2;i < N;i++)
    {
        if(!is_prime[i]) continue;
        for(int j = i + i;j < N;j += i)
            is_prime[j] = false;
    }
    for(int i = 0;i < N;i++)
        if(is_prime[i])
            prime[++ptop] = i;
    return make_pair(is_prime.sum,is_prime.cnt);
}

pair <ll,int> Eratosthenes2()
{
    is_prime.set();
    is_prime[0] = is_prime[1] = false;
    ptop = 0;
    for(int i = 2;i * i < N;i++)
    {
        if(!is_prime[i]) continue;
        for(int j = i * i;j < N;j += i)
            is_prime[j] = false;
    }
    for(int i = 0;i < N;i++)
        if(is_prime[i])
            prime[++ptop] = i;
    return make_pair(is_prime.sum,is_prime.cnt);
}

pair <ll,int> Euler()
{
    is_prime.set();
    for(int i = 2;i < N;i++)
    {
        if(is_prime[i])
            prime[++ptop] = i;
        for(int j = 1;j <= ptop && i * prime[j] < N;j++)
        {
            is_prime[i * prime[j]] = false;
            if(i % prime[j] == 0)
                break;
        }
    }
    return make_pair(is_prime.sum,is_prime.cnt);
}

int main()
{
    printf("start\n");
    auto ans1 = Eratosthenes1(),ans2 = Euler(),ans3 = Eratosthenes2();
    printf("Eratosthenes1: gap sum -> %18lld, times -> %10d, average %.6f\n",
            ans1.first,ans1.second,(double)ans1.first / ans1.second);
    printf("Eratosthenes2: gap sum -> %18lld, times -> %10d, average %.6f\n",
            ans3.first,ans3.second,(double)ans3.first / ans3.second);
    printf("Euler        : gap sum -> %18lld, times -> %10d, average %.6f\n",
            ans2.first,ans2.second,(double)ans2.first / ans2.second);
    return 0;
}

运行结果：

可以看到欧拉筛虽然访问数组的次数（times）比两种埃氏筛都要少，但是访问间隔之和（gap sum）和访问平均间隔（average）都远大于两种埃氏筛。

还可以看到而经过上面“平方 - 平方根优化”的 Eratosthenes2 不仅优化了访问次数，还显著优化了平均访问间隔。

但这跟 bitset 有什么关系呢？

在计算机中，局部的密集访问带来的常数是远小于大范围随机访问的。

我们知道，bitset 单点修改和访问的常数肯定是劣于 bool 数组的。

在欧拉筛中，平均访问间隔太大，bitset 对于空间局部密集访问的常数优化并没有那么明显，即使在内存中的间隔变为了原来的 \frac 1 8，平均间隔还是很大。于是时间消耗的差别便落在了单点修改和访问的常数上。

而在埃氏筛中，平均访问间隔较小，经过 bitset 的空间优化，可以让更多的访问变得非常密集，平均间隔缩小到只有几十个 B，这样对缓存十分友好，极大改善了寻址的常数，这个优化甚至强到盖过了单点修改和访问的常数缺陷，于是 bitset 便显得比 bool 数组快了。 :::

于是你写下了以下小巧清新的代码：

:::success[code]

#include <bits/stdc++.h>
using namespace std;

constexpr int N = 100000010;
bitset <N> is_prime;
int prime[N],ptop;

void premake()
{
    is_prime.set();
    is_prime[0] = is_prime[1] = false;
    ptop = 0;
    for(int i = 2;i * i <= N;i++)
    {
        if(!is_prime[i]) continue;
        for(int j = i * i;j < N;j += i)
            is_prime[j] = false;
    }
    for(int i = 0;i <= N;i++)
        if(is_prime[i])
            prime[++ptop] = i;
}

int main()
{
    premake();
    int n,q;
    scanf("%d %d",&n,&q);
    while(q--)
    {
        int x;
        scanf("%d",&x);
        printf("%d\n",prime[x]);
    }
    return 0;
}

:::

2.63s，还可以。

往后的优化就不太常用了~~因为码量开始增大了~~。

三、分块优化

众所可能周知，埃氏筛是可以分块的，只需要先预处理出小于等于 \sqrt n 的质数，然后对于每一块枚举每个质数进行筛数即可。

需要注意的是，这并不会改善时间复杂度，在大多数情况下也不会劣化时间复杂度（除非你的块太小了），这只是能对内存更加友好，进一步改善改善空间常数。

提供的代码选用块长为 2^{18}，效率会比较优秀。

然后我们顺便手写一个 bitset，为后面的优化做准备；代码中顺便还根据实际的质数上限重新调整了质数列表的大小。

:::success[code]

#include <bits/stdc++.h>
using namespace std;
typedef unsigned long long ull;
typedef long long ll;
constexpr int N = 100000000,P = 6000000,S = 10000,W = 64,B = 18;

struct node
{
    ull data[1 << B - 6];
    void clear() {memset(data,0xff,sizeof data);}
    bool get(int pos) {return data[pos >> 6] & (1ull << (pos & 63));}
    void set(int pos) {data[pos >> 6] &= (ull)-1 ^ 1ull << (pos & 63);}
    void solve(int l);
}tool;

int pl[P],ptop,rtop;

void node::solve(int l)
{
    int e = l + (1 << B);
    clear();
    for(int i = 1;pl[i] * pl[i] < e;i++)
        for(int j = max(pl[i],(l - 1) / pl[i] + 1) * pl[i];j < e;j += pl[i])
            set(j - l);
    for(int i = 0;i < 1 << B;i++)
        if(get(i))
            pl[++ptop] = i + l;
}

void premake()
{
    tool.clear();
    for(int i = 2;i <= S;i++)
        if(tool.get(i))
            for(int j = i * i;j < 1 << B;j += i)
                tool.set(j);
    for(int i = 2;i < 1 << B;i++)
        if(tool.get(i))
            pl[++ptop] = i;
    for(int l = 1 << B;l < N;l += 1 << B)
        tool.solve(l);
}

int main()
{
    premake();
    int n,q;
    scanf("%d %d",&n,&q);
    for(int x;q--;)
    {
        scanf("%d",&x);
        printf("%d\n",pl[x]);
    }
    return 0;
}

:::

2.34s。

四、质数收集优化

注意到我们筛完之后需要重新遍历 bitset，收集所有素数，但我们可以只取为 true 的位——通过 lowbit 或 builtin 系列函数实现。

这样便把收集质数部分的时间复杂度从 O(n) 降低到了 \dis O(\frac {n} {\ln n} + \frac n w)，当然整体的时间复杂度还是没变。

:::success[code]

#include <bits/stdc++.h>
using namespace std;
typedef unsigned long long ull;
typedef long long ll;
constexpr int N = 100000000,P = 6000000,S = 10000,W = 64,B = 18;

struct node
{
    ull data[1 << B - 6];
    void clear() {memset(data,0xff,sizeof data);}
    bool get(int pos) {return data[pos >> 6] & (1ull << (pos & 63));}
    void set(int pos) {data[pos >> 6] &= (ull)-1 ^ 1ull << (pos & 63);}
    void solve(int l);
}tool;

int pl[P],ptop,rtop;

void node::solve(int l)
{
    int e = l + (1 << B);
    clear();
    for(int i = 1;pl[i] * pl[i] < e;i++)
        for(int j = max(pl[i],(l - 1) / pl[i] + 1) * pl[i];j < e;j += pl[i])
            set(j - l);
    for(int i = 0;i < 1 << B - 6;i++)
        while(data[i])
        {
            int b = __builtin_ctzll(data[i]);
            pl[++ptop] = l + (i << 6) + b;
            data[i] ^= 1ull << b;
        }
}

void premake()
{
    tool.clear();
    tool.set(0),tool.set(1);
    for(int i = 2;i <= S;i++)
        if(tool.get(i))
            for(int j = i * i;j < 1 << B;j += i)
                tool.set(j);
    for(int i = 0;i < 1 << B - 6;i++)
        while(tool.data[i])
        {
            int b = __builtin_ctzll(tool.data[i]);
            pl[++ptop] = (i << 6) + b;
            tool.data[i] ^= 1ull << b;
        }
    for(int l = 1 << B;l < N;l += 1 << B)
        tool.solve(l);
}

int main()
{
    premake();
    int n,q;
    scanf("%d %d",&n,&q);
    for(int x;q--;)
    {
        scanf("%d",&x);
        printf("%d\n",pl[x]);
    }
    return 0;
}

:::

1.80s。

五、bitset 进阶应用——小质数提前筛除

考察数字 2，发现我们可以在初始化时将 bitset 初始化为 010101010101\dots 这样，这可以保证所有的偶数都被提前筛除了。

其他的小质数能不能这么做呢？当然是可以的。

对于每个 0 \le i \lt p，我们预处理出一个 64 位的 unsigned long long 变量，其二进制表示下从第 i 位开始，每隔 p 位就有一个 0，其他位都是 1，具体的，运行以下代码：

ull fl[64];
memset(fl,-1,sizeof fl);
for(int i = 0;i < p;i++)
    for(int j = i;j < 64;j += p)
        fl[i] ^= 1ull << j;

这样我们就可以对于 bitset 中的每个 unsigned long long 变量 x，根据其代表的首位数字除以 p 的余数找出一个合适的预处理变量，让这个变量和 x 进行位与操作，筛除 x 所代表的区间内所有 p 的倍数。

到底找哪个预处理变量呢？留给读者自行推敲。（逃

如此我们便可在 \dis O(\frac n w) 的时间复杂度下完成本需 \dis O(\frac n p) 的时间复杂度才能搞定的任务。

但是由于这种操作单次操作常数稍大，提供的代码中只对前 10 个小质数使用了这种技巧。

:::success[code]

#include <bits/stdc++.h>
using namespace std;
typedef unsigned long long ull;
typedef long long ll;
constexpr int N = 100000000,P = 6000000,S = 10000,W = 64,B = 18,Y = 10;

struct node
{
    ull data[1 << B - 6];
    void clear() {memset(data,0xAA,sizeof data);}
    bool get(int pos) {return data[pos >> 6] & (1ull << (pos & 63));}
    void set(int pos) {data[pos >> 6] &= (ull)-1 ^ 1ull << (pos & 63);}
    void solve(int l);
    void fliter(int l,int p);
}tool;

int pl[P],ptop,rtop;

void node::fliter(int l,int p)
{
    ull fl[64];
    memset(fl,-1,sizeof fl);
    for(int i = 0;i < p;i++)
        for(int j = i;j < 64;j += p)
            fl[i] ^= 1ull << j;
    int start = (p - l % p) % p;
    const int x = 64 % p;
    for(int i = 0,j = start;i < 1 << B - 6;i++)
    {
        data[i] &= fl[j];
        if((j -= x) < 0) j += p;
    }
}

void node::solve(int l)
{
    int e = l + (1 << B);
    clear();
    for(int i = 2;i <= Y;i++)
        fliter(l,pl[i]);
    for(int i = Y + 1;pl[i] * pl[i] < e;i++)
        for(int j = max(pl[i],(l - 1) / pl[i] + 1) * pl[i];j < e;j += pl[i])
            set(j - l);
    for(int i = 0;i < 1 << B - 6;i++)
        while(data[i])
        {
            int b = __builtin_ctzll(data[i]);
            pl[++ptop] = l + (i << 6) + b;
            data[i] ^= 1ull << b;
        }
}

void premake()
{
    memset(tool.data,-1,sizeof tool.data);
    tool.set(0),tool.set(1);
    for(int i = 2;i <= S;i++)
        if(tool.get(i))
            for(int j = i * i;j < 1 << B;j += i)
                tool.set(j);
    for(int i = 0;i < 1 << B - 6;i++)
        while(tool.data[i])
        {
            int b = __builtin_ctzll(tool.data[i]);
            pl[++ptop] = (i << 6) + b;
            tool.data[i] ^= 1ull << b;
        }
    for(int l = 1 << B;l < N;l += 1 << B)
        tool.solve(l);
}

int main()
{
    premake();
    int n,q;
    scanf("%d %d",&n,&q);
    for(int x;q--;)
    {
        scanf("%d",&x);
        printf("%d\n",pl[x]);
    }
    return 0;
}

:::

1.15s~~，好的我们现在代码长度来到了 2k~~。

六、快读

终于该轮到快读了，我们可以随意使用一些快读，这里仅使用 fread、fwrite 实现的快读，想要更快？这里推荐浅谈 C++ IO 优化。

:::success[code]

#include <bits/stdc++.h>
using namespace std;
typedef unsigned long long ull;
typedef long long ll;
constexpr int N = 100000000,P = 6000000,S = 10000,W = 64,B = 18,Y = 10;

constexpr int PSIZ = 40,ISIZ = 1 << 16,OSIZ = ISIZ;
char ibuf[ISIZ],*ip1 = ibuf,*ip2 = ibuf,obuf[OSIZ],*outp1 = obuf,*outp2 = obuf + OSIZ,st[PSIZ];
unsigned char ltop;
void flush() {fwrite(obuf,1,outp1 - obuf,stdout);outp1 = obuf;}
char getc() {return ip1 == ip2 && (ip2 = (ip1 = ibuf) + fread(ibuf,1,ISIZ,stdin)),ip1 == ip2?-1:*ip1++;}
void putc(char c) {*outp1++ = c;if(outp1 == outp2) flush();}
void read(int &x) {char c;for(c = getc();c < '0';c = getc());for(x = 0;c >= '0';c = getc())x = x * 10 + (c ^ '0');}
void write(int x) {ltop = 0;do st[ltop++] = x % 10 + '0';while(x /= 10);while(ltop--) putc(st[ltop]);}

struct node
{
    ull data[1 << B - 6];
    void clear() {memset(data,0xAA,sizeof data);}
    bool get(int pos) {return data[pos >> 6] & (1ull << (pos & 63));}
    void set(int pos) {data[pos >> 6] &= (ull)-1 ^ 1ull << (pos & 63);}
    void solve(int l);
    void fliter(int l,int p);
}tool;

int pl[P],ptop,rtop;

void node::fliter(int l,int p)
{
    ull fl[64];
    memset(fl,-1,sizeof fl);
    for(int i = 0;i < p;i++)
        for(int j = i;j < 64;j += p)
            fl[i] ^= 1ull << j;
    int start = (p - l % p) % p;
    const int x = 64 % p;
    for(int i = 0,j = start;i < 1 << B - 6;i++)
    {
        data[i] &= fl[j];
        if((j -= x) < 0) j += p;
    }
}

void node::solve(int l)
{
    int e = l + (1 << B);
    clear();
    for(int i = 2;i <= Y;i++)
        fliter(l,pl[i]);
    for(int i = Y + 1;pl[i] * pl[i] < e;i++)
        for(int j = max(pl[i],(l - 1) / pl[i] + 1) * pl[i];j < e;j += pl[i])
            set(j - l);
    for(int i = 0;i < 1 << B - 6;i++)
        while(data[i])
        {
            int b = __builtin_ctzll(data[i]);
            pl[++ptop] = l + (i << 6) + b;
            data[i] ^= 1ull << b;
        }
}

void premake()
{
    memset(tool.data,-1,sizeof tool.data);
    tool.set(0),tool.set(1);
    for(int i = 2;i <= S;i++)
        if(tool.get(i))
            for(int j = i * i;j < 1 << B;j += i)
                tool.set(j);
    for(int i = 0;i < 1 << B - 6;i++)
        while(tool.data[i])
        {
            int b = __builtin_ctzll(tool.data[i]);
            pl[++ptop] = (i << 6) + b;
            tool.data[i] ^= 1ull << b;
        }
    for(int l = 1 << B;l < N;l += 1 << B)
        tool.solve(l);
}

int main()
{
    premake();
    int n,q;
    read(n);read(q);
    for(int x;q--;)
        read(x),write(pl[x]),putc('\n');
    flush();
    return 0;
}

:::

896ms。

七、奇偶优化

注意到除了 2 以外的偶数都不是质数，所以我们可以直接不存所有偶数，只进行奇数的筛除。

注意，这里代码实现难度增大了，要让不存偶数的 bitset 适用前面几个优化，需要读者仔细思考，精细实现，尤其注意分块时的余数处理。

:::success[code]

#include <bits/stdc++.h>
using namespace std;
typedef unsigned long long ull;
typedef long long ll;
constexpr int N = 100000000,P = 6000000,S = 10000,W = 64,B = 18,Y = 10;
constexpr int prel[Y + 1] = {0,2,3,5,7,11,13,17,19,23,29};

constexpr int PSIZ = 40,ISIZ = 1 << 16,OSIZ = ISIZ;
char ibuf[ISIZ],*ip1 = ibuf,*ip2 = ibuf,obuf[OSIZ],*outp1 = obuf,*outp2 = obuf + OSIZ,st[PSIZ];
unsigned char ltop;
void flush() {fwrite(obuf,1,outp1 - obuf,stdout);outp1 = obuf;}
char getc() {return ip1 == ip2 && (ip2 = (ip1 = ibuf) + fread(ibuf,1,ISIZ,stdin)),ip1 == ip2?-1:*ip1++;}
void putc(char c) {*outp1++ = c;if(outp1 == outp2) flush();}
void read(int &x) {char c;for(c = getc();c < '0';c = getc());for(x = 0;c >= '0';c = getc())x = x * 10 + (c ^ '0');}
void write(int x) {ltop = 0;do st[ltop++] = x % 10 + '0';while(x /= 10);while(ltop--) putc(st[ltop]);}

struct node
{
    ull data[1 << B - 6];
    void clear() {memset(data,-1,sizeof data);}
    bool get(int pos) {return data[pos >> 6] & (1ull << (pos & 63));}
    void set(int pos) {data[pos >> 6] &= (ull)-1 ^ 1ull << (pos & 63);}
    void solve(int l);
    void fliter(int l,int p);
}tool;

int pl[P],ptop,rtop;

void node::fliter(int l,int p)
{
    ull fl[64];
    memset(fl,-1,sizeof fl);
    for(int i = 0;i < p;i++)
        for(int j = i;j < 64;j += p)
            fl[i] ^= 1ull << j;
    /*
    real = l * 2 + 1
    mod = ((real + p) % (2 * p) / 2)
    ans = -mod = p - (real + p) % (2 * p) / 2
    */
    int start = (p - (l * 2 + 1 + p) % (2 * p) / 2) % p;
    const int x = 64 % p;
    for(int i = 0,j = start;i < 1 << B - 6;i++)
    {
        data[i] &= fl[j];
        if((j -= x) < 0) j += p;
    }
}

void node::solve(int l)
{
    int e = l + (1 << B);
    clear();
    for(int i = 2;i <= Y;i++)
        fliter(l,prel[i]);
    for(int i = Y + 1;pl[i] * pl[i] < e * 2 + 1;i++)
    {
        int start = pl[i] * pl[i] / 2;
        int div = (l * 2 + 1 + pl[i] - 1) / (2 * pl[i]) + 1;
        start = max(start,(div * pl[i] * 2 - pl[i]) / 2);
        for(int j = start;j < e;j += pl[i])
            set(j - l);
    }
    for(int i = 0;i < 1 << B - 6;i++)
        while(data[i])
        {
            int b = __builtin_ctzll(data[i]);
            pl[++ptop] = (l + (i << 6) + b) << 1 | 1;
            data[i] ^= 1ull << b;
        }
}

void premake()
{
    memset(tool.data,-1,sizeof tool.data);
    tool.set(0),tool.set(1);
    for(int i = 2;i <= S;i++)
        if(tool.get(i))
            for(int j = i * i;j < 1 << B;j += i)
                tool.set(j);
    for(int i = 0;i < 1 << B - 6;i++)
        while(tool.data[i])
        {
            int b = __builtin_ctzll(tool.data[i]);
            pl[++ptop] = (i << 6) + b;
            tool.data[i] ^= 1ull << b;
        }
    for(int l = 1 << B - 1;l * 2 + 1 < N;l += 1 << B)
        tool.solve(l);
}

int main()
{
    premake();
    int n,q;
    read(n);read(q);
    for(int x;q--;)
        read(x),write(pl[x]),putc('\n');
    flush();
    return 0;
}

:::

687ms。

八、三次方预判

对于一个质数 p，在埃氏筛的过程中如果它筛掉了一个数 kp \lt p^3，那么 k 一定也是一个质数。

原因是显然的，如果 k 是合数，设 k 可以被分解为 xy，那么因为 xy = k \lt p^2，所以有 \min(x,y) \lt p，这说明 kp 之前已经被更小的质数筛掉了。

于是我们可以让每个质数在其三次方内只筛自己的质数倍。

:::success[code]

#include <bits/stdc++.h>
using namespace std;
typedef unsigned long long ull;
typedef long long ll;
constexpr int N = 100000000,P = 6000000,S = 10000,W = 64,B = 18,Y = 10;
constexpr int prel[Y + 1] = {0,2,3,5,7,11,13,17,19,23,29};

constexpr int PSIZ = 40,ISIZ = 1 << 16,OSIZ = ISIZ;
char ibuf[ISIZ],*ip1 = ibuf,*ip2 = ibuf,obuf[OSIZ],*outp1 = obuf,*outp2 = obuf + OSIZ,st[PSIZ];
unsigned char ltop;
void flush() {fwrite(obuf,1,outp1 - obuf,stdout);outp1 = obuf;}
char getc() {return ip1 == ip2 && (ip2 = (ip1 = ibuf) + fread(ibuf,1,ISIZ,stdin)),ip1 == ip2?-1:*ip1++;}
void putc(char c) {*outp1++ = c;if(outp1 == outp2) flush();}
void read(int &x) {char c;for(c = getc();c < '0';c = getc());for(x = 0;c >= '0';c = getc())x = x * 10 + (c ^ '0');}
void write(int x) {ltop = 0;do st[ltop++] = x % 10 + '0';while(x /= 10);while(ltop--) putc(st[ltop]);}

struct node
{
    ull data[1 << B - 6];
    void clear() {memset(data,-1,sizeof data);}
    bool get(int pos) {return data[pos >> 6] & (1ull << (pos & 63));}
    void set(int pos) {data[pos >> 6] &= (ull)-1 ^ 1ull << (pos & 63);}
    void solve(int l);
    void fliter(int l,int p);
}tool;

int pl[P],ptop,rtop,sl[S];

inline void node::fliter(int l,int p)
{
    ull fl[64];
    memset(fl,-1,sizeof fl);
    for(int i = 0;i < p;i++)
        for(int j = i;j < 64;j += p)
            fl[i] ^= 1ull << j;
    int start = (p - (l * 2 + 1 + p) % (2 * p) / 2) % p;
    const int x = 64 % p;
    for(int i = 0,j = start;i < 1 << B - 6;i++)
    {
        data[i] &= fl[j];
        if((j -= x) < 0) j += p;
    }
}

inline void node::solve(int l)
{
    const int e = l + (1 << B),real_end = e << 1 | 1;
    clear();
    for(int i = 2;i <= Y;i++)
        fliter(l,prel[i]);
    {
        int i;
        for(i = Y + 1;pl[i] * pl[i] * pl[i] < real_end;i++)
        {
            if(pl[sl[i]] < pl[i] * pl[i])
            {
                while(pl[sl[i]] < pl[i] * pl[i] && pl[i] * pl[sl[i]] < real_end)
                    set(pl[sl[i]++] * pl[i] / 2 - l);
                if((ll)pl[i] * pl[i] * pl[i] >= real_end) continue;
            }
            int start = pl[i] * pl[i] * pl[i] >> 1;
            int div = (l * 2 + 1 + pl[i] - 1) / (2 * pl[i]) + 1;
            start = max(start,(div * pl[i] * 2 - pl[i]) >> 1);
            for(int j = start;j < e;j += pl[i]) set(j - l);
        }
        for(;pl[i] * pl[i] < real_end;i++)
            while(pl[sl[i]] < pl[i] * pl[i] && pl[i] * pl[sl[i]] < real_end)
                set(pl[sl[i]++] * pl[i] / 2 - l);
    }
    for(int i = 0;i < 1 << B - 6;i++)
    {
        const int cnt = __builtin_popcountll(data[i]);
        for(int j = 0;j < cnt;j++)
        {
            int b = __builtin_ctzll(data[i]);
            pl[++ptop] = (l + (i << 6) + b) << 1 | 1;
            data[i] ^= 1ull << b;
        }
    }
}

inline void premake()
{
    memset(tool.data,-1,sizeof tool.data);
    tool.set(0),tool.set(1);
    for(int i = 2;i <= S;i++)
        if(tool.get(i))
            for(int j = i * i;j < 1 << B;j += i)
                tool.set(j);
    for(int i = 0;i < 1 << B - 6;i++)
        while(tool.data[i])
        {
            int b = __builtin_ctzll(tool.data[i]);
            pl[++ptop] = (i << 6) + b;
            tool.data[i] ^= 1ull << b;
        }
    for(int i = Y + 1;pl[i] * pl[i] < N;i++)
    {
        sl[i] = i;
        while(pl[sl[i]] < pl[i] * pl[i] && pl[sl[i]] * pl[i] < 1 << B)
            sl[i]++;
    }
    for(int l = 1 << B - 1;l * 2 + 1 < N;l += 1 << B)
        tool.solve(l);
}

int main()
{
    premake();
    int n,q;
    read(n);read(q);
    for(int x;q--;)
        read(x),write(pl[x]),putc('\n');
    flush();
    return 0;
}

:::

589ms。

总结

好的，到这里我们已经可以在 60ms 内预处理出 10^8 内的所有质数，稍微改一改上一份代码，就可以在 600ms 内预处理出 10^9 内的所有素数。但可惜这篇文章到这里也该结束了，后面大概就是：

:::warning[不要点开] 你发现五个点输入输出要三百多毫秒而预处理只要二百多毫秒然后开始狂卡输入输出然后搞到总用时四百多毫秒代码长度再翻一倍最后自己不想卡了去看最优解发现都是超标怪输入输出到处是你无法理解的神奇代码但是跑起来飞快然后你一气之下吐血而亡才怪其实是你直接崩溃三百六十度螺旋升天然后决定写下一篇文章引诱大家来卡常。 :::

好的，希望这篇文章中的思想与技巧能够对大家有所启发，虽然没能做到最优解，但相信在卡常的路上会有收获。

后记

这篇文章将会在我的 P7884 【模板】Meissel-Lehmer 题解（可能因为做法逆天没过审）中被引用。

本文目前仅在洛谷平台上发表，如果你在其他平台（包括但不限于 CSDN）上看见了此文章，请尽快联系洛谷 uid1126439 的用户并说明恶意转载情况。

参考资料

大部分比较鬼畜的优化方法是作者自己想的，但是还是感谢有以下资料进行参考。

https://www.luogu.com.cn/discuss/635288

https://www.luogu.com.cn/article/x1qt8l1k