分解質因數

引入

給定一個正整數 \(N \in \mathbf{N}_{+}\)，試快速找到它的一個非平凡因數。

考慮樸素算法，因數是成對分佈的，\(N\) 的所有因數可以被分成兩塊，即 \([2, \sqrt N]\) 和 \([\sqrt N+1,N)\)。只需要把 \([2, \sqrt N]\) 裏的數遍歷一遍，再根據除法就可以找出至少兩個因數了。這個方法的時間複雜度為 \(O(\sqrt N)\)。

當 \(N\ge10^{18}\) 時，這個算法的運行時間我們是無法接受的，希望有更優秀的算法。一種想法是通過隨機的方法，猜測一個數是不是 \(N\) 的因數，如果運氣好可以在 \(O(1)\) 的時間複雜度下求解答案，但是對於 \(N\ge10^{18}\) 的數據，成功猜測的概率是 \(\frac{1}{10^{18}}\), 期望猜測的次數是 \(10^{18}\)。如果是在 \([2,\sqrt N]\) 裏進行猜測，成功率會大一些。我們希望有方法來優化猜測。

樸素算法

最簡單的算法即為從 \([2, \sqrt N]\) 進行遍歷。

C++Python

vector<int> breakdown(int N) {
  vector<int> result;
  for (int i = 2; i * i <= N; i++) {
    if (N % i == 0) {  // 如果 i 能夠整除 N，説明 i 為 N 的一個質因子。
      while (N % i == 0) N /= i;
      result.push_back(i);
    }
  }
  if (N != 1) {  // 説明再經過操作之後 N 留下了一個素數
    result.push_back(N);
  }
  return result;
}

def breakdown(N):
    result = []
    for i in range(2, int(sqrt(N)) + 1):
        if N % i == 0: # 如果 i 能夠整除 N，説明 i 為 N 的一個質因子。
            while N % i == 0:
                N //= i
            result.append(i)
    if N != 1: # 説明再經過操作之後 N 留下了一個素數
        result.append(N)
    return result

我們能夠證明 result 中的所有元素均為 N 的素因數。

證明 result 中均為 \(N\) 的素因數

首先證明元素均為 \(N\) 的素因數：因為當且僅當 N % i == 0 滿足時，result 發生變化：儲存 \(i\)，説明此時 \(i\) 能整除 \(\frac{N}{A}\)，説明了存在一個數 \(p\) 使得 \(pi=\frac{N}{A}\)，即 \(piA = N\)（其中，\(A\) 為 \(N\) 自身發生變化後遇到 \(i\) 時所除的數。我們注意到 result 若在 push \(i\) 之前就已經有數了，為 \(R_1,\,R_2,\,\ldots,\,R_n\)，那麼有 N \(=\frac{N}{R_1^{q_1}\cdot R_2^{q_2}\cdot \cdots \cdot R_n^{q_n}}\)，被除的乘積即為 \(A\)）。所以 \(i\) 為 \(N\) 的因子。

其次證明 result 中均為素數。我們假設存在一個在 result 中的合數 \(K\)，並根據整數基本定理，分解為一個素數序列 \(K = K_1^{e_1}\cdot K_2^{e_2}\cdot\cdots\cdot K_3^{e_3}\)，而因為 \(K_1 < K\)，所以它一定會在 \(K\) 之前被遍歷到，並令 while(N % k1 == 0) N /= k1，即讓 N 沒有了素因子 \(K_1\)，故遍歷到 \(K\) 時，N 和 \(K\) 已經沒有了整除關係了。

值得指出的是，如果開始已經打了一個素數表的話，時間複雜度將從 \(O(\sqrt N)\) 下降到 \(O(\sqrt{\frac N {\ln N}})\)。去篩法處查閲更多打表的信息。

例題：CF 1445C

Pollard Rho 算法

引入

而下面複雜度複雜度更低的 Pollard-Rho 算法是一種用於快速分解非平凡因數的算法（注意！非平凡因子不是素因子）。而在此之前需要先引入生日悖論。

生日悖論

不考慮出生年份（假設每年都是 365 天），問：一個房間中至少多少人，才能使其中兩個人生日相同的概率達到 \(50\%\)?

解：假設一年有 \(n\) 天，房間中有 \(k\) 人，用整數 \(1, 2,\dots, k\) 對這些人進行編號。假定每個人的生日均勻分佈於 \(n\) 天之中，且兩個人的生日相互獨立。

設 \(k\) 個人生日互不相同為事件 \(A\), 則事件 \(A\) 的概率為

\[ P(A)=\prod_{i=0}^{k-1}\frac{n-i}{n} \]

