前言
(Pollard Rho)是一个著名的大数质因数分解算法,它的实现基于一个神奇的算法:(MillerRabin)素数测试(关于(MillerRabin),可以参考这篇博客:初学MillerRabin素数测试)。
期望下,(Pollard Rho)算法可以达到极快的复杂度。
核心思想
在(ZJOI2019Day1)讲课期间,它是被(CQZ)神仙作为随机算法内的一部分来进行介绍的。
由此可见,其核心思想便是随机二字。
操作流程
首先,我们先用(MillerRabin)判断当前数(x)是否为质数,若是,则可直接统计信息并退出函数
否则,我们考虑,如果能找到一个数(s)使得(1<gcd(s,x)<x),则(x)必然可以被分成(gcd(s,x))和(frac x{gcd(s,x)})两部分,接下来递归处理这两部分即可。
那么我们要找的,就是一个符合条件的(s)。
而找的过程利用了随机。
我们随机出一个(1sim x-1)范围内的数(v_0),然后,不断计算(v_i=(v_{i-1}*v_{i-1}+t)\%x)((t)为一个自己设定的常数)。
则我们每次判断(d=gcd(abs(v_i-v_0),x))是否大于(1)且不等于(x)本身,若满足则可直接返回(d)。
由于(v_i)最终肯定会形成环,而在形成环的时候,我们就无法继续操作了,因此当(v_i=v_0),我们就可以退出函数并返回(x)本身,表示分解失败。
对于分解失败的情况,我们可以调整(t)并重新进行分解。
根据生日悖论,形成的期望环长为(O(sqrt x)),因此复杂度得以保证。
但是,光这样依然不够,我们还需要加优化。
优化:路径倍长
考虑到如果对于每个(v_i)都要计算一遍(gcd),显然很慢。
于是,就会想到每隔一段时间将这些数一起进行一次(gcd)。
而隔的时间如何确定呢?
我们可以倍增。
用一个变量(s)统计(abs(v_i-v_0))之积并向(n)取模,因为取模不会改变(gcd)。
如果某一时刻(s)变成了(0),则分解失败,退出函数并返回(x)本身。
然后每隔(2^k)个数,我们计算(gcd(s,x))是否大于(1)且小于(x),然后重新把(v_0)设置为当前的(v_i)。
这样一来,复杂度就得到了大大优化。
代码(板子题)
#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define RL Reg LL
#define Con const
#define CI Con int&
#define CL Con LL&
#define I inline
#define W while
#define LL long long
#define Gmax(x,y) (x<(y)&&(x=(y)))
#define abs(x) ((x)<0?-(x):(x))
#define hl_AK_NOI true
using namespace std;
I LL Qmul(CL x,CL y,CL X)//快速乘
{
RL k=(1.0L*x*y)/X,t=x*y-k*X;
t-=X;W(t<0) t+=X;return t;
}
class MillerRabin//MillerRabin判素数板子
{
private:
#define Pcnt 12
Con int P[Pcnt]={2,3,5,7,11,13,17,19,61,2333,4567,24251};
I LL Qpow(RL x,RL y,CL X)
{
RL t=1;W(y) y&1&&(t=Qmul(t,x,X)),x=Qmul(x,x,X),y>>=1;
return t;
}
I bool Check(CL x,CI p)
{
if(!(x%p)||Qpow(p%x,x-1,x)^1) return false;
RL k=x-1,t;W(!(k&1))
{
if((t=Qpow(p%x,k>>=1,x))^1&&t^(x-1)) return false;
if(!(t^(x-1))) return true;
}return true;
}
public:
bool IsPrime(CL x)
{
if(x<2) return false;
for(RI i=0;i^Pcnt;++i) {if(!(x^P[i])) return true;if(!Check(x,P[i])) return false;}
return true;
}
};
class PollardRho//PollardRho分解质因数
{
private:
#define Rand(x) (1LL*rand()*rand()*rand()*rand()%(x)+1)
LL ans;MillerRabin MR;
I LL gcd(CL x,CL y) {return y?gcd(y,x%y):x;}//求gcd
I LL Work(CL x,CI y)//分解
{
RI t=0,k=1;RL v0=Rand(x-1),v=v0,d,s=1;W(hl_AK_NOI)//初始化随机一个v0
{
if(v=(Qmul(v,v,x)+y)%x,s=Qmul(s,abs(v-v0),x),!(v^v0)||!s) return x;//计算当前v,统计乘积,判断是否分解失败
if(++t==k) {if((d=gcd(s,x))^1) return d;v0=v,k<<=1;}//倍增
}
}
I void Resolve(RL x,RI t)//分解
{
if(!(x^1)||x<=ans) return;if(MR.IsPrime(x)) return (void)Gmax(ans,x);//先进行特判
RL y=x;W((y=Work(x,t--))==x);W(!(x%y)) x/=y;Resolve(x,t),Resolve(y,t);//找到一个因数y,然后分解(注意除尽)
}
public:
I PollardRho() {srand(20050521);}//初始化随机种子
I LL GetMax(CL x) {return ans=0,Resolve(x,302627441),ans;}//求答案
}P;
int main()
{
RI Ttot;RL n,res;scanf("%d",&Ttot);
W(Ttot--) scanf("%lld",&n),(res=P.GetMax(n))^n?printf("%lld
",res):puts("Prime");//输出
return 0;
}