Miller-Rabin素数测试与Pollard Rho大整数分解

Miller-Rabin素数测试

对一个整数n的朴素的素数测试方法是对$[2,\sqrt{n}]$内的数进行试除以查找是否存在因子，时间复杂度为$O(\sqrt{n})$
而Miller-Rabin素数测试是一个以概率为基础的随机算法，时间复杂度可以达到$O(1)$

前置定理

费马小定理：若$p$为素数，则$a^{p-1}\equiv 1(\mathrm{mod}p)$对所有$a\in Z_p^{\ast }$均成立

其中$Z_p^{\ast }$是与素数$p$互素的所有数的集合

若$p$是一个奇素数，则方程$x^2\equiv 1(\mathrm{mod}{p}^2)$仅有两个解$x=1$和$x=-1$

我们定义若一个数$x$满足方程$x^2\equiv 1(\mathrm{mod}n)$，但$x$不等于1或-1，则$x$是一个以$n$为膜的1的非平凡平方根
则我们根据上述定理可以得到推论

如果对膜$n$存在1的非平凡平方根，则$n$为合数

Miller-Rabin的思想

Miller-Rabin素数测试的想法十分简单，即随机一个整数带入上述两个定理检查n是否不是素数

具体过程如下：
随机准备一个小素数$a$，并令$n-1=2^{t}u$
对$x^2=(a^u{)}^{2^k}(k=1,2,3,...,t)$分别判断是否满足$x^2\equiv 1(\mathrm{mod}n)$的解不为$x=1$和$x-1$，若满足则可判定n不是质数
最后再利用费马小定理判断$n$是否满足$a^{n-1}\equiv 1(\mathrm{mod}n)$，若不满足则判定$n$不是质数

若上述检查都未找到$n$是合数的证据，则$n$通过了一次Miller-Rabin测试
若s次测试后仍未找到可判定n未合数的证据，则判定n未素数

算法导论中证明

对任意奇数$n>2$，进行$s$次Miller-Rabin测试，出错的概率至多为$2^{-s}$

代码实现

bool witness(int a, int n, int u, int t)
{	
	int x = qpow(a, u, n); // a^u mod n
	for(int i=1;i<=t;++i)
	{
		int nx = x*x%n; // (a^u) ^ (2^k)
		if(nx==1 && x!=1 && x!=n-1) return true; // x^2为以n为膜的1的非平凡平方根，判断n是合数
		x = nx;
	}
	if(x!=1) return true; // 此时x = (a^u) ^ (2^t) = n-1，费马测试 x!=-1 则判断n是合数
	else return false;
}

bool millerRabin(int n)
{
	if(n==2) return true;
	if(n<2 || (n&1)==0) return false; // n=0或1或n为偶数
	
	int u=n-1; // n-1=u*2^t
	int t=0;
	while(!(u&1))
	{
		u>>=1;
		t++;
	}
	
	// 进行s次Miller-Rabin素数测试
	for(int i=0;i<s&&prime[i]<n;++i){
		if(witness(prime[i], n, u, t)) return false;
	}
	return true;
}

---

Pollard Rho大整数分解

众所周知对一个大整数进行质因数分解是十分困难的
朴素的$O(\sqrt{n})$的试除法即使使用超级计算机也不可能分解千位的大整数，这也是当今密码学的一个重要基础

Pollard给出的Rho启发式算法是一个基于概率的随机算法
算法导论中证明他能在$O(\sqrt{p})$的期望时间内计算n的某个非平凡因子，其中$p$是某个n的最小的因子，满足$p$与$n/p$互质

Rho随机数

一个最简单的随机算法思路就是随机一个数x判断其是否为n的因子
但实际上，若采用"组合随机采样"的方法，满足答案的组合比单个要多很多，生日悖论便是一个例子

那我们不妨产生一组随机数$x_1,x_2,x_3...,x_n$，并每次随机选择一对数$x_i,x_j$，若$d=gcd(|x_i-x_j|,n)>1$则d是n的一个因子
为此Pollard设计了伪随机数列$x_i=(x_i-1+c)\mathrm{mod}n$，其中c是一个随机的常数

这个伪随机数列显然是有循环节的，以$x_i=(x_i-1-1)\mathrm{mod}n$为例，它的随机序列如下图（取自算法导论）
因为形似$\rho$，所以便取名为Rho启发式方法

Floyd优化的Pollard Rho算法

为了不在rho形伪随机数序列中陷入死循环，可以采用floyd判圈法
令p，q同时从rho的“尾巴”出出发，每次p走一步，q走两步
当$x_p=x_q$时，两者走过的路径差恰好是环长的整数倍，此时便可退出循环

