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

trial division algorithm

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

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

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

void factorize(int n) { for (int p=2; p<=n; ++p) // n不斷變小 while (n % p == 0) { cout << p; // 印出質因數 n /= p; // 消滅質因數 } }

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

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

void factorize(int n) { for (int p=2; p*p<=n; ++p) // n不斷變小 while (n % p == 0) { cout << p; // 印出質因數 n /= p; // 消滅質因數 } if (n > 1) cout << n; // 印出質因數 }

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

wheel factorization可以節省更多時間。

void factorize(int n) { while (n % 2 == 0) { cout << 2; n /= 2; } for (int p=3; p*p<=n; p+=2) // 跳二 while (n % p == 0) { cout << p; n /= p; } if (n > 1) cout << n; }

UVa 516 583 10179 10290 10329 10392 10622 10780 10791 10879

亂數產生器

此演算法可以找到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」偵測循環。

#define f(x) ((x * x + c) % n) // 亂數產生器 int pollard_rho(int n) { int a0 = rand() % n; // 最好是+2 int c = rand() % n; // 最好不是0和-2 int x = a0, y = a0; while (true) { // Floyd's algorithm。 // 一旦出現循環，則回傳原本數字n。 x = f(x); y = f(f(y)); if (x == y) return n; int d = gcd(abs(x - y), n); if (d > 1 && d < n) return d; } return -1; // 不會執行到這行 } #define f(x) ((x * x + c) % n) // 亂數產生器 int pollard_rho(int n) { int a0 = rand() % n; // 最好是+2 int c = rand() % n; // 最好不是0和-2 int x = a0, y = a0; while (true) { x = f(x); y = f(f(y)); // if (x == y) return n; // 因為出現循環之時，x-y等於0，gcd(0,n)等於n。 // 所以亦得用d==n判斷是否出現循環，減少一行程式碼。 int d = gcd(abs(x - y), n); if (d > 1 && d <= n) return d; } return -1; // 不會執行到這行 }

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

int pollard_rho(int n) { int a0 = rand() % n; // 最好是+2 int c = rand() % n; // 最好不是0和-2 int x = a0, y = a0; int i = 1, k = 2; while (true) { x = (x * x + c) % n; // if (x == y) return n; int d = gcd(abs(x - y), n); if (d > 1 && d <= n) return d; if (++i == k) y = x, k <<= 1; } return -1; // 不會執行到這行 }

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

一、為何亂數產生器不採用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，就一定可以窮舉出所有可能的因數。我不知道原因。

vector<int> factor; // 存放n的質因數 int pollard_rho(int n, int c) { int x = 2, y = 2, i = 1, k = 2; while (true) { x = (mul(x, x, n) + c) % n; int d = gcd(abs(x - y), n); if (d > 1 && d <= n) return d; if (++i == k) y = x, k <<= 1; } return -1; } void factorize(int n) { if (n == 1) return; // 預先判斷質數。 // 可以使用Miller–Rabin primality test with SPRP。 if (isprime(n)) { factor.push_back(n); return; } // 此處只有合數才能實施Pollard's ρ algorithm， // 否則永遠找不到1與n以外的因數，陷入無窮迴圈。 // 窮舉各種c，直到發現n的因數d，一定能找到。 int d = n; for (int c=0; d == n; ++c) d = pollard_rho(n, c); factorize(d); factorize(n / d); }

UVa 11476 11466

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是不是整數。

void fermat(int n) { for (int a = ceil(sqrt(n)); a<n; ++a) { int bb = a * a - n; int b = sqrt(bb); if (b * b == bb) // 判斷b是不是整數 { int x = a + b, y = a - b; cout << "n =" << x << "*" << y; // 遞迴分解x和y，無法再分解的時候就找到質因數了。 fermat(x); fermat(y); break; } } }

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

窮舉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

請見專著《Algorithmic Number Theory》。

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

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

tuple(int, set<int>, int) refine(int a, int b) { int d = gcd(a, b); if (d == 1) return {a, {}, b}; auto [l1, s1, r1] = refine(a/d, d); auto [l2, s2, r2] = refine(r1, b/d); s1.insert(s2.begin(), s2.end()); // set union s1.insert(l2); return {l1, s1, r2}; // l|a, r|b, s|d, l⊥b, r⊥a } set<int> refine(int a[], int n) { if (n == 1) return {a[0]}; if (n == 2) { auto [l,s,r] = refine(a[0], a[1]); s.insert(l); s.insert(r); s.erase(1); return s; } set<int> s; set<int> t = refine(a, n-1); int r0 = a[n-1]; for (auto ti: t) { auto [l, s, r] = refine(ti, r0); s.insert(s); s.insert(l); r0 = r; } s.insert(r0); s.erase(1); return s; }