[学习笔记]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>\) 如下 :
那么梅森数 \(2^p - 1\) 是素数当且仅当 :
3. AKS素性测试
利用了如下定理 :
对于一个不小于 \(2\) 的整数 \(n\) 当且仅当 :
满足如下同余多项式时
\(n\) 为质数.
具体过程参见维基百科. \(\mathtt{Link.}\)
非确定性素性测试
1.费马素性测试
首先引入费马小定理 :
对于所有的质数 \(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);
}
