[学习笔记] Miller-Rabin质数测试 & Pollard-Rho质因数分解

Miller-Rabin质数测试 & Pollard-Rho质因数分解

Miller-Rabin质数测试 & Pollard-Rho质因数分解

考试遇见卡质因数分解的题了...活久见...毒瘤lun

于是就学了一发qaq

Pollard-Rho分解质因数的话需要依赖另一个算法.

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 \]

那么对于 $x> 0$ 有:

\[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)$.

生日攻击是一种现代密码学攻击方法, 就是利用上面的生日悖论来制造Hash碰撞. 主要攻击对象为数字签名. 如果值域为 $U$ 的话期望生成 $\sqrt U$ 种不同信息即可发生碰撞. 如果产生Hash碰撞则可以进行伪造数字签名等等一系列攻击行为. Cookie也需要防范这种情况.

主要思想

设我们要对 $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$, 那么就一定出环了, 这个时候 $p\mid \gcd (\left|x_a-x_b\right|,n)$, 我们便成功获得了一个 $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...

一般应用没啥问题(吧)...

参考资料

Pollard's Rho Algorithm - Wikipedia
Miller-Rabin primality test - Wikipedia
Birthday problem - Wikipedia
算法导论(第三版), 机械工业出版社; 31.8 素数的测试, 31.9 整数的因子分解

posted @ 2019-03-03 19:48 rvalue 阅读(908) 评论(3) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

rvalue

Τhіnk οncе, сοdе ΝаΝcе;