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(aₓ - a₝),用數字的差值製造所有可能的因數。我不知道如此做的原因。

原論文採用「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。

餘數平方根演算法,請見本站文件「residue」。

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

GCD-free basis

GCD-free basis

請見專著《Algorithmic Number Theory》

質因數分解:一個數字分解成質數相乘。時間複雜度類別不明。

互質因數分解:每個數字分解成互質數字相乘,盡量使用最大公因數。時間複雜度是多項式時間。