演算法：極速搜尋質數
2024/09/09
|
演算法 |
|
極速搜尋質數 |
|
作者：韓特野
|
歷史修訂記錄：
2024/09/09：初稿。
2024/10/31：演算法的效能分析 新增2項。
2024/11/04：新增章節 另一個演算法：建立大質數表 。
2024/11/24: 新增章節 計算機實際測試的數據
2024/12/15: 修改章節 計算機實際測試的數據
2024/12/17: 修改章節 計算機實際測試的數據
前言
質數，指大於1的自然數中，除了1和該數自身以外，無法被其他自然數整除的數。本文件將敘述如何完全不用除法，卻能夠快速逐步搜尋質數。
演算法的初始構思
首先，我們看下圖所示，假設已知質數為2、3、5以及7，每個大於2的自然數所包含的質因數皆做了標記，而未被標記的自然數11、13、17和19便是質數。
|
自然數 |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
10 |
11 |
12 |
13 |
14 |
15 |
16 |
17 |
18 |
19 |
… |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
2 |
|
|
|
＊ |
|
＊ |
|
＊ |
|
＊ |
|
＊ |
|
＊ |
|
＊ |
|
＊ |
|
… |
|
|
3 |
|
|
|
|
|
＊ |
|
|
＊ |
|
|
＊ |
|
|
＊ |
|
|
＊ |
|
… |
|
|
5 |
|
|
|
|
|
|
|
|
|
＊ |
|
|
|
|
＊ |
|
|
|
|
… |
|
|
7 |
|
|
|
|
|
|
|
|
|
|
|
|
|
＊ |
|
|
|
|
|
… |
|
|
. |
. |
. |
. |
. |
. |
. |
. |
. |
. |
. |
. |
. |
. |
. |
. |
. |
. |
. |
. |
|
|
此圖看似埃拉托斯特尼篩法（sieve of Eratosthenes）的應用。然而，當利用連結串列（linked list）把這些屬於同一個合數的質因數連接起來時，將乍然發現另一條搜尋質數之坦道。
如下圖所示，假設已知質數為自然數2，並且將2掛至自然數4（2的平方）的位置。
自2的下一個自然數3開始依序判斷後續的自然數是否為質數。步驟如下：
1. 自然數3沒有質因數，將3加入質數表內，並且將3掛在自然數9（3的平方）的位置。
2. 自然數4有質因數2，將2移到6的位置（原來的位置加上該質數值：4+2）。
3.
自然數5沒有質因數，將5加入質數表內，並且將5掛在自然數25（5的平方）的位置。
4. 自然數6有質因數2 ，將2移到8的位置（原來的位置加上該質數值：6+2）。
5. 自然數7沒有質因數，將7加入質數表內，並且將7掛在自然數49（7的平方）的位置。
6. 自然數8有質因數2 ，將2移到10的位置（原來的位置加上該質數值：8+2）。
7.
自然數9有質因數3 ，將3移到12的位置（原來的位置加上該質數值：9+3）。
8.
自然數10有質因數2 ，將2移到12的位置（原來的位置加上該質數值：10+2）。
9. 自然數11沒有質因數，將11加入質數表內，並且將11掛在自然數121（11的平方）的位置。
10. 自然數12有質因數2 和3。分別將2移到14的位置（原來的位置加上該質數值：12+2），3移到15的位置（原來的位置加上該質數值：12+3）。
11.
自然數13沒有標記，將13加入質數表內，並且將13掛在自然數169（13的平方）的位置。
如上所述的步驟逐步搜尋質數，這就是本文件的探尋質數之路。
演算法：建立質數表
為了便於闡述演算法以及其效能分析，筆者定義了相關的資料結構以及類似C語言的程式碼，僅供參考。
相關的資料結構
為了簡化演算法，使其簡明易懂，假設所有變數的最大值以及陣列變數長度皆沒有上限。
n
陣列變數naturalTable[]: naturalTable[n]代表自然數n的第一個質因數的索引值。起始值皆為0。
n
變數detectingNumber：正受檢測的自然數值。
n
變數primesFound：已搜尋到的質數個數。
n
陣列變數primeTable[]: 質數表，primeTable[n]代表第n個質數值。
n
陣列變數primeNext[]: 當自然數包含多個質因數時，需要陣列變數primeNext[]來指向下一個質因數。起始值皆為0。
公用資料結構初始值設定
primesFound = 1；
primeTable[1]
= 2；
naturalTable[4]
= 1；
detectingNumber
= 3；
演算法程式碼
void BuildPrimeTable() // 生成質數的主程式碼。
{
while(true)
{
if (naturalTable[detectingNumber] == 0)
{
IncreasePrimeTable(); // 將detectingNumber加到質數表, 並且掛到其平方值的位置。
}
else
{
MovePrimeFactors2NextPosition(); // 將所有質因數移至下一個位置。
}
detectingNumber++;
}
}
void IncreasePrimeTable()
{
primeTable[++primesFound] = detectingNumber;
naturalTable[detectingNumber * detectingNumber] = primesFound;
}
void
MovePrimeFactors2NextPosition ()
{
ulong primeIdx = naturalTable[detectingNumber]; // 假設ulong值沒有上限。
while (primeIdx > 0)
{
SetPrimePosition(primeIdx, detectingNumber + primeTable[primeIdx]);
primeIdx = primeNext[primeIdx];
}
naturalTable[detectingNumber] = 0；
}
// 將質因數掛在指定的naturalTable[]位置。
void
SetPrimePosition(ulong idx, ulong position)
{
primeNext[idx] = naturalTable[position];
naturalTable[position] = idx;
}
優化演算法
此演算法非常適合利用計算機來搜尋質數，然而搜尋大質數時必須仰賴硬碟來儲存大量資料結構內容，因此硬碟的讀寫頻率將影響整個運行效率。由於任何自然數只要檢測到存在1個質因數就可判斷不是質數，所以將質因數移到下一個位置時，可額外判斷該位置是否已經存在其它質因數。如果新位置存在其它質因數，則可繼續往下一個位置移動，直到最新位置不存在任何質因數為止。如此可大大降低硬碟的寫入動作，並且降低處理合數的平均時間。譬如質數2每檢測2個自然數就得移動1次，太浪費硬碟寫入動作。如果陣列naturalTable[]的長度夠大，則最理想情況是每個合數只連接1個質因數。
程式碼優化後更改如下：
void
MovePrimeFactors2NextPosition()
{
ulong primeIdx = naturalTable[detectingNumber]; // 假設ulong值沒有上限。
while (primeIdx > 0)
{
ulong newPosition = detectingNumber + primeTable[primeIdx];
while (naturalTable[newPosition] > 0 ) newPosition += primeTable[primeIdx];
SetPrimePosition(primeIdx, newPosition);
primeIdx = primeNext[primeIdx];
}
naturalTable[detectingNumber] = 0；
}
演算法的效能分析
n 完全不需要使用除法器，並且搜尋到的皆是確定型質數。
n 試除法需要試除所有n開根號範圍內的所有質數，才可確認n為質數，相當缺乏效率，且隨著n值越大越耗時間。反之，此演算法只需檢查每個整數的對應值為0，就可直接判斷n為質數。
n 對於非質數的檢測，此演算法的效率亦不遑多讓。以2的64次方範圍為例，2^64 <
2*3*5*7*11*13*17*19*23*29*31*37*41*43*47*53，需要處理（移動）的質因數最多15個。以2的128次方範圍為例，2^128 < 2*3*5*7*11*13*17*19*23*29*31*37*41*43*47*53*59*61*67*71*73*79*83*89*97*101，需要處理（移動）的質因數最多25個。隨著自然數n值遞增，所需檢測的最大質因數個數的比例呈現遞減情況。
n AKS質數測試是一個無仰賴質數的確定型質數測試演算法，並且可以在多項式時間內确定一個被測試數是否為質數。此演算法可以快速判斷超大整數是否為質數，其運算速度為 O(log12n) ,後來被改進到O(log6n)，無數數學家不斷研究如何提高演算法的速度。然而，如果利用它來判斷指定範圍n內的所有質數，其整體運算速度仍不及篩除法。
n
雖然極速搜尋法與埃拉托斯特尼篩法都是以篩法的概念來搜尋給定範圍內的所有質數，但是兩者處理步驟不同，極速搜尋法更適合在計算機上使用。分析結果如下表所示：
|
|
埃拉托斯特尼篩法 |
極速搜尋法 |
|
篩法步驟 |
1. 首先選定所要搜尋的質數範圍n，從2開始將其所有倍數篩除。 2. 將下一個未被篩除的整數（即為質數）的所有倍數篩除。 3. 重複第2個步驟，直到下一個未被篩除的整數值大於等於 “n的方根” 為止。 |
如本文所述。 |
|
運算效能 |
log n |
log n |
|
差異性 |
一次性篩除所有倍數。 |
逐一篩除。 |
|
適用在計算機運算的優劣性 |
1. 只需使用1個位元（bit）來判斷每個整數是否為質數。 2. 由於一次性篩除所有倍數，導致記憶體的再使用率缺乏彈性，能搜尋的質數範圍受到較大限制。 |
1. 使用索引值來判斷每個整數是否為質數，佔用較大的記憶體空間。 2. 由於使用逐一篩除法，可輕易循環重複使用記憶體。在記憶體以及硬碟的容量限制下，較容易定義相關的資料結構，發展更快速的程式碼。 |
另一個演算法：建立大質數表
上述的資料結構架構可依據已建立的質數来搜尋整數區間[a,b]的所有質數，a和b皆可透過人機界面來手動設置。a > primeTable[primesFound]，b < (primeTable[primesFound]+2)2（primeTable[primesFound]+2可能也是質數）。
為了與上述的質數表做區分，本文件稱此演算法所搜尋的質數為大質數，建立的質數表為大質數表。
此演算法有兩大主要用途：
1.
可使用較小的儲存體來建立大質數表。
2.
鑒於計算機的運算能力有限，可分配給不同的計算機同時搜尋不同整數區間[a,b]的所有大質數。
為了便於闡述演算法，筆者定義了相關的資料結構以及類似C語言的程式碼，僅供參考。
相關的資料結構
為了簡化演算法，使其簡明易懂，假設所有變數的最大值以及陣列變數長度皆沒有上限。
基本上，建立質數表以及大質數表的資料結構完全相同，只是宣告的變數以及陣列名稱不同而已。
n
變數a，b：上述的整數區間[a, b]起始值和最大值。
n
陣列變數largeNaturalTable[]: largeNaturalTable[n]代表自然數(n+a)的第一個質因數索引值。起始值皆為0。
n
變數largeDetectingNumber：正受檢測的自然數值。
n
變數largePrimesFound：已搜尋到的大質數個數。
n
陣列變數largePrimeTable[]: 質數表，largePrimeTable[n]代表第n個大質數值。
n
陣列變數largePrimeNext[]: 當自然數包含多個質因數時，需要陣列變數largePrimeNext[]來指向下一個質因數。起始值皆為0。
公用資料結構初始值設定
largePrimesFound = 0；
largeDetectingNumber
= a；
演算法程式碼
建立大質數表分成2個階段：
階段1：設定相關資料結構的起始值。
階段2：開始建立大質數表。
由於階段2的演算法與建立質數表的演算法完全相同，本文不再贅述相關程式碼。階段1只需做一件非常簡單卻特別耗時的事：將b的方根範圍內的所有質數逐一從a值開始計算每個質數的第一個倍數，並且掛在對應的largeNaturalTable[]位置。
void InitLargeData()
// 階段1的主程式碼。
{
ulong primeIdx = 1, primeValue =
primeTable[primeIdx];
ulong rem;
while(primeValue * primeValue <= b)
{
rem = a % primeValue;
if (rem == 0) SetLargePrimePosition(primeIdx, 0);
else SetLargePrimePosition(primeIdx, primeValue - rem);
primeValue = primeTable[++primeIdx];
}
}
// 將質因數掛在指定的largeNaturalTable[]位置。具優化空間。
void SetLargePrimePosition(ulong
idx, ulong position)
{
largePrimeNext[idx] = largeNaturalTable[position];
largeNaturalTable[position]
= idx;
}
計算機實際測試的數據
根據上述的2個演算法，筆者開發一套僅適用在Windows系統應用程式“質數表生成器”，必須配置64位元中央處理器。如果配置2G可用記憶體，則可建立大於2^64的大質數表。經測試結果顯示：在2^40範圍內每秒至少可檢測2千5百萬個數字，包含將搜尋到的質數儲存至硬碟所需的時間。此外，當完成建立2^32範圍內的質數表時，1秒內可檢測2^64範圍內任一數字是否為質數，幾乎等同擁有2^64（18,446,744,073,709,551,616）範圍的質數表。（本章節所謂大數字或大質數意指大於2^64的數字）
試用版應用程式“質數表生成器”下載處：
http://hanteye01.blog.fc2.com/blog-entry-23.html 。
以處理器12th Gen Intel(R) Core(TM)
i5-12400F 2.50 GHz以及1T容量的硬碟為例，以下是實際測試結果：
為了便於敘述，此處定義質數佔有率（Primes Occupancy Rate）函數POR(n,m)：從n值開始，前m個質數的佔用率。
|
質數起始值 |
POR(n, 10000)(以第10000個質數值估算) |
|
10^5 (100,003) |
8.9527% |
|
10^6 (1,000.003) |
7.2167% |
|
10^7 (10,000,019) |
6.2290% |
|
10^8 (100,000,007) |
5.4419% |
|
10^9 (1,000,000,007) |
4.8042% |
|
10^10 (10,000,000,019) |
4.3429% |
|
10^11 (100,000,000,003) |
3.9609% |
|
10^12 (1,000,000,000,039) |
3.6179% |
|
大數字起始值 |
所需硬碟容量 |
初始化資料結構檔所需時間 |
每秒檢測大數字數量 |
每秒搜尋的大質數數量 |
POR(n, 10000)(以第10000個質數值估算) |
|
10^20 |
3.572 GB x 2 |
00:01:52 |
1567 |
32 |
2.1595% |
|
10^21 |
10.69 GB x 2 |
00:06:12 |
595 |
12 |
2.0939% |
|
10^22 |
32.195 GB x 2 |
00:19:47 |
203 |
3 |
1.9851% |
|
10^23 |
97.14 GB x 2 |
01:03:27 |
101 |
1 |
1.8948% |
|
10^24 |
293.87 GB x 2 |
04:59:00 |
67 |
1 |
1.8625% |
當質數值達到10^22時，質數佔有率已經低於2%。很難想象500位數質數的佔有率該是多少?即使最快速的多項式時間確定型質數判定演算法，1天內究竟能找到多少個500位數質數？
此外，根據上述質數佔有率的相關數據，可發現一個規律性：當自然數n值夠大時（至少大於m），隨著m值越大，則POR(n,m)越接近POR(n2,m)的2倍。如果這個規律是完全正確的，則可推論出下列關係式：
《運行效能分析》
整個生成大質數表過程中，各個階段皆有其瓶頸。以下是分析結果：
|
步驟項目 |
主要功能 |
運行瓶頸 |
|
階段1：初始化資料結構檔案 |
1.
預先配置虛擬記憶體所需的檔案大小。 2.
依據生成大質數表起始值，逐步設置每個質因數的起始位置。 |
1.
完全取決於磁碟效能。 2. CPU佔85%，磁碟佔15%。（粗估值） |
|
階段2：生成大質數表 |
根據資料結構檔內容，逐步生成大質數，同時修改資料結構檔內容，並且逐漸增加資料結構檔的長度。由於存取資料結構檔的位置幾乎是隨機的，整個運算效率完全取決於磁碟效能。 如果能保持資料結構檔為連續無分段的情況，可大幅提升效能。 |
完全取決於磁碟效能，CPU的效能無法顯著提升效率。 |
留言
發佈留言