Primality Test

Primality Test

「質數測試」。測試一個數字是否為質數。

質數測試屬於P問題。然而以下介紹的都是指數時間演算法,甚至是不保證結果正確的演算法。

多項式時間、保證結果正確的演算法,請逕行上網搜尋「AKS Primality Test」;但由於方法繁雜,實務上沒人使用。

質數測試,可以直接以篩法建立質數表,再來判斷質數。然而建立質數表需要大量記憶體,因此又發明了其他方法。

Primality Test: Divisibility Primality Test

Divisibility Primality Test

整除性測試法即是試除法。

質數無法以乘法分解,質數沒有任何因數(除了1和自身)。窮舉所有可能的因數,一一試除。時間複雜度O(n)。

因數成雙成對(平方根跟自己一對)。檢查了小於等於平方根的因數,就不必再檢查大於平方根的因數。時間複雜度O(sqrt(n))。

2n+1 Method

因數可以分解成質數。取質數試除,效果等同於取因數試除。反過來說,大可不必取合數試除。

2的更大倍數是合數,不必試除,節省一半時間。

6n±1 Method

2和3的最小公倍數是6。所有數字分為6n、6n+1、6n+2、6n+2、6n+3、6n+4、6n+5六種(n是倍率)。除以六餘零的數字為6n,除以六餘一的數字為6n+1,以此類推。

6n、6n+2、6n+3、6n+4是2或3的倍數,不是質數。只需要拿6n+1和6n+5試除。

從6n+1,跳到6n+5,跳到下一個6n+1,再跳到6n+5,差值是4 2 4 2 4 2...,不斷跳四跳二。

6n+5改寫成6n-1,意義不變。6n+1與6n-1簡寫成6n±1。

2和3預先處理。從5開始試除,下一個是5+2,再下一個是5+2+4,再下一個是5+2+4+2,以此類推。

實作程式碼,用加法加二加四,不用乘法計算6n-1、6n+1,效率較佳。

Wheel Factorization

6n±1法進行推廣,使用頭幾個質數。

取2 3 5的最小公倍數30,建立30n+k法。取2 3 5 7的最小公倍數210,建立210n+k法。取所有質數,等同篩法。

取越多質數,即可剔除越多合數,事前節省時間。

然而,取越多質數,也會建立越多規則,事中浪費時間。

這是取捨兩難的情況,兩者之間必須取得平衡。一般來說,大家只取2 3或者2 3 5。

更有甚者,當你還在煩惱要取多少質數,普通的試除法早就已經執行完畢了。更有甚者,當你還在煩惱要如何實作試除法,網路上已經有現成的質數表了。所謂的取得平衡,端看你的視野高度。

Primality Test: Miller-Rabin Primality Test

Square Root Primality Test

餘數系統的平方根有個特性:

以質數n為模數,0到n-1之間,只有1與n-1,平方之後等於1。
以質數n為模數,1的平方根一定是1與n-1(即±1)。

平方根測試法,運用了此特性:

以n為模數,當2到n-2的每一個數字,平方之後皆不等於1,就推定n是質數。

此演算法的結果不一定正確。通過測試的數字,可能是合數或質數;無法通過測試的數字,一定是合數。

以合數n為模數,0到n-1之間,除了1與n-1,還可能有其它數字,平方之後等於1。
以合數n為模數,1的平方根可能不只是1與n-1(不只±1)。

有些合數會被誤判成質數,例如22會被誤判成質數,2²、…、20²模22之後皆不等於1。

Fermat Primality Test

費瑪小定理:

若n是質數,a小於n,則aⁿ ≡ a (mod n)。
若n是質數,a小於n,a不是零,則aⁿ⁻¹ ≡ 1 (mod n)。

費瑪質數測試法,運用了費瑪小定理:

n是質數,費瑪小定理一定成立:aⁿ⁻¹ % n = 1一定成立。
n是合數,費瑪小定理可能成立:aⁿ⁻¹ % n = 1可能成立。

當aⁿ⁻¹ % n = 1成立,就推定n是質數。

此演算法的結果不一定正確。通過測試的數字,可能是合數或質數;無法通過測試的數字,一定是合數。

使用各式各樣的a來實施測試,那麼判定結果就更準確。

Miller-Rabin Primality Test

