Miller-Rabin 素性测试 & Pollard-Rho 算法 学习笔记
Miller-Rabin 素性测试 & Pollard-Rho 算法 学习笔记
素性测试
判断一个数是否是素数。
试除法
一种确定性算法。枚举 \([1,\sqrt n]\) 的每个数检验能否除 \(n\)。复杂度 \(O(\sqrt n)\)。
Fermat 素性测试
简单的概率性素数检验。根据费马小定理：
对于质数 \(p\)，\(\gcd(a,p)=1\)，有 \(a^{p-1}\equiv 1 \pmod p\)。
我们的思路是，不断选取 \([2,n-1]\) 中的数 \(a\)，判断 \(a^{n-1}\equiv1\pmod n\) 是否总成立。
bool is_prime(int n) {
if(n<3) return n==2;
const int test_time=8;
// test_time 为测试次数,建议设为不小于 8
// 的整数以保证正确率,但也不宜过大,否则会影响效率
for(int i=1;i<=test_time;++i) {
int a=rnd(2,n-1);
if(pow1(a,n-1,n) != 1) return false;
}
return 0;
}
然而存在合数使得对于 \([2,n-1]\) 的数 \(a\) 都有 \(a^{n-1}\equiv 1\pmod n\)，这样的数称为卡迈克尔数，例如 561 就是最小的卡迈克尔数。所以接下来我们考虑 Fermat 素性测试的优化。
Miller–Rabin 素性测试
二次探测定理：对于奇质数 \(p\)，\(x^2\equiv 1\pmod p\) 的解只有 \(x\equiv 1\pmod p\) 和 \(x\equiv p-1 \pmod p\)。证明只需对方程移项，用平方差公式可得 \((x+1)(x-1)\equiv 0\pmod p\)。
依旧选取 \(a\)，将 \(a^{n-1}\equiv 1\pmod n\) 中的 \(n\) 拆成 \(u\times 2^t\)。先求 \(a^u\bmod n\)，之后对其 \(t\) 次平方，如果发现出现为非 \(n-1\) 或 \(1\) 的平方根即可判断其不是素数。
实现上，先判断 \(a^u\) 是否等于 \(1\)，等于 \(1\) 则此轮结束。否则做 \(t\) 次平方，某一次平方前如果 \(a=n-1\) 则 \(a\) 之后总是 \(1\)，此轮结束。否则 \(t\) 次平方 \(a\) 没有变成过 \(n-1\)，那么考虑第一个 \(1\) 之前的 \(a\) 一定是 \(n\) 不满足条件的根，或者 \(a^{p-1}\) 就不是 \(1\)，所以 \(n\) 不是质数。
可以证明，当 \(n<2^{64}\) 时，把 \(A=[2,3,5,7,11,13,17,19,23,29,31,37]\) 这 12 个 37 以内的质数依次作为 \(a\) 时素性测试一定是正确的。实现时要把 \(A_i\bmod n\) 作为基底 \(a\)，且要满足 \(a\ne\pm1,0\)。
计算量是 \(O(|A|\log n)\)。
int A[12]={2,3,5,7,11,13,17,19,23,29,31,37};
ll pow1(ll x,ll y,ll n) {
ll res=1;
for(;y;y>>=1,x=(it128)x*x%n) if(y&1) res=(it128)res*x%n;
return res;
}
bool is_prime(ll n) {
if(n<3||n%2==0) return n==2;
if(n%3==0) return n==3;
ll u=n-1,t=0;
while(u%2==0) u/=2,++t;
fu(i,0,12) {
ll _a=A[i]%n;
if(_a<=1||_a==n-1) continue;
ll a=pow1(_a,u,n);
if(a==1) continue;
int j=0;
while(j<t) {
if(a==n-1) break;
a=(it128)a*a%n;
++j;
}
if(j==t) return 0;
}
return 1;
}
分解质因数
求 \(n\) 的所有质因子。
试除法
复杂度 \(O(\sqrt n)\)，预处理质数后可以做到 \(O(\frac {\sqrt n}{\ln n})\)。
随机法
以下的算法考虑找出 \(n\) 的因数 \(m\) 并递归处理 \(m,\frac nm\) 直到 \(n\) 为素数。
考虑每次随机一些数，将它们作差，如果差 \(d\) 满足 \(\gcd(d,n)>1\)，则 \((d,n)\) 是因数。
最坏情况下，当 \(n=p^2\) 时，在 \([1,n-1]\) 生成随机数，期望只要随机 \(O(\sqrt p)\) 个数就有两个数 \(i,j\) 使得 \(|i-j|\equiv 0 \pmod p\)。然而两两之间作差使得复杂度还是 \(O(\sqrt n)\) 而且还要额外乘上 \(\gcd\) 的复杂度。
Pollard-Rho
考虑另一种随机方式，我们定义一个伪随机数生成方式，\(f(x)=(x^2+c)\bmod n\)，其中 \(c\) 是一个常数。那么初始化 \(x=0\)，将 \(x\) 与 \(f(x)\) 连边会形成 \(\rho\) 的形状，即环加链。考虑到 \(f(x)\) 有性质：对于 \(|i-j|\equiv 0 \pmod p\)，有 \(|f(i)-f(j)|=|i^2-j^2|=|(i+j)(i-j)|\equiv0\pmod p\)，这意味着对于环上一种距离的点对若满足条件，则所有距离相同的点对都满足条件。
考虑初始让 \(a=0,b=0\)，随机一个 \(c\)，每次让 \(a\) 走一步，\(b\) 走两步，即 \(a\gets f(a),b\gets f(f(b))\)，每次判断是否 \(\gcd(|a-b|,n)>1\)。这样 \(a,b\) 每走一次都在检查一个新的距离，如果 \(a,b\) 再次相遇了则这个 \(\rho\) 里没有得到因子，再随机新的 \(c\) 重复这个过程。
如上所述，假设 \(f(x)\) 均匀随机，分解质因子的复杂度为 \(O(n^{0.25}\times \log n)\)，其中 $\log $ 为 \(\gcd\) 的复杂度。实现如下（函数返回一个非平凡因子）：
mt19937_64 rng(chrono::system_clock::now().time_since_epoch().count());
ll rnd(ll l,ll r) {
return rng()%(r-l+1)+l;
}
ll pollard_rho(ll n) {
if(n==4) return 2; // 注意需要特判，不然会进入 1 和 2 的循环然后永远找不到
if(is_prime(n)) return n;
while(1) {
ll c=rnd(1,n-1);
auto f=[&](ll x) {return ((it128)x*x+c)%n;};
ll a=f(0),b=f(a);
while(a!=b) {
ll d=__gcd(abs(a-b),n);
if(d>1) return d;
a=f(a),b=f(f(b));
}
}
}
浙公网安备 33010602011771号