Pollard_Rho 算法

一、前置知识

Miller_Rabin算法


二、概述

\(Pollard\_Rho\) 是一个用来寻找一个合数的因子的算法。

显然,我们可以用试除法在\(\mathcal O(\sqrt{n})\)的复杂度内完成这一操作,但我们对这个复杂度并不满意,我们需要一个更优秀的算法。


三、核心思想

事实上,我们还有一种奇妙的方法来寻找一个合数的因子:

scanf("%d",&n);
vector<int> fac;
for(int i=1;i<=???;++i){
	int a=rand()%(n-1)+1;
	if(n%a==0) fac.push_back(a);
}

这个代码试图通过随机一些数并判断他们是否是\(n\)的因数,毫无疑问,这是一个非常蠢的方法,因为单次操作找到\(n\)的因数的概率实在是太小了。
这个算法非常的劣质,但我们的\(Pollard_Rho\)正是基于这个算法而来的,它的核心思想就是随机。

四、优化

有一个著名的悖论叫做生日悖论,这个悖论的内容大概就是:
如果现在有\(1000\)个数,在他们中寻找到\(37\)的概率非常小,但如果我们要选2个差恰好是37的数,这个概率就会提高一倍,如果我们选择\(k\)个数,使它们中存在2个数的差恰好是37,这个概率又会大大提升,当\(k=30\)时,概率已经超过\(50%\),当\(k=100\)时,概率已经达到\(99.99%\)

这个悖论告诉我们:寻找一些满足答案的组合会比寻找单个个体的概率更大。
我们可以利用这一点来优化上面的那个算法:

因为\(n\)的因数一定是\(n\)与某个数的最大公约数的约数,所以我们考虑寻找一个\(k\)使\(gcd(n,k)>1\),然后对\(gcd(n,k)\)递归处理。
这时候我们就可以用上生日悖论的思想,不直接寻找\(k\),而是寻找一组数\(a_i\),使他们中存在\(gcd(a_i-a_j,n)>1\)

可以证明,这样的做法大约需要选取\(n^{\frac 14}\)个数字,\(\mathcal O(n^2)\)一一对比的话复杂度依然过高。于是我们考虑并不直接选取这么多数,而是依次生成,依次比较。

五、Pollard_Rho

这时,就有人提出了一个函数来生成这样的随机数:

\[f(x)=(x^2+c)\ mod\ n \]

其中c可以随机生成
于是我们先随机生成一个\(x_1\),再代入函数中得到\(x_2=f(x_1),x_3=f(x_2) \dots\)

这个函数生成的随机数效果很好,但它也有一个问题,就是在有些情况下,\(f\)会在嵌套数次之后,在有一个有限数集中无限循环,它的图象大概是这样(图是盗的)

你会发现这个图形很像\(\rho\)(\(rho\)),\(Pollard \_ Rho\)因此而得名。

那么如何解决这个问题呢:直观的想法是直接记录下经过的每个数,出现循环后直接退出。
但这样做消耗的内存太大了,我们可以用\(Floyd\)的方法来判定:
假设我们有2个用\(Pollard\)函数生成的值\(a\)\(b\),假设每嵌套一次函数都是走一步
那么我们让\(a\)每次走一步,\(b\)每次走两步,当\(a\)\(b\)相等时,说明出现了环

六、倍增积累gcd

上面的算法已经非常精妙,但想要通过这道毒瘤的还是有些困难

在计算时,我们将每次产生的\(abs(a-b)\)相乘并积累下来,最后直接判断这个乘积与\(n\)的gcd。
但是,如果某一时刻积累下来的样本的乘积为0了,为了不让样本丢失,我们直接退出循环进行判断即可。
每次计算是的阈值我们可以设为倍增的,并加一个上限,使用\(128\)是不错的选择。

至此,\(Pollard\_Rho\)的算法流程已经结束,这里给出上面那道例题的代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int T;
ll ans;
inline ll add(ll x,ll y,ll mod){return (x+y>=mod)?x+y-mod:x+y;}
inline ll mul(ll x,ll y,ll mod){
	ll z=(long double)x/mod*y;
	ll ret=(x*y-z*mod)%mod;
	return (ret+mod)%mod;
} 
inline ll ksm(ll x,ll y,ll mod){
	ll ret=1;
	for(;y;x=mul(x,x,mod),y>>=1) if(y&1) ret=mul(ret,x,mod);
	return ret;
}
int pr[11]={2,3,5,7,11,13,17,19,23,29,31};
inline bool Miller_Rabin(ll n){
	if(n<2) return false;
	for(int i=0;i<=10;++i){
		if(n==pr[i]) return true;
		if(n%pr[i]==0) return false;
	}
	ll t=n-1,s=0;
	while(!(t&1)) ++s,t>>=1;
	for(int i=0;i<=10;++i){
		ll a=pr[i];ll cur=ksm(a,t,n),nxt;
		for(int j=0;j<s;++j,cur=nxt){
			nxt=mul(cur,cur,n);
			if(nxt==1&&cur!=1&&cur!=n-1) return false;
		}
		if(cur!=1) return false;
	}return true;
}//Miller_Rabin 
inline ll f(ll x,ll c,ll mod){return add(mul(x,x,mod),mod-c,mod);}//用-c在这道题目中的效率貌似更高 
inline ll gcd(ll x,ll y){return y==0?x:gcd(y,x%y);}
inline ll rad(ll x){
	return 1ll*rand()*rand()%x*rand()%x+1;
}
inline ll Pollard_Rho(ll n){
	if(n==4) return 2;//因为一开始我们就将函数迭代了2次,所以我们无法处理n=4时的答案,需要特判 
	ll x=rad(n-1),c=rad(n-1),y=x;
	x=f(x,c,n),y=f(x,c,n);
	for(int lim=1;x!=y;lim=min(128,lim<<1)){//倍增的阈值 
		ll cur=1,nxt;
		for(int i=1;i<lim;++i,cur=nxt){
			nxt=mul(cur,abs(add(x,n-y,n)),n);
			if(!nxt) break;//当nxt=0时,说明已经出现了符合条件的样本,为避免样本丢失直接退出 
			x=f(x,c,n);y=f(f(y,c,n),c,n);//x每次跳1步,y每次跳2步 
		}
		ll d=gcd(cur,n);if(d!=1) return d;
	}
	return n;
}
inline void resolve(ll n){
	ll d=Pollard_Rho(n);
	while(d==n) d=Pollard_Rho(n);//寻找到一个n的因子 
	if(Miller_Rabin(d)) ans=max(ans,d);
	else resolve(d);//如果d是质数,则直接更新答案,否则继续递归处理 
	if(Miller_Rabin(n/d)) ans=max(ans,n/d);
	else resolve(n/d);
}
inline void work(ll n){
	if(Miller_Rabin(n)){
		puts("Prime");
		return ;
	}
	ans=0;
	resolve(n);
	printf("%lld\n",ans);
}
int main(){
	srand(time(0));
	scanf("%d",&T);
	while(T--){
		ll n;
		scanf("%lld",&n);
		work(n);
	}
	return 0;
}
原文地址:https://www.cnblogs.com/tqxboomzero/p/14040987.html