至少有兩個人生日相同的概率為 \(P(\overline A)=1-P(A)\)。根據題意可知 \(P(\overline A)\ge\frac{1}{2}\), 那麼就有

\[ P(A)=\prod_{i=0}^{k-1}\frac{n-i}{n} \le \frac{1}{2} \]

由不等式 \(1+x\le \mathrm{e}^x\) 可得

\[ P(A) \le \prod_{i=1}^{k-1}\exp\left({-\frac{i}{n}}\right)=\exp \left({-\frac{k(k-1)}{2n}}\right) \]

因此

\[ \exp\left({-\dfrac{k(k-1)}{2n}}\right) \le \frac{1}{2}\implies P(A) \le \frac{1}{2} \]

將 \(n=365\) 代入，解得 \(k\geq 23\)。所以一個房間中至少 \(23\) 人，使其中兩個人生日相同的概率達到 \(50\%\), 但這個數學事實十分反直覺，故稱之為一個悖論。

當 \(k>56\)，\(n=365\) 時，出現兩個人同一天生日的概率將大於 \(99\%\)¹。那麼在一年有 \(n\) 天的情況下，當房間中有 \(\frac{1}{2}(\sqrt{8n\ln 2+1}+1)\approx \sqrt{2n\ln 2}\) 個人時，至少有兩個人的生日相同的概率約為 \(50\%\)。

考慮一個問題，設置一個數據 \(n\)，在 \([1,1000]\) 裏隨機選取 \(i\) 個數（\(i=1\) 時就是它自己），使它們之間有兩個數的差值為 \(k\)。當 \(i=1\) 時成功的概率是 \(\frac{1}{1000}\)，當 \(i=2\) 時成功的概率是 \(\frac{1}{500}\)（考慮絕對值，\(k_2\) 可以取 \(k_1-k\) 或 \(k_1+k\)），隨着 \(i\) 的增大，這個概率也會增大最後趨向於 1。

利用最大公約數求出一個約數

\(n\) 和某個數的最大公約數一定是 \(n\) 的約數，即 \(\forall k \in\mathbf{N}_{+},\gcd(k,n) \mid n\)，只要選適當的 \(k\) 使得 \(1<\gcd(k,n)< n\)，就可以求得 \(n\) 的一個約數 \(\gcd(k,n)\)。滿足這樣條件的 \(k\) 不少，\(n\) 有若干個質因子，每個質因子及其大部分倍數都是可行的。

我們通過 \(f(x)=(x^2+c)\bmod n\) 來生成一個序列 \(\{x_i\}\)：隨機取一個 \(x_1\)，令 \(x_2=f(x_1),x_3=f(x_2),\dots,x_i=f(x_{i-1})\)。其中 \(c\) 是一個隨機選取的常數。

舉個例子，設 \(n=50,c=6,x_1=1\)，\(f(x)\) 生成的數據為

\[ 1, 7, 5, 31, 17, 45, 31, 17, 45, 31,\dots \]

可以發現數據在 \(x_4\) 以後都在 \(31,17,45\) 之間循環。如果將這些數如下圖一樣排列起來，會發現這個圖像酷似一個 \(\rho\)，算法也因此得名 rho。

Pollard-rho1

選擇 \(f(x)=(x^2+c)\bmod n\) 這個函數生成序列，是因為它有一個性質：\(\forall x \equiv y \pmod p, f(x) \equiv f(y) \pmod p\)，其中 \(p \mid n\)。

證明

若 \(x\equiv y \pmod p\)，則可以將它們表示為 \(x=k_1p+a\)，\(y=k_2p+a\)，滿足 \(k_1,k_2,a\in \mathbb{Z},a\in \left[0,p\right)\)。

\(f(x)=(x^2+c) \bmod n\)，因此 \(f(x)=x^2+c-kn\)，其中 \(k \in \mathbb{Z}\)。

\[ \begin{aligned} f(x) & = x^2+c-kn\\ & = (k_1p+a)^2+c-kn\\ & = k_1^2 p^2+2k_1pa+a^2+c-kn\\ & \equiv a^2+c \pmod p \end{aligned} \]

同理，\(f(y) \equiv a^2+c \pmod p\)，因此 \(f(x) \equiv f(y) \pmod p\)。

根據生日悖論，生成的序列中不同值的數量約為 \(O(\sqrt{n})\) 個。設 \(m\) 為 \(n\) 的最小非平凡因子，顯然有 \(m\leq \sqrt{n}\)。

