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

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

目录

• Miller-Rabin质数测试 & Pollard-Rho质因数分解
  • Miller-Rabin质数测试
    • 一些依赖的定理
      • 费马小定理
      • 二次探测引理
    • 实现以及正确率
      • 代码实现
  • Pollard-Rho质因数分解
    • 生日悖论与生日攻击
    • 主要思想
    • 具体实现
      • 代码片段
  • 参考资料

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

于是就学了一发qaq

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

Miller-Rabin质数测试

一个多项式时间的基于随机的质数测试算法.

玄学.

一些依赖的定理

费马小定理

若 (p in mathbb P) 则 (forall ain mathbb N^+, a^{p-1} equiv 1pmod p)

它的逆命题基本上是真的, 于是我们可以把这个定理的逆命题作为判断质数的依据.

为啥说"基本上"是真的呢? 因为有些鬼畜的合数 (n) 满足 (forall ain mathbb N^+,a^{n-1}equiv 1 pmod n). 这种鬼畜的合数被称为 Carmichael 数. 前三个 Carmichael 数是 (561), (1105) 和 (1729). OEIS数列编号 A002997 注意到它们可以挺小的...所以只用这个定理的逆命题测试有很大的问题.

二次探测引理

这个名字好像不是很通用...但是我们在这里先这么叫好了...

若 $pin mathbb P $ , 则 (1mod p) 不存在非平凡二次剩余.

我们尝试用它来作为判质数的条件, 那么就可以:

若 (forall a^2 equiv 1 pmod p, aequiv 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) 是合数的证据".

实现以及正确率

具体过程大概是

1. 把乱七八糟的trival情况判掉, 比如 (le 2) 啊, 偶数啊啥的.
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 ight) ]

我们只要让这个值小于 (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 ight) < prod_{i=0}^{k-1}expleft(-frac i n ight) = expleft(-sum_{i=0}^{k-1} frac i n ight) =expleft(frac{-k(k-1)}{2n} ight) ]

那么我们只要让不等式右边小于 (frac 1 2) 就好了. 那么我们有:

[expleft(frac{-k(k-1)}{2n} ight) < frac 1 2 ]

两边取对数解一下就有:

[k^2-k>2nln2 ]

又因为 (ln 2) 是个常数, 于是 (k=O(sqrt n)).

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

主要思想

设我们要对 (n) 进行因数分解, (n) 存在一个非平凡因子 (p) 且 (ple n/p).

我们取一个次数至少为 (2) 且包含常数项的多项式函数 (f(x)) (比如 (f(x)=x^2+c), (c) 随机取), 设 (g(x)=f(x)mod n), 用 (g) 来生成伪随机数. 方法是随机选取一个初始值 (r) 不断将新的 (x) 代入 (g(x)) , 也就是 (x_1=g(r), x_2=g(g(r)),x_3=g(g(g(r)))) . 那么我们可以认为数列 (langle x_k angle) 看上去是随机的. 而且 (x_k=g(x_{k-1})).

也就是说这个伪随机数列只和它的上一个输出有关. 从一个相同的初值可以得到相同的序列.

由于 (pmid n), 所以数列 (langle x_k mod p angle) 也满足这个性质, 即 (x_kequiv g(x_{k-1}) pmod p). 因为这个伪随机数生成器的状态只和上一个输出有关, 那么只要生成的新值命中了之前已经生成的值, 那么这个随机数列就会出现环.

由于 (p) 比较小, (langle x_k mod p angle) 有很大概率比 (langle x_k angle) 更早出现环 (只要新生成的 (x) 对 (p) 取模的值命中以前出现的值就会出现环, 而 (p) 不会超过 (sqrt n)). 根据生日悖论, 在大约 (sqrt p) 次随机后就会有很大概率出现. 如果出现这种情况, 那么我们就获得了两个值 (x_a equiv x_b pmod p). 于是它们之间的差是可以被 (p) 整除的. 我们取 (p=gcd (left|x_a-x_b ight|,n)) 即可.

至于这个判环的过程, 我们可以使用神奇的Floyd判环法. 建立两个哨兵 (x_a) 和 (x_b), 每次令 (x_a) 前进一步, 令 (x_b) 前进两步, 如果 (x_a=x_b), 则说明走到了环上.

但是我们并不知道 (p) 是多少, 所以并不能直接判断 (x_a = x_b) 来判断 (langle x_k mod p angle) 是否进入环. 但是根据上文所述, 如果出现了 (x_aequiv x_b pmod p), 那么就一定出环了, 这个时候 (pmid gcd (left|x_a-x_b ight|,n)), 我们便成功获得了一个 (n) 的非平凡因子.

不过有时候我们可能会手气不佳, 抽到的 (c) 和 (r) 生成的 (langle x_k angle) 和 (langle x_k mod p angle) 可能会同时进入环(显然前者不可能比后者进环更快, 最多同时). 不难发现这个环长也是 (sqrt p) 级别的. 但是由于这时候找到的 (x_a) 和 (x_b) 是相等的, 于是并不能带来什么信息. 在这个时候有两种可能, 第一是 (nin mathbb P), 第二是脸黑.

普通的Pollard-Rho会到此为止而不提供任何信息.

期望的时间复杂度是 (O(n^{1/4}log n)) 的, 那个 (log) 是 (gcd) 的时候来的.

具体实现

首先上文所述的过程只能用来找非平凡因子而非质因子.

我们把它和 Miller-Rabin 质数测试结合起来就可以进行质因数分解了. 大致流程如下:

1. 用 Miller-Rabin 判断当前 (n) 是否是质数. 如果是, 丢进答案集合, 跑路.
2. 随机选择一对 (c) 和 (r), 进行上文所述的测试过程.
3. 如果找到了一个非平凡因子 (p), 递归求解 (p) 和 (n/p), 然后跑路.
4. 否则, 承认自己脸黑, 返回第 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 imes 10^{30})...跑路了QAQ...

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

参考资料

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