最近看了一些关于如何找出 1 ~ n 中所有质数的方法，跟大家分享一下。
最基本的方法是 sieve of Eratosthenes
令 is_prime 为一长度为 n + 1 的 bool 阵列 (index 0 起算)，
所有元素初值为 true 。
扫描 i = 2 ~ sqrt(n) 的整数，
如果 is_prime[i] 是 true ， 则 i 必为质数，
可以把 i 的所有小于 n 的倍数 j 设定 is_prime[j] = false
亦即 扫描 k = i ~ n / i 的整数, 令 j = k * i, 设定 is_prime[j] = false
而常见 sieve of Eratosthenes 的改进法有几种
1. 使用 bitset 而不是用 bool 阵列来节省空间。
2. 忽略 2 的倍数，可以再省去一半的空间和计算量，同样的道理可以忽略
3 的倍数、5的倍数等等，这概念就延伸成 wheel sieve [3] 。
3. 减少 is_prime[j] = false 被重复设定的次数。如果每一个非质数 j，
is_prime[j] 都只被设定一次，就是 linear sieve [4] 。
4. 增加 cache 的效能。
当然这些改进法有些是可以混合使用的，我这篇文章主要是在讨论 2 和 4 。
同时实际测试了 12 组不同的方法。
首先介绍减少 cache miss rate 的两个技巧，
第一个技巧比较简单，在 扫描 k = i ~ n / i 的整数 这一步中，
其实只需要考虑 k 是质数的情况即可，只有当 k 是质数才设定
is_prime[k * i] = false ，这样可以减少内存写入的次数。
当然在计算过程中，我们没办法知道 k 是否是质数，但是我们可以利用
is_prime[k] ，如果 is_prime[k] 是 false ，那 k 必不可能为质数，
此时就不用去设定 is_prime[k * i] = false 。
乍看之下这样作并没有改变复杂度，内层循环还是得执行 n / i - i + 1 次，
而且还多了一个条件判断，似乎不会变快。但是实际上因为存取 is_prime[k] 是
非常规则的，而且 k 的范围跟 j = k * i 的范围比起来小的多，更容易
被放在 cache 里面，所以利用小范围 is_prime[k] 的读取 + 判断，来
避免 cost 较高的 memory write ，是划得来的。另外内存的读取比
写入要来的有效率。所以使用这种技巧，执行时间可以变成普通
sieve of Eratosthenes 执行时间的 40% 左右。
另一个技巧叫做 segmented sieve [1] ，程式码可以参考这网页：
http://primesieve.org/segmented_sieve.html
是把 is_prime 分成多个 segment ，而每个 segment 的大小接近于
L1 data cache 的大小，然后每个 segment 都各筛一次。
复杂度本身是没变的，但是因为大幅降低了 cache miss rate ，速度
会比前一个技巧还要再更快。
接下来介绍如何跳过更多不可能为质数的数字 (wheel sieve)。
如果忽略所有偶数，在 bitset 中只需要表示所有奇数。
如果忽略所有 2 和 3 的倍数，在 bitset 中就只需要表示 6k+1, 6k-1 的数字。
同样的道理可以再忽略 5 的倍数。这么做的好处除了可以减少工作量之外，
还可以减少内存的使用。
虽然概念很简单，但是实作起来有点麻烦，因为要计算 index 。
如果只跳过偶数还算容易，如果还要跳过 3 的倍数和 5 的倍数，
就需要一些功夫。在实作上需要解决的问题有
1. 给定一个非 2, 3, 5 倍数的整数 i ，如何找出 i 在 bitset 中的 index 。
(bitset 中只表示非 2, 3, 5 的倍数）。
2. 给定一个非 2, 3, 5 倍数的整数 i 及其 index，
如何找出 > i 中最小的非 2, 3, 5 的倍数 i'。
3. 给定一个非 2, 3, 5 倍数的整数 i ，一个非 2, 3, 5 倍数的整数 k > i，
及 i, k, k * i 的 indices ，
如何找出 > k 中最小的非 2, 3, 5 倍数的整数 k' 及 k' * i 的 indices 。
而且这些计算都必须要尽量的使用 +- 运算和查表，避免 / 及 % ，因为 / 和 %
速度都很慢，甚至比 L2 cache miss 还慢。使用太多 / 和 % 会很没效率。
不过 / 和 % 一个常数是可以的，因为 compiler 可以把除法转成乘以倒数。
只跳过偶数的称为 2-wheel ，跳过 2 和 3 的倍数称为 6-wheel ，而跳过
2、3 和 5 的倍数称为 30-wheel 。同样的道理也可以跳过 2、3、5 和 7 的倍数
（210-wheel)，只是我实验的结果 30-wheel 比 2-wheel、6-wheel 和 210-wheel都
来的有效率（可能跟我的实作方法有关），所以接下来的实验部分就只
测试 30-wheel 了。
实验设计
我测试了 6 种不同的方法
sieve of Eratosthenes
sieve of Eratosthenes + 增加 cache 效能技巧 1
sieve of Eratosthenes + 增加 cache 效能技巧 2 (segmented sieve)
30-wheel sieve
30-wheel sieve + 增加 cache 效能技巧 2 (segmented sieve)
Linear sieve
及其相对应使用 bitset 的版本，所以总共有 12 个方法。
这 12 种方法全部都已预先排除偶数。
实验结果
实验结果中使用非 bitset 版本的 sieve of Eratosthenes 为基准，令其
执行时间为 1 单位。
而我得到的结果是当 n 比较小时（n <= 2^21）， 30-wheel sieve 最快。
当 n = 2^21 时，30-wheel sieve 的执行单位时间为 0.4 。
而 n 稍微大时（2^21 <= n <= 2^26），30-wheel sieve + bitset 最快。
当 n = 2^26 时，30-wheel sieve + bitset 的执行单位时间为 0.14 。
而 n 更大时，segmented 30-wheel sieve + bitset 。
当 n = 2^30 时， 执行的单位时间为 0.13 。
bitset
在我的实验中，使用 bitset 的版本未必会比未使用 bitset 的版本快，
当然这应该是因为 n 不够大（我只测到 2^30 ）所导致的。
在我的实验中，只有 sieve of Eratosthenes、30-wheel sieve和
segmented 30-wheel sieve 在使用 bitset 之后可以加速。
增加 cache 效能技巧 1
当 n = 2^30 时， sieve of Eratosthenes + 技巧 1 的执行单位时间为 0.4 。
Segmented sieve
当 n = 2^30 时， segmented sieve of Eratosthenes 的执行单位时间为 0.17。
Linear sieve
当 n = 2^30 时， linear sieve 的执行单位时间为 0.87
结论
文中介绍的这些技巧，其实程式码复杂度与 sieve of Eratosthenes 的程式码
复杂度差异不大，但是都是很有效的加速技巧。
基于 segmented sieve ，还有其他的改进法 [2] ，但是实作会比较复杂些。
[1] Carter Bays and Richard H. Hudson,
The segmented sieve of eratosthenes and primes in arithmetic
progressions to 10^12
BIT Numerical Mathematics June 1977, Volume 17, Issue 2, pp 121–127
[2] Emil Vatai, Sieving in primality testing and factorization
[3] Paul Pritchard,
Explaining the wheel sieve,
Acta Informatica 17 (1982), 477–485.
[4] Paul Pritchard
Linear prime-number sieves: A family tree
Science of Computer Programming Volume 9, Issue 1, August 1987,
Pages 17-35