对于循环中每一对$x_p,x_q$都计算$d=gcd(|x_p,x_q|,n)$，若d>1则找到了n的一个因数d
若直到推出循环都未找到满足条件的d，可以重新随机一个c再次进行算法

int f(int x,int c, int n)
{
	return (x*x + c) % n;
}

int pollardRho(int n)
{
	int c = randint(1, n);
	int t1 = randint(0, n-1);
	int t2 = f(t1, c, n);
	
	while(t1!=t2)
	{
		int d = gcd(abs(t1-t2), n);
		if(d!=1 && d!=n) return d;
		t1 = f(t1, c, n);
		t2 = f(f(t2, c, n), c, n);
	}
	return n; // 没有找到因数
}

---

例题

洛谷P4718 【模板】Pollard-Rho算法
题目比较卡常，pollardRho做了一些减少gcd次数的常数优化，这部分可以翻看洛谷题解的解释

#include<iostream>
#include<cstdlib> 
#include<algorithm>
#include<cstring>
#include<cstdio>
using namespace std;
typedef __int128 ll128;
typedef long long lt;
typedef double db;
typedef unsigned long long ull;
#define eps 1e-8
#define abs(x) (x>=0 ?x :-(x))
#define cmax(a, b) (a>=b?a:b)
#define cmin(a, b) (a<=b?a:b)

lt read()
{
    lt f=1,x=0;
    char ss=getchar();
    while(ss<'0'||ss>'9'){if(ss=='-')f=-1;ss=getchar();}
    while(ss>='0'&&ss<='9'){x=x*10+ss-'0';ss=getchar();}
    return x*f;
}

lt prime[11]={2,3,5,7,11,13,19,61,1103,2333,24251};
lt ans;

lt gcd(lt a, lt b){
	return b==0 ?a :gcd(b, a%b);
}

lt qpow(lt a, lt k, lt mod)
{
	lt res=1;
	while(k){
		if(k&1) res=(ll128)res*a%mod;
		a=(ll128)a*a%mod; k>>=1;
	}
	return res;
}

bool witness(lt a, lt n, lt u, int t)
{	
	lt x = qpow(a, u, n);
	for(int i=1;i<=t;++i)
	{
		lt nx = (ll128)x*x%n;
		if(nx==1 && x!=1 && x!=n-1) return true;
		x = nx;
	}
	if(x!=1) return true;
	else return false;
}

bool millerRabin(lt n)
{
	if(n==2) return true;
	if(n==3 || n==7 || n==61 || n==2333 || n==24251) return true;
	if(n<2 || (n&1)==0 || n==46856248255981ll) return false;
	
	lt u=n-1; 
	int t=0;
	while(!(u&1))
	{
		u>>=1;
		t++;
	}
	
	for(int i=0;i<11&&prime[i]<n;++i){
		if(witness(prime[i], n, u, t)) return false;
	}
	return true;
}

lt randint(lt mi, lt mx)
{
	return rand() * (mx-mi) + mi;
}

lt f(lt x,lt c, lt n)
{
	return ((ll128)x*x + c) % n;
}

//lt pollardRho(lt n)
//{
//	lt c = randint(1, n);
//	lt t1 = randint(0, n-1);
//	lt t2 = f(t1, c, n);
//	
//	while(t1!=t2)
//	{
//		lt d = gcd(abs(t1-t2), n);
//		if(d!=1 && d!=n) return d;
//		t1 = f(t1, c, n);
//		t2 = f(f(t2, c, n), c, n);
//	}
//	return n;
//}

lt pollardRho(lt n)
{
	if(n==4) return 2;
	lt c = randint(1, n);
	lt t1 = randint(0, n-1), t2 = t1;
	
	t1 = f(t1, c, n);
	t2 = f(f(t2, c, n), c, n);
	
	for(int lim=1; t1!=t2; lim=cmin(128, lim<<1))
	{
		lt pro = 1;
		for(int i=1; i<=lim; ++i)
		{
			lt tmp = (ll128)pro*abs(t1-t2)%n;
			if(!tmp) break;
			pro = tmp;
			
			t1 = f(t1, c, n);
			t2 = f(f(t2, c, n), c, n);
		}
		lt d = gcd(pro, n);
		if(d!=1) return d;
	}
	return n;
}

void find(lt x)
{
	if(x<=ans || x<2) return;
	if(millerRabin(x))
	{
		ans = cmax(ans, x);
		return;
	}
	
	lt d = x;
	while(x==d){
		d = pollardRho(x);
	}
	
	while(x%d==0){
		x /= d;
	}
	
	find(x); find(d);
}

int main()
{
	int T=read();
	while(T--)
	{
		lt n=read();
		
		ans=0;
		find(n);
		if(ans==n) printf("Prime\n");
		else printf("%lld\n", ans);
	}
	return 0;
}