[心得] 篩法的加速技巧

[FRAXIS (喔喔)]

最近看了一些關於如何找出 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://www.ptt.cc/bbs/Prob_Solve/M.1470524853.A.577.html
推 prisonbreak: 推推！
※ 編輯: FRAXIS (24.23.200.172), 08/07/2016 22:49:28
推 yr: 只能跪了
推 cocoyan: 滿有趣的！
推 NaiveRed: 推!
推 hioska: 推
推 obelisk0114: 記得最快的是 Sieve of Atkin,不過那三條公式一直不
→ obelisk0114: 知道怎麼推出來的
推 jxzhe: 推~