[学习笔记]Miller-Rabin 素性测试与 Pollard-Rho 算法

前言#

新年新气象!

sto 神 $\mathrm{{\color{black}{\_}}{\color{red}{slb}}}$ orz

他爆卷 Miller Rabin 和 Pollard Pho.

素性测试#

素性测试是检验给定数是否是素数的测试.

确定性素性测试#

0. 朴素枚举#

根据素数的定义,可以朴素地枚举 $[2,n)$ 中的每一个数并判断其是否是 $N$ 的因数,时间复杂度 $\mathcal{O}(n)$.

1. 试除法#

显然朴素的枚举连一种算法也谈不上,最多就是一种方法.

于是考虑因数的性质,其成对出现.

那么只需要枚举 $[2,\sqrt{n}]$ 的每个数然后判断其是否为 $n$ 的因数即可.

inline bool IsPrime(int x) {
    if(x < 2) return 0;
    for(int i = 2; i * i <= x;++i)
        if(!(x % i)) return 0;
    return 1;
}

2. 卢卡斯-莱默检验法#

这种素性测试仅能检验梅森数.

梅森数 : 形如 $2^p - 1$ 的数,其中 $p$ 为质数.

首先定义序列 $<s_n>$ 如下 :

$$\large s_n = \left\{\begin{matrix} &4\ (n = 0)\\ &s_{n - 1}^{2} - 2\ \mathrm{otherwise.}\\ \end{matrix}\right.$$

那么梅森数 $2^p - 1$ 是素数当且仅当 :

$$\large s_{p - 2} \equiv 0 \pmod{2^p - 1}$$

3. AKS素性测试#

利用了如下定理 :

---

对于一个不小于 $2$ 的整数 $n$ 当且仅当 :

$$\large \forall a \perp n,a \in \mathbb{Z}$$

满足如下同余多项式时

$$\large (x + a)^n \equiv (x^n + a) \pmod n$$

$n$ 为质数.

---

具体过程参见维基百科. $\mathtt{Link.}$

非确定性素性测试#

1.费马素性测试#

首先引入费马小定理 :

对于所有的质数 $p$,

$$\large a^p \equiv a \pmod p$$

对于互质的情况,将其变形为 :

$$\large a^{p - 1} \equiv 1 \pmod p$$

那么我们不断选取整数 $a$,判断等式是否成立.

inline int QuickPower(int a,int b,int p) {
	int res = 1;
	while(b) {
		if(b & 1) res = (ll)res * a % p;
		a = (ll)a * a % p;
		b >>= 1;
	}
	return res;
}

bool IsPrime(int n) {
	if(n < 3) return n == 2;
	rep(i,1,16) {
		int a = rand() % (n - 2) + 2;
		if(QuickPower(a,n - 1,n) != 1)
			return 0;
	}
	return 1;
}

这个过程的时间复杂度是 $\mathcal{O} (T \log n)$ 的,其中 $T$ 为检验次数,次数过少则不能保证一定的准确性.

但是费马小定理的逆定理并不具有正确性.

对于 卡迈克尔数 而言,满足以上性质的同时,其还为一个合数.

更遗憾的是,卡迈克尔数是可证明无穷多的,且并不是非常稀疏,这导致不能靠一张卡迈克尔数表来完善测验过程.

2.米勒-拉宾素性测试#

Miller Rabin素性测试最初是基于广义黎曼猜想,由 Gary Lee Miller 提出.

但因为广义黎曼猜想至今未得到证明,Michael Oser Rabin 优化了这个方法,使得其不许依赖该假设,发明了 Miller Rabin 素性测试.

二次探测定理 : 对于奇素数 $p$ ,方程 $x^2 \equiv 1 \pmod p$ 的两解为 $x \equiv 1 \pmod p$ 和 $x \equiv p - 1 \pmod p$

然后把二次探测定理与费马素性测试的过程结合.

首先偶数除 $2$ 以外的偶数都是合数,那么在检测的时候 $n - 1$ 是个偶数.

考虑把 $n - 1$ 分解为 $u \times 2^t$

设选取的底数为 $a$,令 $d=p-1$.

判断 $a^d\bmod{p}$ 是否为 $1$,若非 $1$ 则 $p$ 为合数.

否则,重复 $d \leftarrow \frac{d}{2}$ 直到 $d$ 为奇数或 $a^d\not\equiv 1\pmod{p}$

最后，若 $a^d \not\equiv\pm 1\pmod{p}$，则 $p$ 为合数,否则基本可以确认 $p$ 是素数.

然后是关于误判 :

选定一个底并进行一次测试,误判概率最大为 $\frac{1}{4}$.

那么随机选择 $k$ 个底的误判概率就是 $(\frac{1}{4})^k$ 显然来个几十次准确率已经非常高了.

但是依然要1选择几十个,考虑被检测数值域对检测的影响.

实际上只使用前 $12$ 个质数就可以保证检测 unsigned long long 范围内均有确定性了.

这个方法最大的问题可能就是快速幂容易溢出,在检验的数较大时不要吝啬使用 __int128.

constexpr ll base[13] = {0,2,3,5,7,11,13,17,19,23,29,31,37};

inline ll qpow(ll a,ll b = MOD - 2,ll p = MOD) {
	ll res = 1;
	while(b) {
		if(b & 1) res = res * a % p;
		a = a * a % p;
		b >>= 1;
	}
	return res;
}

inline bool check(ll a,ll p) {
	ll d = p - 1,pw = qpow(a,d,p);
	if(pw != 1) return 0;
	while(!(d & 1)) {
		d >>= 1,pw = qpow(a,d,p);
		if(pw ==  p - 1) return 1;
		else if(pw != 1) return 0;
	}
	return 1;
}

bool IsPrime(ll p) {
	if(p < 40) {
		rep(i,1,12) if(p == base[i])
			return 1;
		return 0;
	}
	rep(i,1,12) if(!check(base[i],p))
		return 0;
	return 1;
}

Pollard-Rho 算法#

用于在 $\mathcal{O} (n^{\frac{1}{4}})$ 期望时间内求合数 $n$ 的某个非平凡因子.

(更加正确的写法应该是 : Pollard $\rho$ 算法)

(卡常算法毁我青春)

先讲明白什么是非平凡因子 : 平凡因子就是最普适的因子,指 $1$ 和 $n$ 本身,非平凡因子就是其余的因子.

假设你拿到一个 $(10^{100},10^{101})$ 这个范围的很平凡的一个数,有人拿着枪指着你的头,让你在一年内不靠外界工具,不使用纸笔,求出其任意一个非平凡因数,你可以任意次给他答案直到猜出来或者时间结束,你该怎么办?

还能怎么办,猜呗,反正猜着猜着肯定能猜到一个,这就是 Pollard-Rho.

首先我们要找的数一定在 $n$ 以内,尝试构造一个伪随机序列 :

$f_n = (f_{n - 1}^2 + c) \bmod n$

然后可以发现这个序列是循环的,循环节大小大约 $\mathcal{O}(\sqrt{n})$.

令 $n$ 的最小非平凡因子为 $m$ ,那么有 $m \in [2,\sqrt{n}]$,现在可以生成一个新序列 $g_i = f_i \bmod m$,现在序列中不同数字个数为 $\mathcal{O} (n^{\frac{1}{4}})$ 的了,且推一下式子发现 $g_n$ 这个序列满足类似的递推关系 : $g_n = (g_{n - 1}^{2} + c) \bmod m$

然后考虑 $m$ 是如何求出来的.

首先由于 $m < n$ , 一定存在一对 $(i,j)$ 使得 $f_i \ne f_j$ 且 $g_i = g_j$.

那么这时候 $n \not | \ |f_i - f_j|$ 且 $m | \ |g_i - g_j|$,那么求出来 $\gcd(n,|f_i - f_j|)$ 就可以得到一个 $n$ 的非平凡质因子了.

那么该如何快速求出这样的 $(i,j)$ 呢?求其实并不重要,都是 $\mathcal{O} (n^{\frac{1}{4}})$ 的级别了,求这个数对不如靠枚举.

首先看两个序列,对 $n$ 取模的,循环节一定比对 $m$ 取模的长,且就像同样速度跑在跑道内外圈一样会相遇.

那么就枚举这个伪随机的序列直到出现环,如果没有在途中得到一个非平凡因子,就调整伪随机数的常数 $c$ 然后继续尝试.

具体过程就是每次求一个 $\gcd(|f_i - f_j|,n)$ , 如果得到一个 $n$ 的平凡因子则完成过程,否则继续,然后判环其实就是个双指针,一个每次走一步,另一个每次走两步,相遇就找到环了.

但是众所周知辗转相除是带一个 $\log$ 的,于是我们求乘积和 $n$ 的 $\gcd$ 累计 $2^k - 1$ 次后求一次 $\gcd$ 这就是倍增优化.

然后再来一个黑科技(?) gcd()

众所周知有一个求 $\gcd$ 的算法是 更相减损法,这个算法求 $\gcd(a,b)$ 第一步是尽量将 $a,b$ 除以 $2$, 最后将 $\gcd$ 乘以 $2$ 的对应次幂.

那么这个过程可以很快地实现,首先是判断 $a,b$ 能除以 $2$ 几次,众所周知整除 $2$ 就是二进制末位为 $1$ 的时候二进制右移一位,那么求出 $a$ 按位或 $b$ 的二进制末尾连续 $0$ 个数即可,这个过程通过跑得飞快的内置函数 __builtin_ctzll() 实现.

模板 : $\mathtt{Link.}$

并不想写龟速乘,于是使用了 __int128.

但是用了 __int128 做乘法就不能使用 Barrett Reduction 了.有得必有失

据说 C++ 的编译器对常数取模是有优化的,所以固定模数一定要用 const / constexpr.

Code :

inline ll qpow(ll a,ll b = MOD - 2,ll p = MOD) {
	ll res = 1;
	for(;b;b >>= 1) {
		if(b & 1) res = (I128)res * a % p;
		a = (I128)a * a % p;
	}
	return res;
}

inline ll gcd(ll a,ll b) {
	if(!a || !b) return a | b;
	ll k = __builtin_ctzll(a | b);
	a >>= __builtin_ctzll(a);
	b >>= __builtin_ctzll(b);
	ll p = a % b;
	while(p)
		a = b,b = p,p = a % b;
	return b << k;
}

constexpr ll base[13] = {0,2,3,5,7,11,13,17,19,23,29,31,37};

inline bool check(ll a,ll p) {
	ll d = p - 1,pw = qpow(a,d,p);
	if(pw != 1) return 0;
	while(!(d & 1)) {
		d >>= 1,pw = qpow(a,d,p);
		if(pw ==  p - 1) return 1;
		else if(pw != 1) return 0;
	}
	return 1;
}

inline bool IsPrime(ll p) {
	if(p < 40) {
		rep(i,1,12) if(p == base[i])
			return 1;
		return 0;
	}
	rep(i,1,12) if(!check(base[i],p))
		return 0;
	return 1;
}

ll ans;

ll PollardRho(ull x) {
	ll s = 0,t = 0;
	ll c = rand() % (x - 1) + 1;
	ll val = 1;
	for(ll fac = 1;;fac <<= 1,s = t,val = 1) {
		for(ll i = 1;i <= fac;++i) {
			t = ((I128)t * t + c) % x;
			val = (I128)val * std::abs(t - s) % x;
			if(!(i % 127)) {
				ll d = gcd(val,x);
				if(d > 1) return d;
			}
		}
		ll d = gcd(val,x);
		if(d > 1) return d;
	}
	return 1;
}

void resolve(ll x) {
	if(x <= ans || x < 2) return ;
	if(IsPrime(x)) {
		ans = std::max(ans,x);
		return ;
	}
	ll p = x;
	while(p >= x) p = PollardRho(x);
	while(!(x % p)) x /= p;
	resolve(x),resolve(p);
}