3.3 大數質因數分解#

概念說明#

把整數 n 分解成質因數的乘積，聽起來是中小學程度的工作，但當 n 有幾百個位數時，它變成了世界級的難題：目前沒有已知的多項式時間經典演算法。這個「分解困難」假設撐起了整個 RSA 密碼系統的公鑰——公開的是合成數 n = p·q，私鑰則藏著質數 p 與 q。只要分解不出 n，第三人就算拿到公鑰也解不開密文。

有趣的不對稱在於：測質數很快（2.1 節米勒–拉賓，多項式時間），分解很慢。所以我們知道一個 2048 位元的數幾乎肯定是質數（用測試），卻可能完全無法分解一個 2048 位元的合成數。本節按難度遞增介紹三種方法：

試除法：vérifier質因數繼續試到 √n，O(√n)。n 到 10¹⁶ 以內可用。
Pollard 的 ρ 演算法：用隨機行走產生「模 n 下的偽隨機數列」，透過差值求 gcd 撈出因數，期望 O(n^(1/4))。可以把 n = p**2 之類的數字分解到幾十位以上。
（概念介紹）二次篩法與一般數體篩法（GNFS）：目前最快的亞指數演算法，GNFS 複雜度約 exp(1.9·(ln n)^(1/3)·(ln ln n)^(2/3))，RSA-250（829 位元）正是用 GNFS 分解出來的，動用上千核年。

Pollard ρ 的核心是費馬偽隨機函數 f(x) = x² + c (mod n)：從種子 x₀ 出發的迭代 x, f(x), f(f(x)), … 在模 n 的世界裡隨機漫步。若某條鏈軌跡的差值 δ = xᵢ − xⱼ 與 n 有公因數 d = gcd(δ, n)，且 1 < d < n，這個 d 就是 n 的非平凡因數。因為差值不屬於質因數 p 的集合。由於大數模 n 不可巡遊，我們用 Floyd 循環偵測或 Brent 優化，以「單位時間兩步」的方式抓出軌跡的接近值。

演算法步驟#

Pollard ρ（Brent 版）： 1. 若 n 為偶數，回傳 2；先用米勒–拉賓檢查 n 是否本身為質數，是則回 n。 2. 隨機選 c ∈ [1, n−1)、起點 x。令 y = x，d = 1。 3. 反複執行：（a） x = f(x)，y = f(f(y))——y 每輪前進兩步；（b） d = gcd(|x − y|, n)；（c）若 d = n（找到循環但無因數），換 c 重來；若 1 < d < n，回傳因數 d。 4. 對得到的 d 與 n/d 遞迴分解。

試除法： 從 2 開始試除，除盡就收，直到試商大於 √n。每除掉一個質因數就縮小 n，加速後續。

Python 程式碼範例#

import math, random
from sympy import factorint, isprime

def factor_small(n):
    """試除法：n <= 10^16"""
    factors = {}
    d = 2
    while d * d <= n:
        while n % d == 0:
            factors[d] = factors.get(d, 0) + 1
            n //= d
        d += 1 if d == 2 else 2        # 質數之外跳過偶數
    if n > 1:
        factors[n] = factors.get(n, 0) + 1
    return factors

def pollard_rho(n):
    """回傳 n 的一個非平凡因數（n 為合成奇數）"""
    if n % 2 == 0:
        return 2
    while True:
        x = random.randrange(2, n)
        c = random.randrange(1, n)
        y = x
        d = 1
        while d == 1:
            x = (x * x + c) % n
            y = (y * y + c) % n
            y = (y * y + c) % n
            d = math.gcd(abs(x - y), n)
        if d != n:
            return d                    # 否則換 c 重來

def factorize(n):
    """遞迴分解：回傳 {質因數: 次數}"""
    if n == 1:
        return {}
    if isprime(n):                      # 米勒–拉賓（sympy 提供）
        return {n: 1}
    d = pollard_rho(n)
    f1 = factorize(d)
    f2 = factorize(n // d)
    for k, v in f2.items():
        f1[k] = f1.get(k, 0) + v
    return f1

# 試除法示範
print(factor_small(2**10 * 3**5 * 257))
print(factor_small(987654321987))

# Pollard ρ 示範（質因數超過 sqrt 可處理範圍也能解）
n = 10**12 + 39                 # 合成數測試
print(factorize(n), "乘回去 =", math.prod(k**v for k, v in factorize(n).items()))
print(factorize(2**67 - 1))     # 梅森數：分解出 193707721 × 761838257287

# 與 sympy 比對
print(factorint(2**67 - 1) == factorize(2**67 - 1))

# 單向性示範：2048 位元的 n = p*q，乘法毫秒、分解則天文時間
p = 2**1024 + 1     # 這裡僅示例；真實 RSA 需安全質數
q = 3**650 + 2
n = p * q
print(len(str(n)))  # n 的位數，乘法瞬間完成；分解 n 則遠超宇宙年齡

複雜度分析#

試除法 O(√n)：1024 位元的 n 需要試約 2^512 次——絕望。Pollard ρ 的期望時間 O(n^(1/4))（來自 ρ 形狀循環的生日悖論式碰撞論證），實務上可分解約 50~60 位元內的數。GNFS 亞指數級：RSA-240（795 位元）分解動用約 900 核年；RSA-250 花了約 2700 核年；這些數字隨位數爆增但非爆炸——亞指數增長比多項式快、比指數慢，處在密碼學的最痛區間。Pollard ρ 是隨機化演算法：同個 n 每次跑可能回不同的因數（因為 c 與起點隨機），但「期望時間有界」的保證是不變的。實作上 Brent 的迴圈偵測是先把 gcd 的計算「湊批」再算（一次 gcd 處理多步差值的乘積），能省掉大部分 gcd 開銷；本節程式為了可讀性逐步算 gcd，速度稍慢但邏輯等價。

練習#

把 Pollard ρ 改成 Brent 優化版：累積 100 步的差值乘積 mod n 後一次 gcd，並量測分解 2**67−1 的加速比。
用 factorize 分解 n = 1012 + 39、1012 + 457、2**64 + 1，並與 sympy.factorint 比對；記錄各自時間，觀察 Pollard ρ 的期待效率。
「測質快、分解慢」實驗：產生兩個約 30 位元的隨機安全質數 p、q（用 isprime 搜尋），計算 n = p·q 並列印位數；分別量測 (a) isprime(n)、(b) factorize(n) 的時間，解釋為何 (a) 秒殺而 (b) 可能要幾秒到幾分鐘。