最近看了一些关於如何找出 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 --

※ 发信站: 批踢踢实业坊(ptt.cc), 来自: 24.23.200.172 ※ 文章网址: https://webptt.com/cn.aspx?n=bbs/Prob_Solve/M.1470524853.A.577.html ※ 编辑: FRAXIS (24.23.200.172), 08/07/2016 22:49:28