miller-rabin 和 pollard-rho

发表于 更新于 分类于 阅读次数：

miller-rabin

用于快速测试一个数 \(n\) 是否为质数，有概率出错。

不妨假设 \(n\) 是奇数，此外 miller-rabin 需要随意选取一个小质数 \(p\) 。

理论

miller-rabin 用到的两个基本定理：

  • 费马小定理：如果 \(n\) 为质数，一定有 \(p^n \equiv p \pmod{n}\) ，即 \(p^{n-1} \equiv 1\) 。

  • 二次探测定理：如果 \(n\) 为质数，对于 \(a^2 \equiv 1 \pmod{n}\) 一定有 \(a \equiv -1\) 或 \(a \equiv n - 1\) 。

如果 \(p^{n-1} \equiv 1\) 不成立，根据费马小定理， \(n\) 一定不为质数，
否则令 \(x = n - 1\) ，\(x\) 一定是偶数，那么可以计算 \(y = p^{\frac{x}{2}}\) ，
由于 \(p^x \equiv 1\) ，那么根据二次探测定理，如果 \(y \equiv 1\) 和 \(y \equiv n - 1\) 都不成立，\(n\) 一定不为质数。
否则如果 \(y \equiv 1\) ，令 \(x = \frac{x}{2}\) 继续判断，
直到 \(y \equiv n - 1\) 或者 \(x\) 为奇数，此时 \(n\) 通过了测试，但无法保证 \(n\) 是质数还是合数。

如果仅仅通过一次测试就认为 \(n\) 是质数的话错误的概率很大，
实践中往往选取多个质数 \(p\) 测试多次，全部通过才认为 \(n\) 是质数。
实践证明，只要取前 9 个质数，就能准确判断 \(10^{18}\) 内的所有数。

实现

参考实现（经过测试）：

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
inline ll mul(ll x, ll y, ll mod) {
    ll tmp = x * y - ll((long double)x / mod * y + 0.5) * mod;
    return tmp < 0 ? tmp + mod : tmp;
}

inline ll power(ll x, ll k, ll mod) {
    ll res = 1;
    while (k) {
        if (k & 1)
            res = mul(res, x, mod);
        x = mul(x, x, mod);
        k >>= 1;
    }
    return res % mod;
}

inline bool miller(ll n, int p) {
    if (power(p, n - 1, n) != 1)
        return false;
    ll x = n - 1;
    while (not(x & 1)) {
        x >>= 1;
        ll y = power(p, x, n);
        if (y == n - 1)
            return true;
        else if (y != 1)
            return false;
    }
    return true;
}

bool is_prime(ll n) {
    static const int len = 9;
    static const int p[len] = { 2, 3, 5, 7, 11, 13, 17, 19, 23 };
    for (int i = 0; i < len; i++)
        if (n == p[i])
            return true;
    for (int i = 0; i < len; i++)
        if (not miller(n, p[i]))
            return false;
    return true;
}

优化

事实上 miller-rabin 中进行的多次快速幂是没有必要的。

将 \(n\) 表示为 \(n=2^s t + 1\) 的形式，其中 \(t\) 为奇数，
那么 miller-rabin 中可能用到的所有快速幂，都是 \(p^{2^it} (0 \leq i < s)\) 的形式。
原本 \(x\) 是从 \(2^s t\) 到 \(t\) ，考虑倒着进行，\(x\) 从 \(t\) 到 \(2^s t\) ，
这样只需要一次快速幂计算 \(y = p^t\) ，每次 \(x\) 扩大一倍时将 \(y\) 平方即可。

参考实现（经过测试）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
bool miller(ll n, int p) {
	ll t = n - 1;
	while(!(t & 1)) t >>= 1;
	ll y = power(p, t, n), ny;
	while((t << 1) < n) {
		ny = mul(y, y, n);
		if(ny == 1 and !(y == 1 or y == n - 1)) return 0;
		t <<= 1;
		y = ny;
	}
	return y == 1;
}

