大整数质因数分解

大整数质因数分解

Miller–Rabin 负责快速判断一个大整数是否为质数；
Pollard–Rho 负责从大合数中随机找到一个非平凡因子。

1 从试除法到随机分解

试除法分解一个整数 nn 需要枚举到 n\sqrt n。

当 nn 接近 101810^{18} 时，n\sqrt n 接近 10910^9，单次试除已经难以承受。此时需要把“逐个寻找因子”改成：

本章主线：安全模乘 → Miller–Rabin 素性测试 → Pollard–Rho 寻找非平凡因子 → 递归分解。

两部分算法承担不同任务：

  • Miller–Rabin 只能判断当前整数是否为质数；
  • Pollard–Rho 只能尝试找到当前合数的一个非平凡因子；
  • 递归过程负责继续拆分因子，直到所有部分都成为质数。

2 大整数模乘

2.1 为什么普通乘法会溢出

若 a,b,na,b,n 都接近 long long 上限，即使最终只需要

ab mod n,ab\bmod n,

中间的 abab 也可能远远超过 long long 的表示范围。

在支持 __int128_t 的 C++ 环境中，可以先将乘积提升到 128 位，再对模数取模。

1
2
3
4
5
6
7
8
9
10
ll mul( ll a, ll b, ll mod ) { 
    return (ll)( (i128)a * b % mod ); 
}

ll qpow( ll a, ll k, ll mod ) {
    ll res = 1 % mod; a %= mod;
    for( ; k; k >>= 1, a = mul( a, a, mod ) )
        if( k & 1 ) res = mul( res, a, mod );
    return res;
}

这里的 mul 是 Miller–Rabin 与 Pollard–Rho 的共同底座。

3 从 Fermat 测试到 Miller–Rabin

3.1 Fermat 测试的局限

费马小定理说明：若 nn 是质数，且 gcd⁡(a,n)=1\gcd(a,n)=1，则

an−1≡1(modn).a^{n-1}\equiv1\pmod n.

因此，若某个底数 aa 不满足这一同余式，可以确定 nn 是合数。

但反过来不成立：某些合数也可能对一个甚至许多个底数通过 Fermat 测试。因此，Fermat 测试只能快速排除一部分合数，不能单独作为可靠的素性判定。

3.2 二次探测

若 pp 是奇质数，则

x2≡1(modp)x^2\equiv1\pmod p

只有两个解：

x≡1(modp)或x≡−1(modp).x\equiv1\pmod p \qquad\text{或}\qquad x\equiv-1\pmod p.

这是因为

x2−1=(x−1)(x+1),x^2-1=(x-1)(x+1),

而质数 pp 整除乘积时，必然整除其中一个因子。

因此，在不断平方的过程中，若一个数在模 pp 下最终变成 11，它此前只能从 11 或 −1-1 平方得到。

3.3 强伪素数测试

对待测奇数 nn，将

n−1=d⋅2sn-1=d\cdot2^s

写成 dd 为奇数的形式。

对一个底数 aa，先计算

x=ad mod n.x=a^d\bmod n.

若 nn 是质数，则必须满足以下情况之一：

  • x=1x=1；
  • x=−1x=-1；
  • 在之后至多 s−1s-1 次平方中出现 −1-1。

若这些情况都没有出现，则 nn 一定是合数。

3.4 固定底数

一般的 Miller–Rabin 是概率算法：随机选择若干底数，合数通过所有测试的概率会快速降低。

在 long long 范围内，可以使用一组固定底数完成确定性判定。当前模板使用前十二个质数：

2,3,5,7,11,13,17,19,23,29,31,37.2,3,5,7,11,13,17,19,23,29,31,37. 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
bool miller_rabin( ll n ) {
    if( n < 2 ) return 0;
    if( n == 2 || n == 3 ) return 1;
    if( n % 2 == 0 || n % 3 == 0 ) return 0;

    ll d = n - 1, s = 0;
    while( d % 2 == 0 ) d >>= 1, s ++;

    static const ll primes[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 };
    for( ll p : primes ) {
        if( n == p ) return 1;

        ll x = qpow( p, d, n );
        if( x == 1 || x == n - 1 ) continue;

        bool ok = 0;
        for( int i = 0; i < s - 1; ++ i ) {
            x = mul( x, x, n );
            if( x == n - 1 ) {
                ok = 1;
                break;
            }
        }

        if( !ok ) return 0;
    }
    return 1;
}

