如果你想品尝“美味”的算法讲解，请点击 这里 。

由于$60 = 2^2 \times 3 \times 5$，因此所有大于 5 的素数都满足 $\gcd(n, 60) = 1$。

60 以下的与 60 互素的剩余类共有 $\phi(60) = 16$ 个：

Atkin 把这 16 个数分成三组：

为什么这样分呢？ 因为每组恰好对应一个二次方程在模 60 下能表示的数。

我们列出三个二次方程：

我们需要了解算法主要使用的定理：

阿特金筛法定理：对于 $\gcd(n,60)=1$ 且不含平方因子的 $n$：

• 若 $n \bmod 60 \in S_1$，则 $n$ 是素数 $\iff$ 方程①的正整数解个数为奇数
• 若 $n \bmod 60 \in S_2$，则 $n$ 是素数 $\iff$ 方程②的正整数解个数为奇数
• 若 $n \bmod 60 \in S_3$，则 $n$ 是素数 $\iff$ 方程③的正整数解个数为奇数

千万别想着证明，这需要大学代数数论。其实是我不会讲 但不过我们可以验证一下：

例如 $n=13$ （$13 \bmod 60 = 13 \in S_1$，应该用方程①）：

• $x=1$: $4+y^2=13 => y^2=9 => y=3$ => 解 $(1,3)$
• $x=2$: $16+y^2=13$ => 无解
• $x \geq 2$ 后 $4x^2 > 13$，停止

解的个数为 1 ， 1 是奇数， 得出 13 是素数，与事实符合。

再例如 $n=49$（$49 \bmod 60 = 49 \in S_1$ ）：

• $x=1$: $y^2=45$ => 无整数解
• $x=2$: $y^2=33$ => 无整数解
• $x=3$: $y^2=13$ => 无整数解
• $x=4$: $4x^2=64>49$ ，结束

解的个数为 0 ，正确判断为非素数。

注意到 $49 = 7^2$ 含平方因子。有些含平方因子的数解的个数确实是奇数，这就是为什么算法需要额外剔除含平方因子的数子。定理只对无平方因子数保证等价性。

代码

#include<bits/stdc++.h>
using namespace std;
constexpr int M = 5e8 + 10;

// b: 素性标记数组，只存储奇数
// b[i] 代表数字 (2*i+1) 是否为素数候选
bitset<((M + 1) >> 1)> b;

inline void Atkin(int n) {//O(N/log^2N) 
    assert(n <= M);//这是报错代码，如果n<=M为false，程序就会发出错误信息。
	 
    // 这里用 (n+1)>>1 作为循环上界(不含)
    n = ((n + 1) >> 1);

    // 二次型翻转阶段：利用 XOR 翻转统计解的奇偶性
    // 所有循环都做了两重优化：
    //   1. if(var%3): 模60预过滤，跳过使 n 被3整除的情况
    //   2. 内层用二阶等差数列增量代替乘法，避免重复计算二次型

    // 【二次型 3x^2+y^2, x为奇数的情形】
    // y 从1开始步长2(只取奇数)，初始 m=(y^2+36)/2 是该y下第一个合法x对应的下标
    // 内层 k+=36 是 x 步长为6时的二阶差分(/2后)
    // m+=(k+=36)+18 生成增量序列 18,54,90,... 对应原始增量 36,108,180,...
    for (int y = 1, m; (m = (y * y + 36) >> 1) < n; y += 2)
        if (y % 3)  // y不被3整除 => n=3x2+y^2不被3整除
            for (int k = 0; m < n; m += (k += 36) + 18)
                b.flip(m);  // XOR翻转，累计解的奇偶性

    // 【二次型 4x^2+y^2】
    // x 从1开始遍历，y固定从1起步(y必为奇数)
    // 初始 m=(4x2+1)/2 对应 y=1
    // k+=4 是 y 步长2时 4y+4 的二阶差分(/2后=4)
    // m+=(k+=4) 生成增量 4,8,12,16,... 对应原始增量 8,16,24,32,...
    for (int x = 1, m; (m = (4 * x * x + 1) >> 1) < n; x++)
        if (x % 3)  // x不被3整除 => n=4x2+y^2不被3整除
            for (int k = 0; m < n; m += (k += 4))
                b.flip(m);

    // 【二次型 3x^2+y^2, y为偶数(x为奇数)的情形】
    // 与第一个循环互补，共同覆盖 3x^2+y^2 的全部合法解
    // y 从2开始步长2(只取偶数)
    // k+=12 是 x 步长2时 12x+12 的二阶差分
    for (int y = 2, m; (m = (y * y + 3) >> 1) < n; y += 2)
        if (y % 3)
            for (int k = 0; m < n; m += (k += 12))
                b.flip(m);

    // 【二次型 3x2-y^2 (x>y)】
    // 初始 m=((2y+6)*y+3)/2 = (2y2+6y+3)/2
    //   对应 x=y+1 时的最小合法值: 3(y+1)2-y^2 = 2y2+6y+3
    // k 初始为 6y，然后 k+=12
    //   对应 x 步长2时增量 12x+12 在下标空间的表达
    for (int y = 1, m; (m = ((2 * y + 6) * y + 3) >> 1) < n; y++)
        if (y % 3)
            for (int k = 6 * y; m < n; m += (k += 12))
                b.flip(m);

    // 平方剔除阶段：清除含平方因子的误报合数
    // 只对已标记为True的p操作(真素数)，清除p^2的所有倍数
    // p从5开始(2,3的平方倍数不与60互素，从未被标记)
    for (int p = 5, q; (q = p * p) >> 1 < n; p += 2)
        if (b[p >> 1])  // p是素数才需要剔除
            // m 从 p^2/2 开始，步长 q=p^2
            // 下标空间中 p^2 的倍数间隔恰好是 p^2 (不是p^2/2)
            // 因为 (k*p^2-1)/2 - ((k-1)*p^2-1)/2 = p^2/2... 
            // 实际清除的是 p^2, 3p^2, 5p^2...
            // 偶数倍 p^2 对应偶数，不在奇数bitset中，自然跳过
            for (int m = q >> 1; m < n; m += q)
                b.reset(m);  // 设为false

    // 手动标记数字3 (下标1) 为素数
    // 前面的二次型未覆盖到3，且n>1时才标记
    if (n>1) b.set(1);
}
int l, a[1000010];
int N, A, B;

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);
    
    cin >> N >> A >> B;
    Atkin(N);
    
    int c = 0;// c: π(N)的计数器，同时用于系数索引

    // 单独处理素数2（偶数，不在奇数bitset中）
    // 若2<=N，则c++，并检查是否为常数项B的元素
    // c++%A==B：该素数是第(c%A)个系数元素，恰好等于目标常数项B
    if (2 <= N && c++ % A == B)
        a[l++] = 2;

    // 遍历所有奇数
    for (int n = 3; n <= N; n += 2) {
        if (b[n >> 1] && c++ % A == B) // b[n>>1]: 查bitset判断素性
            a[l++] = n;
    }
    // 循环结束后 c = π(N)
    cout << c << " " << (c + A - 1 - B) / A << "\n";
    for (int i = 0; i < l; i++)
        cout << a[i] << " ";
}

算法运行极快，时间复杂度为 $O(\frac{N}{\log^2N})$ 。