結合Square Root Primality Test與Fermat Primality Test。

此演算法的結果不一定正確。通過測試的數字,可能是合數或質數;無法通過測試的數字,一定是合數。

一、選定一個底數a,大於1、小於n,用來進行費瑪測試。
  待測數字n,會是費瑪測試的模數。
二、令 n-1 = u × 2ᵗ,其中t盡量大(u為奇數)。
三、觀察a^u這個數字:
  若等於±1,
  則表示步驟四,所有數字都是1,永不出現平方根測試的反例。
  也表示步驟五,最終將通過費瑪測試。
  推定n是質數。
四、依序觀察a^(u × 2¹)、a^(u × 2²)、…、a^(u × 2ᵗ⁻¹)這些數字,
  每個數字正好是前一個數字的平方:
 甲、一旦發現有個數字的平方等於+1,
   則表示無法通過平方根測試。
   (但是步驟五,最終將通過費瑪測試。)
   確定n是合數。
 乙、一旦發現有個數字的平方等於-1,
   則表示接下來的數字都是1,永不出現平方根測試的反例。
   也表示步驟五,最終將通過費瑪測試。
   推定n是質數。
五、觀察a^(u × 2ᵗ)這個數字:
 甲、若等於+1,表示通過費瑪測試,推定n是質數。
 乙、若不等於+1,表示無法通過費瑪測試,確定n是合數。
 回、也就是說,a^(u × 2ᵗ⁻¹)必須等於±1,平方之後才會等於+1。
   步驟五才能通過費瑪測試。
 回、由於步驟四已經檢查過a^(u × 2ᵗ⁻¹)是否為±1,
   所以步驟五可以省略,直接確定n是合數。

步驟二:不斷分解因數2,最多log(n)步,需時O(log(n))。

步驟三:餘數次方,Divide-and-Conquer Method,需時O(log(n))。

步驟四:根據步驟二,最多log(n)-1個數字,需時O(log(n))。

總時間複雜度O(log(n))。

Strong Probable-prime Base(SPRP)

以過濾合數的角度來看,多取幾個相異的底數a實施測試,判定結果就更準確。

事實上已經有熱心人士,找出特定的底數組合,仔細檢查了各種數字的判定結果,保證判定結果正確。例如底數組合{2, 7, 61}就能正確判斷2³²以下的數字是不是質數。

http://primes.utm.edu/prove/prove2_3.html

http://miller-rabin.appspot.com/

UVa 10956

Primality Test: Baillie-PSW Primality Test

Lucas Primality Test

https://en.wikipedia.org/wiki/Lucas_primality_test

Baillie-PSW Primality Test

https://faculty.lynchburg.edu/~nicely/misc/bpsw.html

https://eprint.iacr.org/2018/749.pdf

目前仍然可以正確判斷質數,尚未找到反例。

Integer Factorization

Factor(Divisor)

一個數字的「因數」,乘法分解而得的數字。通常有許多個。

0: 0 1 2 ...   5: 1 5       10: 1 2 5 10
1: 1           6: 1 2 3 6   11: 1 11
2: 1 2         7: 1 7       12: 1 2 3 4 6 12
3: 1 3         8: 1 2 4 8   13: 1 13
4: 1 2 4       9: 1 3 9     14: 1 2 7 14

Prime Factor

一個數字的「質因數」。又是因數、又是質數。

0: 2 3 5 ...   5: 5         10: 2 5
1: 1           6: 2 3       11: 11
2: 2           7: 7         12: 2 3
3: 3           8: 2         13: 13
4: 2           9: 3         14: 2 7

Integer Factorization

「整數分解」。一個正整數分解成質因數的連乘積。

               5 = 5        10 = 2 × 5
1 = 1          6 = 2 × 3    11 = 11
2 = 2          7 = 7        12 = 2² × 3
3 = 3          8 = 2³       13 = 13
4 = 2²         9 = 3²       14 = 2 × 7

正式名稱稱作「整數分解」,著眼於分解對象。中學課本譯作「質因數分解」,著眼於分解結果。

質因數分解屬於NP問題,但是目前還不確定它究竟是P問題或是NP-complete問題,相當特別。

Fundamental Theorem of Arithmetic

「算術基本定理」。凡是正整數,都可以藉由質因數分解,成為一個獨一無二的式子。

