Miller Rabin素数测试和Pollard Rho算法

翻了好多博客和题解，感觉都讲得不是很清晰qwq，很多地方就一个显然轻飘飘地带过，自己想了好久才想通。

目录

• $Miller\ Rabin$素性测试
  • 前置定理
    • 费马小定理
    • 二次探测定理
  • 算法分析
  • Code View
• $Pollard\ Rho$算法
  • 算法分析
  • Code View

$Miller\ Rabin$素性测试

$MillerRabin$算法是一种高效的单个质数判定方法。虽然是一种不确定的质数判断法，但是在选择多种底数的情况下，正确率是可以接受的。它可以判定的数字范围较大，速度也比较优秀，所以是一种比较实用的算法。

前置定理

费马小定理

若$p$是质数，则$a^p\equiv a\mod p$，如果$a$不是$p$的倍数，还可以写成$a^{p-1}\equiv 1\mod p$，这种写法更常见一些。

二次探测定理

若$p$是素数且$x^2\equiv 1\mod p$,则满足$x\equiv 1\mod p$或$x\equiv p-1\mod p$

证明：

因为$x^2\equiv 1\mod p$，所以$x^2-1\equiv 0\mod p$，则$p|(x+1)(x-1)$

又因为$p$是质数，所以$(x+1)$或$(x-1)$有因子$p$，则$p|(x+1)$或$p|(x-1)$

$x+1\equiv p \mod p$或$x-1\equiv p \mod p$

$x\equiv p-1 \mod p$或$x\equiv 1 \mod p$

$Q.E.D.$

算法分析

设我们要判定的数为$x$，我们用一个素数$p$来进行判定。

首先，如果$x==p$，那么$x$是素数；如果$p|x$，那么$x$不是素数。可以特判掉。

然后，先用费马小定理进行测试（这一步也叫做费马测试），如果$p^{x-1}\%x!=1$，那么$x$不是质数。

否则，我们用二次探测定理进行测试。

令$k=x-1$,如果$2|k$，显然有$(p^{\frac k 2})^2\%x==1$,因为这个式子等价于$p^{x-1}\%x==1$，就是费马小定理，刚才已经判断过了。

令$t=p^{\frac k 2}\%x$，根据二次探测定理，如果$t!=1||t!=x-1$，那么$x$不是素数。

如果$t==1$，那么把$p^{\frac k 2}\%x==1$看作是新的一个条件，如果$2 | k$,将$k/2$,继续重复刚才的内容，判定$p^{\frac k 4}$；当然，如果$k$已经是奇数，那么无法继续判定，所以认定$x$是素数。

如果$t==x-1$，不符合二次探测定理的那个条件式，那么就没有办法继续判定，所以认定$x$是素数。

以上就是$MillerRabin$的算法流程了。

事实上存在很少一部分强伪素数是没有办法被$MillerRabin$算法筛掉，所以可以多选几个底数$p$进行判定，它能逃脱所有底数的筛选的概率很小，正确率是在可接受范围内的。

经过大佬的经验传授，如果$x<=10^{12}$，$p$取$\{2,13,23,1662803\}$就可以。

如果$x<=10^{18}$$p$取$\{2,3,5,7,11,13,17,19,23\}$就可以。

Code View

const int P[]={2,3,5,7,11,13,17,19,23},pn=9;
int ksm(int a,int b,int MOD)
{
	int res=1;
	while(b)
	{
		if(b&1) res=1ll*res*a%MOD;
		a=1ll*a*a%MOD;
		b>>=1;
	}
	return res;
}
bool check(int x,int p)
{
	if(ksm(p%x,x-1,x)!=1) return 0;
	int k=x-1,t;
	while(!(k&1))
	{
		k>>=1;
		t=ksm(p%x,k,x);
		if(t!=1&&t!=x-1) return 0;
		if(t==x-1) return 1;//不符合二次探测的条件式 没有办法继续判定 
	}
	return 1;//k变成了奇数 仍然没有筛出来 
}
bool Miller(int x)
{
	for(int i=1;i<=pn;i++)
	{
		if(x==P[i]) return 1;
		if(x%P[i]==0) return 0;
		if(!check(x,P[i])) return 0;
	}
	return 1;
} 

在快速幂会$T$的情况下，可以先把$2$提出来，然后逆着做，倒着乘过去。

$Pollard\ Rho$算法

$Pollard\ Rho$算法是一个大数质因数分解算法，它的实现是基于$Miller\ Rabin$素性测试。它是一种比较玄学的随机化算法，《算法导论》给出的时间复杂度是$O(\sqrt p)$的，$p$是$n$的一个最小因子。

算法分析

设我们要分解的数是$x$。

首先，我们用$Miller\ Rabin$判断一下$x$是否为质数，如果是，那么就可以统计信息然后返回。

接下来，我们考虑如果可以找到一个数$y$使得$1<gcd(x,y)<x$，$gcd(x,y)$就是$x$的一个非平凡因子，然后可以把$x$分成$gcd(x,y)$和$\frac x{gcd(x,y)}$两部分进行递归计算。

然后考虑怎么求这个$y$。首先随机化一个数$v_0∈[0,x-1]$，然后生成一个序列$v_i=f(v_{i-1})%x$，其中$f(n)$是一个伪随机数函数，例如$f(n)=n^2+t$,$t$是常数。

$v[]$是会形成一个环的，环的最长长度是$x$，因为最长到第$x+1$个数时就会重复，而这个映射关系是唯一的，所以就会成环。