將 \(\{x_i\}\) 中每一項對 \(m\) 取模，我們得到了一個新序列 \(\{y_i\}\)（當然也可以寫作 \(\{x_i \bmod m\}\)），並且根據生日悖論可以得知新序列中不同值的個數約為 \(O(\sqrt{m})\leq O(n^{\frac{1}{4}})\)。

因此，我們可以期望在 \(O(n^{\frac{1}{4}})\) 的時間內找到兩個位置 \(i,j\)，使得 \(x_i\neq x_j\wedge y_i=y_j\)，這意味着 \(n \nmid |x_i−x_j| \wedge m \mid |x_i−x_j|\)，我們可以通過 \(\gcd(n, |x_i-x_j|)\) 獲得 \(n\) 的一個非平凡因子。

同理，任何滿足 \(\forall x \equiv y \pmod p, f(x) \equiv f(y) \pmod p\) 的函數 \(f(x)\)（例如多項式函數）都可以用在此處，且我們希望它對大部分初始值儘可能快地進入循環、循環週期儘可能短（但對於一個合數與其任一質因子，循環週期長度不能相同，否則 Pollard Rho 算法將無法分離出這個因子）。

一般選擇 \(x^2+c\) 生成序列，原因是它滿足上述條件，且運行效率較高。同時，每次調用函數都將 \(c\) 設置為隨機常數可以保證分解失敗後重新分解的成功率。

性質

我們期望枚舉 \(O(\sqrt{m})\) 個 \(i\) 來分解出 \(n\) 的一個非平凡因子 \(\gcd(|x_i−x_j|,n)\)，因此。Pollard-rho 算法能夠在 \(O(\sqrt{m})\) 的期望時間複雜度內分解出 \(n\) 的一個非平凡因子，通過上面的分析可知 \(O(\sqrt{m})\leq O(n^{\frac{1}{4}})\)，那麼 Pollard-rho 算法的總時間複雜度為 \(O(n^{\frac{1}{4}})\)。

下面介紹兩種實現算法，兩種算法都可以在 \(O(\sqrt{m})\) 的時間複雜度內完成。

實現

Floyd 判環

假設兩個人在賽跑，A 的速度快，B 的速度慢，經過一定時間後，A 一定會和 B 相遇，且相遇時 A 跑過的總距離減去 B 跑過的總距離一定是圈長的 n 倍。

設 \(a=f(1),b=f(f(1))\)，每一次更新 \(a=f(a),b=f(f(b))\)，只要檢查在更新過程中 a、b 是否相等，如果相等了，那麼就出現了環。

我們每次令 \(d=\gcd(|x_i-x_j|,n)\)，判斷 d 是否滿足 \(1< d< n\)，若滿足則可直接返回 \(d\)。由於 \(x_i\) 是一個偽隨機數列，必定會形成環，在形成環時就不能再繼續操作了，直接返回 n 本身，並且在後續操作裏調整隨機常數 \(c\)，重新分解。

基於 Floyd 判環的 Pollard-Rho 算法

C++Python

ll Pollard_Rho(ll N) {
  ll c = rand() % (N - 1) + 1;
  ll t = f(0, c, N);
  ll r = f(f(0, c, N), c, N);
  while (t != r) {
    ll d = gcd(abs(t - r), N);
    if (d > 1) return d;
    t = f(t, c, N);
    r = f(f(r, c, N), c, N);
  }
  return N;
}

import random
def Pollard_Rho(N):
    c = random.randint(1, N - 1)
    t = f(0, c, N)
    r = f(f(0, c, N), c, N)
    while t != r:
        d = gcd(abs(t - r), N)
        if d > 1:
            return d
        t = f(t, c, N)
        r = f(f(r, c, N), c, N)
    return N

倍增優化

使用 \(\gcd\) 求解的時間複雜度為 \(O(\log N)\)，頻繁地調用會使算法運行地很慢，可以通過乘法累積來減少求 \(\gcd\) 的次數。如果 \(1< \gcd(a,b)\)，則有 \(1< \gcd(ac,b)\)，\(c\in \mathbf{N}_{+}\)，並且由歐幾里得算法有 \(1< \gcd(ac \bmod b,b)=\gcd(ac,b)\)。

我們每過一段時間將這些差值進行 \(\gcd\) 運算，設 \(s=\prod|x_0-x_j|\bmod n\)，如果某一時刻得到 \(s=0\) 那麼表示分解失敗，退出並返回 \(n\) 本身。每隔 \(2^k-1\) 個數，計算是否滿足 \(1< \gcd(s, n) < n\)。此處取 \(k=7\)，可以根據實際情況進行調節。

注意到在環上更容易分解出因數，我們可以先迭代一定的次數。