n = 2n₁ × 3n₂ × 5n₃ × 7n₄ × 11n₅ × …

換句話說,不同的n對應不同的(n₁, n₂, …)。反方向亦然。

n  ⟷  (n₁, n₂, n₃, n₄, n₅, …)

根據算術基本定理,凡是牽扯到乘法、除法、因數、倍數的數學運算,可以變成比較簡單的運算。

分解前    	| 分解後
n         	| (n₁, n₂, …)
m         	| (m₁, m₂, …)
--------------------------------------------------
乘除法    	| 加減法
n × m     	| (n₁ + m₁, n₂ + m₂, …)
n ÷ m     	| (n₁ - m₁, n₂ - m₂, …)
          	|
整除      	| 大於等於
n % m = 0 	| (n₁, n₂, …) ≥ (m₁, m₂, …)
m | n     	| n₁≥m₁ and n₂≥m₂ and …
          	|
最大公因數	| 最小值
gcd(n, m) 	| (min(n₁, m₁), min(n₂, m₂), …)
          	|
最小公倍數	| 最大值
lcm(n, m) 	| (max(n₁, m₁), max(n₂, m₂), …)
          	|
互質      	| 對應項必須有零
n ⊥ m     	| min(n₁, m₁) = 0 and min(n₂, m₂) = 0 and …
          	| n₁×m₁ = 0 and n₂×m₂ = 0 and …

算術基本定理闡述了另一種世界觀,把數字看作是質數的結合。質數的英文prime有著「原始就有」的意思,便是指質數是所有數字的根本。

UVa 11889

Integer Factorization: Trial Division Algorithm

Trial Division Algorithm

把所有可能的因數拿來試除。由小到大試除,即時消滅因數,因此每次整除的是質因數。時間複雜度O(n)。

n不是質數的情況,遠小於O(n)。精確的寫法是O(max(pₖ) + sum(eₖ))。pₖ和eₖ來自質因數分解n = p₁e₁ × p₂e₂ × ...。

平均時間複雜度,似乎可以精確用n表示,但是我沒有研究。讀者可以逕行上網搜尋「最大質因數分布」、「質因數數量分布」。

因數成雙成對(平方根跟自己一對),窮舉一半足矣。乘法分解的一半:平方根。時間複雜度降為O(sqrt(n))。

精確的寫法是O(sqrt(max(pₖ)) + sum(eₖ))。

預先處理質因數2,然後每次跳二,節省一半時間。

Wheel Factorization可以節省更多時間。

UVa 516 583 10179 10290 10329 10392 10622 10780 10791 10879

Integer Factorization: Pollard's ρ Algorithm

亂數產生器

此演算法可以找到n的其中一個因數。

使用一個簡單的亂數產生器f(aᵢ₊₁) = aᵢ² + c (mod n),嘗試各種a₀與c,製造所有可能的因數,一一試除即可。

以此亂數產生器公式,依序枚舉a₀ a₁ a₂ …,模n的情況下,最終必定循環。圖示習慣畫成一個ρ的形狀,演算法因此而得名。

希臘字母ρ,唸做/ro/,可以寫做rho。

運用最大公因數找到因數

因為亂數產生器製造的數字a,a恰是n的因數的機會較小,而a與n有共同因數的機會較大,所以改用d = gcd(a, n)來找到n的因數d。最大公因數有著極快的演算法,對執行速度影響不大。

偵測循環

亂數產生器最終必會出現重覆數字,產生循環。一旦遇到循環,立刻結束枚舉,不再進行重覆運算。

另外,此演算法改用abs(ax - ay),用數字的差值製造所有可能的因數。我不知道如此做的原因。

原論文採用「Floyd's Algorithm」偵測循環。

亦可採用「Brent's Algorithm」偵測循環,效率較佳。

讀者可以思考一下這些問題

一、為何亂數產生器不採用f(aᵢ₊₁) = aᵢ + c (mod n)這條更加簡單的式子?

二、為何a₀至少要是+2?(經過實測,a₀最好是+2。)

三、為何c最好不是0和-2?試試看將亂數產生器公式,代入到x和y之中,計算一下x-y,然後計算一下gcd(abs(x-y), n)。

四、為何x等於y的時候,就要馬上結束迴圈呢?