根据生日悖论，环长的期望是$\sqrt x$，所以复杂度可以保证。不过这个结论我不知道怎么证它，所以大概可以用一个期望$dp$来验证的方法来说明它是对的？

定义$f[i]$表示序列已经排到了第$i$个数，$[1,i-1]$中的数都互不相同形成的环长的期望。

那么$f[x+1]=x,f[i]=\frac {i-1} x*\frac{(1+(i-1))*i}{2*i}+\frac{x-i+1}{x}*f[i+1]$

前面是和$[1,i-1]$的数一样产生的期望，和$v[i-1]$一样环长是$1$,和$v[i-2]$一样环长是$2$...概率是一样的，所以求一个平均值就是期望；后面是和前面的数都不一样的产生的期望。

我们可以写出代码：

double f[N];
void work()
{
	int x=rd();
	f[x+1]=x;
	for(int i=x;i>=1;i--)
		f[i]=1.0*(i-1.0)/x*i/2+1.0*(x-i+1)/x*f[i+1];
	printf("%.9f\n%.9f",f[1],f[1]*f[1]);
}

实际结果的话，比$\sqrt x$小，$\sqrt x$大概是结果的$2$到$3$倍。

---

对于求出来的$v_i$，算出$d=gcd(|v_i-v_0|,x)$，并判断$1<d<x$，如果满足就记录$d$并继续递归计算。

如果已经成环了，就没用了，就分解失败，可以调整$t$的值并重新分解。

关于如何探测环的出现，一个稍微有点暴力的方法是每次都存下来，判断是否有出现过，但这样做在数据范围比较大的时候会$MLE$。一种比较有意思的做法是$Floyd$提出来的（怎么又是他）。在一个很长的圆形轨道上行走，要判断什么时候至少走了一圈，我们可以让$A,B$同时从同一起点开始走，$B$的速度是$A$的两倍，当$B$第一次赶上$A$的时候，$B$就走了至少一圈了（准确地来说，应该是2圈？$t=\frac L v,S_B=2v*t=2L$）。

优化：我们每求出一个$v_i$就求了一个$gcd$，$gcd$求得太过频繁，我们完全可以把很多个数累在一起求$gcd$，不会影响正确性。因为如果$gcd(p,x)>1$，那么$gcd(p*q,x)>1$

Code View

洛谷板题：P4718 【模板】Pollard-Rho算法

#include<cstdio>
#include<algorithm>
#include<queue>
#include<cstring>
#include<iostream>
#include<cstdio>
#include<cmath>
#include<map>
using namespace std;
#define INF 0x3f3f3f3f
#define LL long long
LL rd()
{
	LL x=0,f=1;char c=getchar();
	while(c<'0'||c>'9'){if(c=='-')f=-1; c=getchar();}
	while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+(c^48); c=getchar();}
	return f*x;
}
LL ans;
LL Abs(LL x)
{
	if(x>=0) return x;
	return -x;
}
LL gcd(LL a,LL b)
{
	if(b==0) return a;
	return gcd(b,a%b);
}
LL qmul(LL x,LL y,LL MOD)
{//快速乘 
	return (x*y-(long long) ((long double) x/MOD*y)*MOD+MOD)%MOD;
}
const int P[]={0,2,3,5,7,11,13,17,19,23},pn=9;
int ksm(int a,int b,int MOD)
{
	int res=1;
	while(b)
	{
		if(b&1) res=1ll*res*a%MOD;
		a=1ll*a*a%MOD;
		b>>=1;
	}
	return res;
}
bool check(LL x,LL p)
{
	if(ksm(p%x,x-1,x)!=1) return 0;
	LL k=x-1,t;
	while(!(k&1))
	{
		k>>=1;
		t=ksm(p%x,k,x);
		if(t!=1&&t!=x-1) return 0;
		if(t==x-1) return 1;
	}
	return 1;
}
bool Miller(LL x)
{
	for(int i=1;i<=pn;i++)
	{
		if(x==P[i]) return 1;
		if(x%P[i]==0) return 0;
		if(!check(x,P[i])) return 0;
	}
	return 1;
}
void Rho(LL x)
{//
	if(Miller(x))
	{
		ans=max(x,ans);
		return ;
	}
	LL t1=rand()%(x-1)+1;
	LL t2=t1,b=rand()%(x-1)+1;
	t2=(qmul(t2,t2,x)+b)%x;
	LL p=1,i=0;
	while(t1!=t2)
	{
		i++;
		p=qmul(p,Abs(t1-t2),x);
		if(p==0)
		{
			LL t=gcd(Abs(t1-t2),x);
			if(t!=1&&t!=x)
			{
				Rho(t);
				Rho(x/t);
			}
			return ;
		}
		if(i%127==0)//为什么是127...玄学
		{
			p=gcd(p,x); 
			if(p!=1&&p!=x)
			{
				Rho(p);
				Rho(x/p);
				return ;
			}
			p=1;
		} 
		t1=(qmul(t1,t1,x)+b)%x;
		t2=(qmul(t2,t2,x)+b)%x;
		t2=(qmul(t2,t2,x)+b)%x;
	}
	p=gcd(p,x);
	if(p!=1&&p!=x)
	{
		Rho(p);
		Rho(x/p);
		return ;
	}
}
int main()
{
	int T=rd();
	while(T--)
	{
		LL x=rd();
		if(Miller(x))
		{
			puts("Prime");
			continue;
		}
		ans=0;
		while(ans==0)
			Rho(x);
		printf("%lld\n",ans); 
	}
    return 0;
}

（MR还好,PR就是真的脑壳大