设固定测试底数的数量为常数，单次判定主要进行若干次快速幂，复杂度可以看作 O(log⁡n)O(\log n) 次模乘。

4 Pollard–Rho

4.1 寻找因子的基本想法

设 pp 是 nn 的某个未知质因子。

如果构造一个模 nn 的伪随机序列

xk+1=f(xk),x_{k+1}=f(x_k),

那么它同时也在模 pp 意义下形成一个序列。

当两个位置满足

xi≡xj(modp),x_i\equiv x_j\pmod p,

就有

p∣(xi−xj).p\mid(x_i-x_j).

因此

gcd⁡(∣xi−xj∣,n)\gcd(|x_i-x_j|,n)

可能暴露出 nn 的一个非平凡因子。

Pollard–Rho 常用的迭代函数为

f(x)=x2+c(modn).f(x)=x^2+c\pmod n.

由于有限剩余类上的迭代最终一定进入循环，序列形状类似希腊字母 ρ\rho，算法因此得名。

4.2 为什么能够较快发生碰撞

若只考虑模未知因子 pp 后的序列，根据生日悖论，随机取值大约经过 O(p)O(\sqrt p) 次就可能发生碰撞。

当 pp 接近 n\sqrt n 时，这对应约 O(n1/4)O(n^{1/4}) 的期望规模。

这只是随机算法的期望表现，不是严格的最坏复杂度保证。

4.3 批量计算最大公因数

逐步计算

gcd⁡(∣xi−xj∣,n)\gcd(|x_i-x_j|,n)

会产生大量 gcd 调用。

模板将一批差值先乘起来：

V=∏∣xi−xj∣(modn),V=\prod |x_i-x_j|\pmod n,

再计算

gcd⁡(V,n).\gcd(V,n).

只要其中某个差值含有非平凡因子，这个乘积也会含有该因子。代码每累计一段后检查一次 gcd，以减少常数。

如果得到的结果恰好是 nn，说明这一轮没有分离出有效因子，外层会重新随机选择参数再运行。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
mt19937_64 rng( time(0) );

ll pollard_rho( ll n ) {
    if( n % 2 == 0 ) return 2;
    if( n % 3 == 0 ) return 3;

    ll x = 0, y = 0, c = rng( ) % ( n - 1 ) + 1;
    ll val = 1, d = 1;

    for( int k = 1; ; k <<= 1, y = x, val = 1 ) {
        for( int i = 0; i < k; ++ i ) {
            x = ( mul( x, x, n ) + c ) % n;
            val = mul( val, abs( x - y ), n );

            if( i % 127 == 0 ) {
                d = gcd( val, n );
                if( d > 1 ) return d;
            }
        }

        d = gcd( val, n );
        if( d > 1 ) return d;
    }
}

pollard_rho(n) 的返回值可能是一个非平凡合因子，也可能在失败时返回 nn，所以它的结果仍需要继续判断与拆分。

5 递归质因数分解

5.1 递归结构

对当前整数 nn：

  1. 若 n<2n<2，停止；
  2. 若 Miller–Rabin 判断 nn 为质数，记录它；
  3. 否则用 Pollard–Rho 找到因子 pp；
  4. 递归分解 pp 与 n/pn/p。

这与算术基本定理完全对应：每次把一个合数拆成更小的因子，最终所有叶子都会成为质数。

5.2 去除重复因子

当前模板的目标是求最大质因子。

找到因子 pp 后，代码先不断令 n←n/pn\gets n/p，一次去除 pp 的全部重复贡献，再递归处理剩余部分与 pp。

这样即使原数中含有较高次幂，也不需要对同一个因子重复递归。

5.3 最大质因子接口

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
ll max_factor;

void _find_fac( ll n ) {
    if( n <= max_factor || n < 2 ) return;

    if( miller_rabin( n ) ) {
        max_factor = max( max_factor, n );
        return;
    }

    ll p = n;
    while( p >= n ) p = pollard_rho( n );

    while( n % p == 0 ) n /= p;

    _find_fac( n );
    _find_fac( p );
}