bool is_prime(ll n) {
	static const int len = 9, p[len] = {2, 3, 5, 7, 11, 13, 17, 19, 23};
	if(n == 1) return 0;
	for(int i = 0; i < len; i ++) if(n == p[i]) return 1;
	if(!(n & 1)) return 0;
	for(int i = 0; i < len; i ++) if(!miller(n, p[i])) return 0;
	return 1;
}

pollard-rho

用于快速找到一个数 \(n\) 的一个非平凡因子，常应用于大数的质因数分解。

这玩意大家似乎都是背代码，很难找到关于原理以及复杂度的证明，学习资料很少，如果有误，还请指正。

非平凡因子

关于非平凡因子：如果 \(d\) 能整除 \(n\) 且满足 \(1 < d < n\) ，则称 \(d\) 是 \(n\) 的非平凡因子。

从定义可以看出质数是没有非平凡因子的，事实上质数直接拿来做 pollard-rho 会陷入死循环，
因此 pollard-rho 常常需要配合 miller-rabin 使用，在一个数为质数时停止分解。

因此以下讨论假定 \(n\) 是合数。

随机实验

随机一个正整数 \(d(d < n)\) ，如果 \(d\) 是 \(n\) 的因子，那么就找到了 \(n\) 的一个非平凡因子 \(d\) 。
然而这样的概率很小，最坏情况下是 \(\frac{1}{n}\) 的。

随机一个正整数 \(v(v < n)\) ，求出 \(d = \gcd(v, n)\) ，如果 \(d > 1\) 那么就找到了 \(n\) 的一个非平凡因子 \(d\) 。
这样成功的概率就大得多，最坏是 \(\frac{1}{\sqrt{n}}\) 的。

证明：
设 \(p = minp(n)\) 是合数 \(n\) 的最小质因子，那么 \(p \leq \sqrt{n}\) ，而 \(p\) 的倍数有 \(\lfloor\frac{n}{p}\rfloor\) 个，
因此概率的一个下界是 \(\frac{1}{p}\) ，这是取到 \(p\) 的倍数的概率。
而这个下界是可以达到的，当 \(n = p^2\) 且时，概率恰为 \(\frac{1}{p}\) ，也就是 \(\frac{1}{\sqrt{n}}\) 。

以下讨论中均有 \(p\) 是 \(n\) 的最小质因子。

生日悖论

然而这个概率还是很小，期望下需要随机 \(O(p)\) 次才能成功（这是个上界）。

根据生日悖论可以提供一个优化思路，
随机 \(\sqrt{p}\) 个数 \(v\) ，其中有 \(\frac{1}{2}\) 的概率存在两个数 \(v_1, v_2\) 模 \(p\) 同余，
也就是说存在两个数的差是 \(p\) 的倍数。
但是如果两两作差与 \(n\) 取 gcd ，期望复杂度仍然是 \(O(p)\) 的。

rho

pollard-rho 用到了一个特殊的伪随机数列 \(\{a\}\) ，对于参数 \(c\) ，其满足：

\[ a_{i+1} = (a_i^2 + c) \bmod n \]

这个数列是个 \(\rho\)(rho) 形，这也是 pollard-rho 的名字由来。
换言之，存在链长 \(T\) 和环长 \(M\) ，满足 \(\forall i,k>0: a_{T+i} = a_{T+i+kM}\) 。
如果把 \(a\) 两两不同的前缀部分近似看做真随机数列，那么根据生日悖论，这个长度 \(T + M\) 是期望 \(O(\sqrt{n})\) 的。
也就是说 \(T\) 和 \(M\) 期望都是在 \(O(\sqrt{n})\) 以内的。

令 \(b_i = a_i \bmod p\) ，当然由于 \(p\) 未知，\(b\) 这个数列实际上是未知的。
那么同样的道理，\(b\) 的链长和环长 \(t, m\) 是期望 \(O(\sqrt{p})\) 的。