五、如果略去abs(),改成gcd(x-y, n),對結果有影響嗎?

質因數分解

似乎只要嘗試各種c,就一定可以窮舉出所有可能的因數。我不知道原因。

UVa 11476 11466

Integer Factorization: Quadratic Sieve Algorithm

Fermat's Factorization Method

一個數字n,分解成兩個數x y的乘積。

運用平方差公式,令n = a² - b² = (a+b) (a-b) = x y。尋找剛好相差n的兩個平方數a²與b²,就能分解n。

實作程式碼時,窮舉整數a,然後推導出b = sqrt(a² - n),判斷b是不是整數。

另一種更直覺的解讀方式:

窮舉a²,計算a² - n,判斷是不是平方數。如果是平方數,便得到了剛好相差n的兩個平方數a²與b²。

a²由小到大依序窮舉。由於a²必須大於等於n,才能相減得到b²,所以由等於n、略大於n的平方數開始窮舉。

n = 52   a₀ = ceil(sqrt(52)) = 8

 a |  a² | a² - n | ² ? |
---|-----|--------|-----|
 8 |  64 |   12   |  X  |
 9 |  81 |   29   |  X  |
10 | 100 |   48   |  X  |
11 | 121 |   69   |  X  |
12 | 144 |   92   |  X  |
13 | 169 |  117   |  X  |
14 | 196 |  144   |  O  |

52 = 196 - 144 = 14² - 12² = (14+12) (14-12) = 26 × 2
n  = a²  - b²  = a²  - b²  = (a+b) (a-b)     = x y

Fermat's Factorization Method推廣至餘數系統

本來是尋找剛好相差n的兩個平方數,現在是尋找相差kn的兩個平方數,k是任意整數倍率。放寬限制,機會更多。

令kn = a² - b² = (a+b) (a-b)。想要分解n,可以計算d = gcd(n, a+b),從a+b當中擷取n的因數。或者a-b亦可。

運氣不好時,d可能等於1或n,沒有達到分解效果;此時再嘗試其他平方數組合即可。

回到一開始。n變kn,可以想成是模n。原問題變成:尋找兩個平方數,模n同餘。

n = a² - b²           差n
n ≡ a² - b² (mod n)   推廣成餘數。模n,就變成差kn。
0 ≡ a² - b² (mod n)   顯然
a² ≡ b² (mod n)       移項。式子好看又好用。

當d = gcd(n, a+b) = 1 or n,此時a ≡ ±b (mod n)。

Dixon's Factorization Method

窮舉a²,可是a² - n一直不是平方數,怎麼辦?

各個a² - n相乘,湊出平方數,如此便得到了剛好相差kn的兩個平方數a²與b²。

n = 26   a₀ = ceil(sqrt(26)) = 6

 a |  a² | a² - n | ² ? |
---|-----|--------|-----|
 6 |  36 |   10   |  X  | ✓
 7 |  49 |   23   |  X  | ✓
 8 |  64 |   38   |  X  |
 9 |  81 |   55   |  X  |
 : |   : |    :   |  :  |
16 | 256 |  230   |  X  | ✓
 : |   : |    :   |  :  |