ll get_max_factor( ll n ) {
    max_factor = 0;
    _find_fac( n );
    return max_factor;
}

其中剪枝条件 n <= max_factor 表示：即使当前 nn 本身是质数，也不可能更新已经得到的最大质因子。

当前接口只返回最大质因子。如果题目要求完整质因数集合或质因数指数，需要在质数递归终点记录因子，再对结果排序和计数。

6 完整模板

下面汇总当前模板中的全部实现。代码保留原有命名与接口，不再增加额外封装。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
mt19937_64 rng( time(0) );

ll mul( ll a, ll b, ll mod ) { 
    return (ll)( (i128)a * b % mod ); 
}

ll qpow( ll a, ll k, ll mod ) {
    ll res = 1 % mod; a %= mod;
    for( ; k; k >>= 1, a = mul( a, a, mod ) )
        if( k & 1 ) res = mul( res, a, mod );
    return res;
}

bool miller_rabin( ll n ) {
    if( n < 2 ) return 0;
    if( n == 2 || n == 3 ) return 1;
    if( n % 2 == 0 || n % 3 == 0 ) return 0;

    ll d = n - 1, s = 0;
    while( d % 2 == 0 ) d >>= 1, s ++;

    static const ll primes[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 };
    for( ll p : primes ) {
        if( n == p ) return 1;

        ll x = qpow( p, d, n );
        if( x == 1 || x == n - 1 ) continue;

        bool ok = 0;
        for( int i = 0; i < s - 1; ++ i ) {
            x = mul( x, x, n );
            if( x == n - 1 ) {
                ok = 1;
                break;
            }
        }

        if( !ok ) return 0;
    }
    return 1;
}

ll pollard_rho( ll n ) {
    if( n % 2 == 0 ) return 2;
    if( n % 3 == 0 ) return 3;

    ll x = 0, y = 0, c = rng( ) % ( n - 1 ) + 1;
    ll val = 1, d = 1;

    for( int k = 1; ; k <<= 1, y = x, val = 1 ) {
        for( int i = 0; i < k; ++ i ) {
            x = ( mul( x, x, n ) + c ) % n;
            val = mul( val, abs( x - y ), n );

            if( i % 127 == 0 ) {
                d = gcd( val, n );
                if( d > 1 ) return d;
            }
        }

        d = gcd( val, n );
        if( d > 1 ) return d;
    }
}

ll max_factor;

void _find_fac( ll n ) {
    if( n <= max_factor || n < 2 ) return;

    if( miller_rabin( n ) ) {
        max_factor = max( max_factor, n );
        return;
    }

    ll p = n;
    while( p >= n ) p = pollard_rho( n );

    while( n % p == 0 ) n /= p;

    _find_fac( n );
    _find_fac( p );
}

ll get_max_factor( ll n ) {
    max_factor = 0;
    _find_fac( n );
    return max_factor;
}

7 方法选择

问题规模与目标 方法
单个整数较小 试除法
批量处理 11 到 NN 埃氏筛或线性筛
只判断 long long 范围内的大整数是否为质数 Miller–Rabin
分解 long long 范围内的大合数 Miller–Rabin + Pollard–Rho
只求最大质因子 get_max_factor
求完整质因数分解 在递归终点收集所有质因子

需要注意：

  • mul 依赖 i128，用于避免模乘溢出；
  • Pollard–Rho 是随机算法，不保证每轮都找到有效因子；
  • pollard_rho 的结果不一定是质数；
  • 得到一个因子后，必须继续递归分解；
  • 多次调用 get_max_factor 时，每次都会重新清空 max_factor；
  • 当前模板面向正整数分解，00、负数与符号需要由题目外围单独处理。

8 与前面内容的联系

  • 算术基本定理保证递归分解最终得到唯一的质因数集合与指数；
  • 快速幂和模乘来自模运算的公共底座；
  • Miller–Rabin 使用费马小定理与二次探测；
  • Pollard–Rho 使用最大公因数从伪随机序列的碰撞中提取因子；
  • ExLucas 在模数较大且难以试除分解时，可以使用 Pollard–Rho 获得质数幂分解。

相关内容：

数学 > 数论
#数论/质因数分解
大整数质因数分解
http://miqiu.space/posts/b73241e7/
作者
Miqiu
发布于
2026年7月1日
许可协议