现在要做的是通过对 \(a\) 作差来间接对 \(b\) 作差，目标是选到两个下标 \(i, j\) 满足 \(b_i = b_j\) ，
这样的话，\(|a_i - a_j|\) 就是 \(p\) 的倍数，也就是说 \(\gcd(n, |a_i - a_j|) \geq p\) 。
关键在于环长，如果只选取下标 \((i, 2i)\) 作差，那么当 \(i \geq t\) 且 \(i=km\) 时，
根据 \(b_i = b_{t+(i-t)} = b_{t+(i-t)+km} = b_{t+(i-t)+i} = b_{2i}\) ，
可以确定此时一定有 \(b_i = b_{2i}\) ，
而上面提到过 \(t, m\) 期望都是 \(O(\sqrt{p})\) 的，那么第一个这样的 \(i\) 的大小也是期望 \(O(\sqrt{p})\) 的。

一个期望复杂度 \(O(\sqrt{p})\) 的 pollard-rho 基本思想大抵如此。

实践

算法流程就是随机一个生成序列的参数 \(c\) 。
然后枚举 \(i\) ，同时通过递推维护 \(a_i\) 和 \(a_{2i}\) ，
当 \(\gcd(n, |a_i - a_{2i}|) > 1\) 时就得到了一个非平凡因子。
直接做考虑 gcd 期望复杂度 \(O(\sqrt{p} \log n)\) 。
但这个 gcd 很好优化，设 \(v_i = |a_i - a_{2i}|\) ，把 \(v\) 相邻 \(k\) 个乘起来再和 \(n\) 做 gcd ，
根据 \(\gcd(x, n) = \gcd(x \bmod n, n)\) ，乘积可以是模 \(n\) 意义下的。
理论上只要 \(k > O(\log n)\) ，复杂度就能做到 \(O(\sqrt{p} + k + \log n)\) ，
一般认为 \(p\) 足够大，此时复杂度就是 \(O(\sqrt{p})\) ，最坏情况下就是 \(O(n^{\frac{1}{4}})\) 。
实践中往往采用倍增的方式，起初 \(k\) 不断成倍增大，增大到 128 时保持不变。
这样可以保证在 \(p\) 很小的时候不用算 \(\log n\) 次就能得到非平凡因子，减小常数。

有一些需要注意的地方：

  • 有小概率 \(M\) 与 \(m\) 恰好相等，在 \(b_i = b_{2i}\) 的时候恰好有 \(a_i = a_{2i}\) ，此时应该更换参数 \(c\) 重新进行算法。

  • 当 \(v\) 的乘积模 \(n\) 为 0 时需要及时退出，否则算法会返回 \(n\) ，而 \(n\) 不是 \(n\) 的非平凡因子。

  • \(n = 4\) 时 无论怎么取参数 \(c\) 和序列初值，\(M = m\) 恒成立，需要特判。

参考实现（经过测试，上接 miller-rabin 的部分）：

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
inline ll nxt(ll x, ll c, ll n) {
	ll res = mul(x, x, n) + c;
	return res >= n ? res - n : res;
}

ll pollard(ll n) {
	if(n == 4) return 2;
	ll c = 1ll * rand() * rand() % n;
	ll a = nxt(rand(), c, n);
	ll a2 = nxt(a, c, n);

	int lim = 2;
	while(a != a2) {
		ll v = 1;
		for(int i = 0; i < lim and a != a2; i ++) {
			ll tov = mul(v, std::abs(a - a2), n);
			if(!tov) return std::__gcd(v, n);
			v = tov;
			a = nxt(a, c, n);
			a2 = nxt(nxt(a2, c, n), c, n);
		}
		ll d = std::__gcd(v, n);
		if(d > 1) return d;
		if(lim < 128) lim <<= 1;
	}

	return pollard(n);
}

应用

求最小质因数，最大质因数，质因数分解等，
都只需要每次找到 \(n\) 的非平凡因子 \(d\) 然后不断令 \(n\) 除以 \(d\) 转换为子问题，
需要注意的是 \(d\) 不一定是质数，也要递归处理。
终止条件是 \(n\) 为质数。