PTT數位生活區 / Prob_Solve (計算數學 Problem Solving)

[心得] 篩法的加速技巧

看板Prob_Solve (計算數學 Problem Solving)作者FRAXIS (喔喔)時間8年前 (2016/08/07 07:07)推噓7(7推 0噓 1→)

留言8則, 7人參與討論串1/1

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