實現

C++Python

ll Pollard_Rho(ll x) {
  ll t = 0;
  ll c = rand() % (x - 1) + 1;
  // 加速算法，這一步可以省略
  for (int i = 1; i < 1145; ++i) t = f(t, c, x);
  ll s = t;
  int step = 0, goal = 1;
  ll val = 1;
  for (goal = 1;; goal <<= 1, s = t, val = 1) {
    for (step = 1; step <= goal; ++step) {
      t = f(t, c, x);
      val = val * abs(t - s) % x;
      // 如果 val 為 0，退出重新分解
      if (!val) return x;
      if (step % 127 == 0) {
        ll d = gcd(val, x);
        if (d > 1) return d;
      }
    }
    ll d = gcd(val, x);
    if (d > 1) return d;
  }
}

from random import randint
from math import gcd
def Pollard_Rho(x):
    c = randint(1, x-1)
    s = t = f(0, c, x)
    goal = val = 1
    while True:
        for step in range(1, goal+1):
            t = f(t, c, x)
            val = val * abs(t - s) % x
            if val == 0: 
                return x #如果 val 為 0，退出重新分解
            if step % 127 == 0:
                d = gcd(val, x)
                if d > 1:
                    return d
        d = gcd(val, x)
        if d > 1:
            return d
        s = t
        goal <<= 1
        val = 1

例題：P4718【模板】Pollard-Rho 算法

對於一個數 \(n\)，用 Miller Rabin 算法判斷是否為素數，如果是就可以直接返回了，否則用 Pollard-Rho 算法找一個因子 \(p\)，將 \(n\) 除去因子 \(p\)。再遞歸分解 \(n\) 和 \(p\)，用 Miller Rabin 判斷是否出現質因子，並用 max_factor 更新就可以求出最大質因子了。由於這個題目的數據過於龐大，用 Floyd 判環的方法是不夠的，這裏採用倍增優化的方法。

實現

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

int t;
long long max_factor, n;

long long gcd(long long a, long long b) {
  if (b == 0) return a;
  return gcd(b, a % b);
}

long long quick_pow(long long x, long long p, long long mod) {  // 快速幂
  long long ans = 1;
  while (p) {
    if (p & 1) ans = (__int128)ans * x % mod;
    x = (__int128)x * x % mod;
    p >>= 1;
  }
  return ans;
}

bool Miller_Rabin(long long p) {  // 判断素数
  if (p < 2) return 0;
  if (p == 2) return 1;
  if (p == 3) return 1;
  long long d = p - 1, r = 0;
  while (!(d & 1)) ++r, d >>= 1;  // 将d处理为奇数
  for (long long k = 0; k < 10; ++k) {
    long long a = rand() % (p - 2) + 2;
    long long x = quick_pow(a, d, p);
    if (x == 1 || x == p - 1) continue;
    for (int i = 0; i < r - 1; ++i) {
      x = (__int128)x * x % p;
      if (x == p - 1) break;
    }
    if (x != p - 1) return 0;
  }
  return 1;
}

long long Pollard_Rho(long long x) {
  long long s = 0, t = 0;
  long long c = (long long)rand() % (x - 1) + 1;
  int step = 0, goal = 1;
  long long val = 1;
  for (goal = 1;; goal *= 2, s = t, val = 1) {  // 倍增优化
    for (step = 1; step <= goal; ++step) {
      t = ((__int128)t * t + c) % x;
      val = (__int128)val * abs(t - s) % x;
      if ((step % 127) == 0) {
        long long d = gcd(val, x);
        if (d > 1) return d;
      }
    }
    long long d = gcd(val, x);
    if (d > 1) return d;
  }
}

void fac(long long x) {
  if (x <= max_factor || x < 2) return;
  if (Miller_Rabin(x)) {              // 如果x为质数
    max_factor = max(max_factor, x);  // 更新答案
    return;
  }
  long long p = x;
  while (p >= x) p = Pollard_Rho(x);  // 使用该算法
  while ((x % p) == 0) x /= p;
  fac(x), fac(p);  // 继续向下分解x和p
}

int main() {
  scanf("%d", &t);
  while (t--) {
    srand((unsigned)time(NULL));
    max_factor = 0;
    scanf("%lld", &n);
    fac(n);
    if (max_factor == n)  // 最大的质因数即自己
      printf("Prime\n");
    else
      printf("%lld\n", max_factor);
  }
  return 0;
}

參考資料與鏈接

https://en.wikipedia.org/wiki/Birthday_problem#Reverse_problem ↩

本页面最近更新：，更新历史
发现错误？想一起完善？在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用