{ 6²  ≡ 10  (mod 26)
{ 7²  ≡ 23  (mod 26)   ===>   (6×7×16)² ≡ 230² (mod 26)
{ 16² ≡ 230 (mod 26)

等號左邊相乘 6² × 7² × 16² = (6×7×16)² 平方數相乘仍是平方數
等號右邊相乘 10×23×230 = 230²          人工湊得平方數
得到一組相差 kn 的兩個平方數 a² 與 b²!

為了方便湊出平方數,運用算術基本定理:a² - n實施質因數分解,取次方值,形成向量。平方數的次方值皆為偶數。

為了快速湊出平方數,運用線性代數:向量模2,挑其中幾個向量,XOR等於零。線性組合等於零。解一次方程組。高斯消去法。

窮舉所有組合,簡化成高斯消去法,時間複雜度從指數變成線性,大幅改進!此為精髓!

 a |  a² | a² - n | ² ? | factor   | vector (mod 2)  |
---|-----|--------|-----|----------|-----------------|
 6 |  36 |   10   |  X  | 2 5      | [1 0 1 0 0 ...] | ✓
 7 |  49 |   23   |  X  | 23       | [0 0 0 0 0 ...] | ✓
 8 |  64 |   38   |  X  | 2 19     | [1 0 0 0 0 ...] | 
 9 |  81 |   55   |  X  | 5 11     | [0 0 1 0 1 ...] | 
 : |   : |    :   |     | :        |      :          | 
16 | 256 |  230   |  X  | 2 5 23   | [1 0 1 0 0 ...] | ✓
 : |   : |    :   |     | :        |      :          | 
24 | 576 |  550   |  X  | 2 5 5 11 | [1 0 0 0 1 ...] | 
 : |   : |    :   |     | :        |      :          | 

                               [ x₆  ]
      [ 1 0 1 0 .. 1 .. 1 .. ] [ x₇  ]
      [ 0 0 0 0 .. 0 .. 0 .. ] [ x₈  ]
      [ 1 0 0 1 .. 1 .. 0 .. ] [ x₉  ]
solve [ 0 0 0 0 .. 0 .. 0 .. ] [  :  ] = 0
      [ 0 0 0 1 .. 0 .. 1 .. ] [ x₁₆ ]
      [ : : : :    :    :    ] [  :  ]
      [ : : : :    :    :    ] [ x₂₄ ]
                               [  :  ]

矩陣通常很大,於是採用遞推法:a² - n實施質因數分解,逐步增加質因數;每多一個質因數,就解一次方程組。

每回合只拿已經徹底分解的a² - n,來解一次方程組;也就是除到剩下1的。

用數學術語來說:逐步增加B,每次只拿B-smooth的a² - n,來解一次方程組。

n = 26
   |     |        | round 1      | round 2      | round 3      
 a |  a² | a² - n |  ÷2   factor |  ÷3   factor |  ÷5   factor 
---|-----|--------|--------------|--------------|--------------
 6 |  36 |   10   |   5 | 2      |   5 | 2      |   1 | 2 5    
 7 |  49 |   23   |  23 |        |  23 |        |  23 |        
 8 |  64 |   38   |  19 | 2      |  19 | 2      |  19 | 2      
 9 |  81 |   55   |  55 | 5      |  55 | 5      |  11 | 5      
10 | 100 |  100   |  37 | 2      |  37 | 2      |  37 | 2      
 : |   : |    :   |   : | :      |   : | :      |   : | :      
24 | 576 |  576   | 225 | 2      |  25 | 2      |   1 | 2 5 5  
25 | 625 |  625   | 599 |        | 599 |        | 599 |        

最後介紹一個加速技巧:剩下的數字,如果不為1,但是數字一樣,此時可以融合得到一組新的B-smooth。

有時候剩下的數字很大、仍需數回合才能除盡。運用此技巧,得以提早納入一些a² - n,提早找到正確答案。

n = 26
   |     |        || round 5      |
 a |  a² | a² - n || ÷11   factor |
---|-----|--------||--------------|
 : |   : |     :  ||   : | :      |
14 | 196 |   170  ||  17 | 2 5    |
 : |   : |     :  ||   : | :      |
20 | 400 |   374  ||  17 | 2 11   |

剛好都剩17,可以相乘,去掉17,得到一組新的11-smooth。
甚至可以考慮消除原本那兩組。

   |     |        || round 5      |
 a |  a² | a² - n || ÷11   factor |
---|-----|--------||--------------|
280| *** | ****** ||   1 | 2 5 11 |

Quadratic Sieve Algorithm

整除質因數2的間隔似乎是2。似乎是篩法。仔細推導一下:

已知某個 a² - n 整除質因數 p。現在令間隔為p:
  (a+p)² - n
= a² + 2ap + p² - n
= (a² - n) + 2ap + p²
= (a² - n) + p (2a+p)   仍整除質因數p!
             ~~~~~~~~

在0到p-1之間,找到所有的a,滿足a² - n整除質因數p,作為篩法的起點;跳躍間隔皆是p。如此一來便無一遺漏。

a² - n ≡ 0 (mod p)
a² ≡ n (mod p)      已知n p,解a,找出所有的a。

餘數平方根演算法,請參考「Cipolla's Algorithm」、「Tonelli-Shanks Algorithm」。

經過簡單的觀察和推導,一一試除改為篩法,節省不少時間。此即當今世上第二快的質因數分解演算法。