[学习笔记] Miller-Rabin质数测试 & Pollard-Rho质因数分解
目录
Miller-Rabin质数测试 & Pollard-Rho质因数分解
考试遇见卡质因数分解的题了…活久见…毒瘤lun
于是就学了一发qaq
Miller-Rabin质数测试
一个多项式时间的基于随机的质数测试算法.
玄学.
一些依赖的定理
费马小定理
若 \(p \in \mathbb P\) 则 \(\forall a\in \mathbb N^+, a^{p-1} \equiv 1\pmod p\)
它的逆命题基本上是真的, 于是我们可以把这个定理的逆命题作为判断质数的依据.
为啥说”基本上”是真的呢? 因为有些鬼畜的合数 \(n\) 满足 \(\forall a\in \mathbb N^+,a^{n-1}\equiv 1 \pmod n\). 这种鬼畜的合数被称为 Carmichael 数. 前三个 Carmichael 数是 \(561\), \(1105\) 和 \(1729\). OEIS数列编号 A002997 注意到它们可以挺小的…所以只用这个定理的逆命题测试有很大的问题.
二次探测引理
这个名字好像不是很通用…但是我们在这里先这么叫好了…
若 $p\in \mathbb P $ , 则 \(1\bmod p\) 不存在非平凡二次剩余.
我们尝试用它来作为判质数的条件, 那么就可以:
若 \(\forall a^2 \equiv 1 \pmod p, a\equiv \pm1 \pmod p\), 则 \(p \in \mathbb P\)
这个命题也基本上是真的…
但是也有些强伪质数满足上面这个条件…
不过好像两个合起来用就没毛病了?
具体用法是先令 \(n-1=2^tu\), 其中 \(u\) 是奇数. 然后随机一个值 \(a\), 把 \(a^u\) 平方 \(t\) 次. 如果某次平方出 \(1\) 了, 那么判断一下平方前的值是否是 \(\pm 1\). 如果是那么就算通过了二次探测. 通过二次探测后刚好求出了 \(a^{2^tu}\) 也就是 \(a^{n-1}\), 判一下是不是 \(1\) 就可以验证费马测试. 如果都是, 那么就称 \(n\) 通过了 \(a\) 的测试, 或者说 “\(a\) 不是 \(n\) 是合数的证据”.
实现以及正确率
具体过程大概是
- 把乱七八糟的trival情况判掉, 比如 \(\le 2\) 啊, 偶数啊啥的.
- 进行 \(s\) 轮测试, 每次随机一个底数 \(a\) 进行上文所述的合并测试.
正确率大概要靠下面这个结论:
若 \(n\) 是一个奇合数, 那么 \(n\) 为合数的证据数量至少为 \(\frac{n-1}2\).
具体证明不详细展开了…算法导论上证了…
那么也就是说, 对于一个合数我们随机选择一个 \(a\) 然后刚好命中不是证据的 \(a\) 的概率小于 \(\frac 1 2\). 进行 \(s\) 轮测试后均命中非证据的概率小于 $2^{-s} $.
实际上不是证据的数的数量远远达不到 \(\frac {n-1} 2\). 能证明的最好结论是非证据的 \(a\) 最多有 \(\frac {n-1} 4\) 个. 而且这个界是能达到的.
所以实际正确率大概是 \(1-4^{-s}\). 一般取 \(s=10\) 效果就已经很好了.
一般用的时候要快速乘. 龟速乘(倍增加法)多半会GG.
代码实现
下面是板子题 LOJ #143. 质数判定 的AC代码
#include <bits/stdc++.h>
typedef long long intEx;
bool MillerRabin(intEx);
inline intEx RandInt(intEx,intEx);
inline intEx Mul(intEx,intEx,intEx);
inline intEx Pow(intEx,intEx,intEx);
int main(){
for(intEx x;scanf("%lld",&x)!=EOF;){
if(MillerRabin(x))
puts("Y");
else
puts("N");
}
return 0;
}
bool MillerRabin(intEx n){
static const int ROUND=10;
if(n<2)
return false;
if(n==2)
return true;
if(!(n&1))
return false;
int p=0;
intEx m=n-1;
while(!(m&1))
++p,m>>=1;
for(int k=0;k<ROUND;k++){
intEx pw=Pow(RandInt(1,n-1),m,n);
intEx last=pw;
for(int i=0;i<p;i++){
pw=Mul(pw,pw,n);
if(pw==1&&last!=1&&last!=n-1)
return false;
last=pw;
}
if(pw!=1)
return false;
}
return true;
}
inline intEx Mul(intEx a,intEx b,intEx p){
intEx t=a*b-(intEx)((long double)a*b/p+0.5)*p;
return t<0?t+p:t;
}
inline intEx RandInt(intEx l,intEx r){
static std::mt19937_64 mt(int(new int));
return std::uniform_int_distribution<intEx>(l,r)(mt);
}
inline intEx Pow(intEx a,intEx n,intEx p){
intEx ans=1;
while(n>0){
if(n&1)
ans=Mul(ans,a,p);
a=Mul(a,a,p);
n>>=1;
}
return ans;
}
Pollard-Rho质因数分解
这个算法也是个概率算法. 不过它运行时间是期望的…
首先它也有个依赖的结论
生日悖论与生日攻击
首先是生日悖论. 一个 \(23\) 人的团体存在两人生日相同的概率要大于 \(50\%\). 一个推论是在 \(n\) 个值中随机选取若干个值, \(O(\sqrt n)\) 次后就会有很大概率产生某个值与之前选的值重复的情况.
大概证明可以长这样:
我们补集转化一下, 转而求选出的 \(k\) 个值都不同的概率. 显然它应该长这样:
\[
P_n(k)=\prod_{i=0}^{k-1} \left(1-\frac i n\right)
\]
我们只要让这个值小于 \(\frac 1 2\) 就好了. 而由泰勒展开可得:
\[
\exp(x)=1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\cdots
\]
那么有:
\[
1+x < \exp(x)
\]
于是就有:
\[
P_n(k)=\prod_{i=0}^{k-1}\left(1-\frac i n\right) < \prod_{i=0}^{k-1}\exp\left(-\frac i n\right) = \exp\left(-\sum_{i=0}^{k-1} \frac i n\right) =\exp\left(\frac{-k(k-1)}{2n}\right)
\]
那么我们只要让不等式右边小于 \(\frac 1 2\) 就好了. 那么我们有:
\[
\exp\left(\frac{-k(k-1)}{2n}\right) < \frac 1 2
\]
两边取对数解一下就有:
\[
k^2-k>2n\ln2
\]
又因为 \(\ln 2\) 是个常数, 于是 \(k=O(\sqrt n)\).
主要思想
设我们要对 \(n\) 进行因数分解, \(n\) 存在一个非平凡因子 \(p\) 且 \(p\le n/p\).
我们取一个次数至少为 \(2\) 且包含常数项的多项式函数 \(f(x)\) (比如 \(f(x)=x^2+c\), \(c\) 随机取), 设 \(g(x)=f(x)\bmod n\), 用 \(g\) 来生成伪随机数. 方法是随机选取一个初始值 \(r\) 不断将新的 \(x\) 代入 \(g(x)\) , 也就是 \(x_1=g(r), x_2=g(g(r)),x_3=g(g(g(r)))\) . 那么我们可以认为数列 \(\langle x_k\rangle\) 看上去是随机的. 而且 \(x_k=g(x_{k-1})\).
也就是说这个伪随机数列只和它的上一个输出有关. 从一个相同的初值可以得到相同的序列.
由于 \(p\mid n\), 所以数列 \(\langle x_k \bmod p\rangle\) 也满足这个性质, 即 \(x_k\equiv g(x_{k-1}) \pmod p\). 因为这个伪随机数生成器的状态只和上一个输出有关, 那么只要生成的新值命中了之前已经生成的值, 那么这个随机数列就会出现环.
由于 \(p\) 比较小, \(\langle x_k \bmod p \rangle\) 有很大概率比 \(\langle x_k \rangle\) 更早出现环 (只要新生成的 \(x\) 对 \(p\) 取模的值命中以前出现的值就会出现环, 而 \(p\) 不会超过 \(\sqrt n\)). 根据生日悖论, 在大约 \(\sqrt p\) 次随机后就会有很大概率出现. 如果出现这种情况, 那么我们就获得了两个值 \(x_a \equiv x_b \pmod p\). 于是它们之间的差是可以被 \(p\) 整除的. 我们取 \(p=\gcd (\left|x_a-x_b\right|,n)\) 即可.
至于这个判环的过程, 我们可以使用神奇的Floyd判环法. 建立两个哨兵 \(x_a\) 和 \(x_b\), 每次令 \(x_a\) 前进一步, 令 \(x_b\) 前进两步, 如果 \(x_a=x_b\), 则说明走到了环上.
但是我们并不知道 \(p\) 是多少, 所以并不能直接判断 \(x_a = x_b\) 来判断 \(\langle x_k \bmod p\rangle\) 是否进入环. 但是根据上文所述, 如果出现了 \(x_a\equiv x_b \pmod p\), 那么就一定出环了, 这个时候 \(\gcd (\left|x_a-x_b\right|,n)\neq 1\), 我们便成功获得了一个 \(n\) 的非平凡因子.
不过有时候我们可能会手气不佳, 抽到的 \(c\) 和 \(r\) 生成的 \(\langle x_k \rangle\) 和 \(\langle x_k \bmod p\rangle\) 可能会同时进入环(显然前者不可能比后者进环更快, 最多同时). 不难发现这个环长也是 \(\sqrt p\) 级别的. 但是由于这时候找到的 \(x_a\) 和 \(x_b\) 是相等的, 于是并不能带来什么信息. 在这个时候有两种可能, 第一是 \(n\in \mathbb P\), 第二是脸黑.
普通的Pollard-Rho会到此为止而不提供任何信息.
期望的时间复杂度是 \(O(n^{1/4}\log n)\) 的, 那个 \(\log\) 是 \(\gcd\) 的时候来的.
具体实现
首先上文所述的过程只能用来找非平凡因子而非质因子.
我们把它和 Miller-Rabin 质数测试结合起来就可以进行质因数分解了. 大致流程如下:
- 用 Miller-Rabin 判断当前 \(n\) 是否是质数. 如果是, 丢进答案集合, 跑路.
- 随机选择一对 \(c\) 和 \(r\), 进行上文所述的测试过程.
- 如果找到了一个非平凡因子 \(p\), 递归求解 \(p\) 和 \(n/p\), 然后跑路.
- 否则, 承认自己脸黑, 返回第 2 步再来一次.
代码片段
void PollardRho(intEx n,std::vector<intEx>& factor){
if(MillerRabin(n))
factor.push_back(n);
else while(true){
intEx a,b,c=RandInt(0,n-1);
a=b=RandInt(0,n-1);
auto f=[=](intEx x){
return (Mul(x,x,n)+c)%n;
};
b=f(b);
while(a!=b){
intEx delta=std::abs(a-b);
delta=std::__gcd(delta,n);
if(delta!=1){
PollardRho(delta,factor);
PollardRho(n/delta,factor);
return;
}
a=f(a);
b=f(f(b));
}
}
}
洛谷板子过于SXBK地卡常于是没过TwT…我活该大常数…
LOJ板子更加SXBK地把数据出到了 \(1\times 10^{30}\)…跑路了QAQ…
一般应用没啥问题(吧)…