起因是看到了这个帖子如果高中生能证明哥德巴赫猜想，会被清华北大数学系保送吗？ - 知乎，emmm…
本文的目标是筛选 32 位无符号整数（约 40 亿）中的所有素数，并在此范围内验证哥德巴赫猜想。
试除法
一开始并没有打算寻找所有素数，而是暴力判断是否为素数。一个小技巧是将自然数表示为 6k + b，当中只有 6k + 1 和 6k + 5 的形式可能为素数。也可以用 30k + b 等，简单起见不再优化了。
后来发现我还是 naive，验证所有的 32 位无符号整数耗时太长，验证到 1000 万都要 10 秒左右。
bool prime(uint32_t number) {
if (number < 6)
return number == 2 || number == 3 || number == 5;
if (number % 6 != 1 && number % 6 != 5)
return false;
for (uint32_t i = 5, sq = (uint32_t) sqrt(number); i <= sq; i += 6)
if (!(number % i) || !(number % (i + 2)))
return false;
return true;
}埃拉托斯特尼筛法
如果能先求得所有素数，然后查表，速度会有很大提升。
对于每个整数使用一个 bit 来标识是否为素数，32 位无符号整数需要 2^32 / 8 = 512 MiB 空间，剔除所有的偶数，需要一半即 256 MiB。如果剔除 6k + 1 和 6k + 5 以外的数，只需要 1/3 的空间，但是不利于映射时使用位操作，故只剔除偶数。
经了解，筛选素数常见的有两种方法：埃拉托斯特尼筛法和欧拉筛法。前者思路简单，且无需保存素数序列，时间复杂度为 $O(n \log \log n)$；后者需要保存素数序列，时间复杂度为 $O(n)$。
所谓埃氏筛法，就是筛去 $\le \sqrt{n}$ 的素数 $2, 3, 5, …$ 的倍数，维基百科的图片清晰易懂：
bool prime(char *data, uint32_t n) {
if (n == 2)
return true;
if (!(n & 1))
return false;
return !(data[n >> 4] & (1 << (n >> 1 & 7)));
}
void sieve_of_eratosthenes(char *data) {
data[0] = 1;
uint32_t n, sq, start, c, step;
for (n = 3, sq = (uint32_t) sqrt(UINT32_MAX); n <= sq; n += 2) {
if (data[n >> 4] & (1 << (n >> 1 & 7)))
continue;
for (start = n * n, step = n << 1, c = start; c >= start; c += step)
data[c >> 4] |= 1 << (c >> 1 & 7);
}
}这里除了位操作外，还用到一些优化：
- 当筛选素数 $n$ 的倍数时，$2n, 3n, …, (n-1)n$ 必然已经筛选过了，故只需从 $n^2$ 开始，于是筛选序列为 $n^2, n^2 + n, n^2 + 2n…$
- 当 $k$ 为奇数时，由于 $n$ 为奇数，则 $n^2 + kn$ 必为偶数，算法只对奇数保存标识，故 $k$ 只需取偶数，筛选序列只剩下 $n^2, n^2 + 2n, n^2 + 4n…$
在 i7-7700HQ + 16 GiB 内存的 MacBook Pro 上，开启 O3 优化后耗时约 22 秒：
Sieve of Eratosthenes cost 21.723315 seconds.
Same story.
Verification cost 63.079718 seconds.
这个算法的主要缺陷是，不利于 CPU 利用程序的局部性原理。筛选的时候 data 数组并不是按顺序读写的，而是类似列优先读写二维数组那样，cache 命中率太低，大量耗时都用在了读写内存上。
分段埃氏筛法
cache 命中率对程序速度有很大影响，于是有了分段筛选。
在 32 位无符号整数范围内，埃氏筛法需要筛掉 sqrt(UINT32_MAX) 以内的素数的倍数。根据素数计数函数，在此范围内的素数约有 8000 个，保存只需占用约 32 KiB，不再像上面那样查询 data 数组。
将 32 位无符号整数划分为若干个等长的段，先使用埃氏筛法筛选第一个段，并收集 sqrt(UINT32_MAX) 以内的素数，对其它段筛选前若干个素数的倍数。数据段不宜过大，否则 cache 命中率会降低，但也不宜过小，否则复杂度也会增大。
i7-7700HQ 的 L1 数据 cache 大小为 32 KiB，L2 cache 为 256 KiB，L3 cache 为 6 MiB。在实践中，发现选择 128 KiB 作为一个段大小较为合适。
// Note: the size must be power of 2.
#define SEGMENT_SIZE 0x200000
// 32 KiB L1 data cache = 2^15*2^3 = 2^18 bits, identifying 2^19 numbers.
// #define SEGMENT_SIZE 0x80000
// 256 KiB L2 cache
// #define SEGMENT_SIZE 0x400000
// 6 MiB L3 cache, use 2 MiB of it
// #define SEGMENT_SIZE 0x2000000
bool prime(char *data, uint32_t n) {
if (n == 2)
return true;
if (!(n & 1))
return false;
return !(data[n >> 4] & (1 << (n >> 1 & 7)));
}
/*
* Run sieve of Eratosthenes in [0, max(SEGMENT_SIZE, sqrt(UINT32_MAX))),
* and store prime numbers in [3, sqrt(UINT32_MAX)].
*/
uint32_t initial_sieve(char *data, uint32_t *primes) {
const uint32_t SQ = (uint32_t) sqrt(UINT32_MAX), LIMIT = SEGMENT_SIZE > SQ ? SEGMENT_SIZE : SQ;
data[0] = 1;
uint32_t n, sq, c, step, count = 0;
for (n = 3, sq = (uint32_t) sqrt(LIMIT); n <= sq; n += 2) {
if (data[n >> 4] & (1 << (n >> 1 & 7)))
continue;
for (c = n * n, step = n << 1; c <= LIMIT; c += step)
data[c >> 4] |= 1 << (c >> 1 & 7);
}
for (n = 3; n <= SQ; n += 2) {
if (!(data[n >> 4] & (1 << (n >> 1 & 7))))
primes[count++] = n;
}
return count;
}
void segmented_sieve_of_eratosthenes(char *data) {
// store prime numbers in [3, sqrt(UINT32_MAX)]
uint32_t *primes = (uint32_t *) malloc(8192 * sizeof(uint32_t));
const uint32_t count = initial_sieve(data, primes);
uint32_t sq, i, p, c, step;
const uint32_t SEGMENTS = (uint32_t)(1 << 31) / (SEGMENT_SIZE >> 1);
// first segment has been sieved
for (uint32_t s = 1, low = SEGMENT_SIZE, high = (SEGMENT_SIZE << 1) - 1;
s < SEGMENTS; ++s, low += SEGMENT_SIZE, high += SEGMENT_SIZE) {
sq = (uint32_t) sqrt(high);
for (i = 0; i < count && (p = primes[i]) <= sq; ++i) {
// c = p^2 + 2k * p
c = low / p * p;
if (c < low) c += p;
if (!(c % 2)) c += p;
// uint32t overflow in last segment
for (step = p << 1; low <= c && c <= high; c += step)
data[c >> 4] |= 1 << (c >> 1 & 7);
}
}
free(primes);
}速度有了明显的提升：
Segmented sieve of Eratosthenes cost 5.576287 seconds.
Same story.
Verification cost 63.089097 seconds.
这里对每个素数都计算其倍数 c，也可以将所有素数的倍数保存到一个数组中，与 primes 数组一样只占 32 KiB，这样每次无需重新计算。但在实践中发现性能反而下降了一点，推测是仅仅几千个数的乘除和加法 * 段数，耗时并不多，而内存读写已经超过了这个耗时。
欧拉筛法
在埃氏筛法中，一个合数会被多次筛选到，因此它不是线性的算法。欧拉筛法使每个合数只被其最小质因数筛去一次，其复杂度为 $O(n)$，在数量很大时有很大优势。
bool prime(char *data, uint32_t n) {
if (n == 2)
return true;
if (!(n & 1))
return false;
return !(data[n >> 4] & (1 << (n >> 1 & 7)));
}
void sieve_of_euler(char *data) {
uint32_t *primes = (uint32_t *) malloc(8192 * sizeof(uint32_t));
const uint32_t SQ = (uint32_t) sqrt(UINT32_MAX), LIMIT_N = UINT32_MAX / 3;
uint32_t count = 0, n, p, c, i, LIMIT_P;
data[0] = 1;
for (n = 3; n <= LIMIT_N; n += 2) {
if (n <= SQ && !(data[n >> 4] & (1 << (n >> 1 & 7))))
primes[count++] = n;
LIMIT_P = UINT32_MAX / n;
for (i = 0; i < count && (p = primes[i]) <= LIMIT_P; ++i) {
c = p * n;
data[c >> 4] |= 1 << (c >> 1 & 7);
if (!(n % p))
break;
}
}
free(primes);
}关键在当 n 是素数 p 的倍数时直接 break。
算法的正确性（所有合数都能被筛去）：对于任意合数 $c$，将其写成质因数乘积 $c = p_1^{k_1} p_2^{k_2}… p_j^{k_j}$，$c$ 一定能在 $n = p_1^{k_1 - 1} p_2^{k_2}… p_j^{k_j}$ 被筛去，因为 $n \ge p_1$，即素数 $p_1$ 必然已经被收集。
算法的复杂度（合数只被最小质因数筛去一次）：假设合数 $c = p_1^{k_1} p_2^{k_2}… p_x^{k_x}… p_j^{k_j}$ 也被 $p_x$ 筛去，此时 $n = p_1^{k_1} p_2^{k_2}… p_x^{k_x - 1}… p_j^{k_j}$。若 $n < p_x$，此时 $p_x$ 还没被收集，不符合；若 $n \ge p_x$，之前的循环中已经枚举了 $p_1, p_2…, p_i$，在 $p_1$ 的时候就 break 了，故不符合。
虽然是线性复杂度，实际运行好于埃氏筛法，但比起分段埃氏筛法还是慢很多：
Sieve of Euler cost 15.921141 seconds.
Same story.
Verification cost 63.412436 seconds.
而欧拉筛法本身又不适合改写成分段算法，看来 cache 的威力还是相当可观的。
总结
试除法适用于数据较小的情况，埃氏筛法和欧拉筛法适用于数据较大的情况。如果对性能要求较高，可以考虑分段埃氏筛法。
