Pollard_Rho
Pollard_Rho
\(Pollard Rho \)(在此简称PR)可以用来在 \(O(N^{\frac{1}{4}})\) 的时间内分解质因数.
(这个算法是\(Pollard\)提出来的;算法中会涉及到一个环,它的形状为\(''\rho''\),所以叫\(Pollard Rho\) )
题面
求一个数的最大质因数.
这题不需要卡常,不需要卡常,不需要卡常!!!
前置知识:
Miller_Rabin,快速幂,快速乘,gcd
1.Miller_Rabin
Miller_Rabin(在此简称MR)是一个高效(\(O(log_2 N)\))判素数的随机算法,是PR的基础,十分重要
而MR也有两个前置知识:费马小定理,二次探测定理
(1)费马小定理
这个应该比较简单吧.
若 \(p\) 为质数,有
我们现在要验证的数为\(p\) ,那么可以选取一些数作为\(a\),带入上式,分别检验是否满足费马小定理.
只要有一个\(a\) 不满足 \(a^p \equiv a\pmod{p}\) ,那么p为合数.
但是这只是必要不充分条件.存在这样一类强伪素数\(p\),满足
这也就意味着无法使用费马小定理将它判断为合数.
(2)二次探测
若\(x^2\equiv1\pmod{p}\),若\(p\)为质数,则\(x\equiv1\pmod{p}\)或\(x^2\equiv p-1\pmod{p}\)
证明:
因为\(x^2\equiv1\pmod{p}\)
所以\(p\mid(x-1)(x+1)\)
所以$p\mid x-1 \quad or\quad p\mid x+1 $
所以\(x\equiv1\pmod{p}\)或\(x\equiv p-1\pmod{p}\)
证毕.
那么要如何利用它呢?
我们要检测的数仍然为\(p\),然后再选定一个数\(x\).
此时\(p\)已经满足了费马小定理(否则会被直接判合数),即:\(x^{p-1} \equiv 1\pmod{p}\)
(注意,这个式子中的模数在递归求解的过程中是始终不变的,但指数会改变)
这个式子的形式是不是和二次探测定理中的形式很相似?
实际上,我们可以利用二次探测来继续判断质数.
若\(2\mid p-1\) ,则\((x^\frac{p-1}{2})^2\equiv1\pmod{p}\),判断\(x^\frac{p-1}{2}\) 模\(p\)意义下的值.
a. 若既不是1,也不是p-1,那么说明p是合数,返回false.
b. 若是1,则继续递归 \(x^\frac{p-1}{2}\)
c.若为p-1,那么不能利用二次探测继续递归,说明目前无法验证p为合数,返回true.
多取几个x来测试就可以了.
Code:
IL int qpow(int x,int a,int mod)//快速幂
{
int ans = 1;
while(a)
{
if(a&1) ans=ans*x%mod;
a>>=1,x=x*x%mod;
}
return ans;
}
IL bool test(int x,int p)//费马小定理
{
return qpow(x,p-1,p)==1;
}
IL bool Miller_Rabin(int x,int p)
{
if(!test(x,p)) return false;
int k=p-1;
while(!(k&1))
{
//k/=2;
k>>=1;
RG int t=qpow(x,k,p);
if(t!=1&&t!=p-1) return false;
if(t==p-1) return true;
}
return true;
}
IL bool Test(int p)
{
if(p==1) return false;
if(p==2||p==3||p==5||p==7||p==11) return true;
return Miller_Rabin(2,p)&&Miller_Rabin(3,p)&&Miller_Rabin(5,p)&&Miller_Rabin(7,p);
//取不同的x来测试,提高正确率
}
2.快速幂,快速乘
前者就不说了,主要是快速乘.(不过一般快速幂里面的乘法也要用到快速乘)
它的用处是计算两个longlong级别的数的模意义下的乘法
原理:
网上搜\(O(1)\)快速乘,这里不解释了.
Code:
IL int qmul(int x,int y,int mod)
{
return (x*y-(long long)((long double)x/mod*y)*mod+mod)%mod;
}
3.gcd
我用的是辗转相除,但听说用二进制更快?但是好像也只是常数级别的优化(后文会提到一个很重要的东东)
IL int gcd(int x,int y)
{
return y?gcd(y,x%y):x;
}
了解了上述知识后,你就可以开始做这题了QAQ
做法
方法一:试除法
这个很简单,就不说了.
方法二:rand
我要分解的数为\(N\),我在区间\([1,N]\)中rand一个数,看它是不是N的因数.
(很搞笑吧)
方法三:Birthday Trick 优化的rand(正解)
嗯,我们从\([1,N]\)中rand两个数,那么它们的差为N的因数的可能性会大大提升.
但是因为要两两作差,所以没有优化.
但是我们可以退而求其次,使这两个数的差不一定要满足与N的因数完全相等,而是它们的最大公因数不等于一.
那么这个时候我们成功的几率就很高了.
原因大概如下:
若一个数\(N=p*q\) (p,q均为质数)
那么满足\((N,x) \ne 1(x<N)\)的个数是\(p+q-2\).
所以找一次找出\(p,q\) 成功的概率为\(\frac{p+q-1}{N}\).
推广到求\(N=\prod_{i=1}^{n}prime_i^{a_i}\)的非平凡因子(不是1与本身的因数),找\(i\)次能找到的概率.
\(p_i=1-\prod_{j=0}^i \frac{\phi(N)-j}{N}\)
-----口胡结束-----
如果这样做,据说要选\(N^{\frac{1}{4}}\)个数,无法存储.
因此,我们必须每次运算只记录两个数.
我们构造一个伪随机数,推荐以下代码
t1=(qmul(t1,t1,x)+b);
就是\(x_i=x_{i-1}^2+C\)(C为常数).
比较\(x\)数列的相邻两项.
但是,会出现环.因为我们的\(x\)数列是模意义下生成的,所以可能\(\exists j\ne i,x_i=x_j\).
那么如何解决呢?
用hash吗?不不不,那太麻烦了.
以下内容为引用
那么,如何探测环的出现呢?
一种方法是记录当前产生过的所有的数$x_1 , x_2 , ... ... x_n \(,并检测是否存在\)x_l = x_n (l < n)$。
在实际过程中,当 n 增长到一定大小时,可能会造成的内存不够用的情况。
另一种方法是由Floyd发明的一个算法,我们举例来说明这个聪明而又有趣的算法。
假设我们在一个很长很长的圆形轨道上行走,我们如何知道我们已经走完了一圈呢?
当然,我们可以像第一种方法那样做,但是更聪明的方法是让 A 和 B 按照 B 的速度是
A 的速度的两倍从同一起点开始往前走,当 B 第一次敢上 A 时(也就是我们常说的套圈) ,
我们便知道,B 已经走了至少一圈了。while ( b != a )
a = f(a); // a runs once
b = f(f(b)); // b runs twice as fast.
p = GCD( | b - a | , N);
就这样,我们可以较快的找出\(N\)的两个非平凡因子\(p,q\).
递归求解,直到本身为素数返回即可.
代码如下:
void Pollard_Rho(int x)
{
if(Test(x))//素数测试
{
Ans=max(x,Ans);
return;
}
int t1=rand()%(x-1)+1;
int t2=t1,b=rand()%(x-1)+1;
t2=(qmul(t2,t2,x)+b)%x;//要用快速乘
int i=0;
while(t1!=t2)
{
++i;
int t=gcd(abs(t1-t2),x);
if(t!=1&&t!=x)
{
Pollard_Rho(t),Pollard_Rho(x/t);
}
t1=(qmul(t1,t1,x)+b)%x;
t2=(qmul(t2,t2,x)+b)%x;
t2=(qmul(t2,t2,x)+b)%x;
}
}
正解优化
我在这里提一个比较重要的优化,加上后基本不需要卡常就可以跑进\(2500ms\).
我们要频繁地求gcd,可不可以更快地求呢?
可以!
我们可以将若干个差值累积乘到一起,再求gcd.这与分别求是等价的.
设我们得到的差值为\(\{a_i\}\)
因为:
\(gcd(\prod_{i=1}^{n}a_i,N)=\prod_{i=1}^n gcd(a_i,N)\).(下文会有解释) ①
又有\(gcd(\prod_{i=1}^{n}a_i,N)=gcd((\prod_{i=1}^{n}a_i)\%N,N)\)
所以\(gcd(\prod_{i=1}^{n}(a_i\%N),N)=\prod_{i=1}^n gcd(a_i,N) \%N\)
这样就可以少求很多gcd.
这里的\(a_i\) 在推导过程中只考虑两两互质(\(gcd(k,N)\)不是完全积性,若不加以说明,①式不恒成立).
试想\(gcd(a_i,a_j) = x(x \ne 1)\),则\(a_i=x*b_i,a_j=x*b_j\)
(i)\(gcd(x,N)=1,\)那么\(gcd(a_i,N)=gcd(b_i,N),gcd(a_j,N)=gcd(b_j,N)\)
所以\(gcd(a_i,N)*gcd(a_j,N)=gcd(b_i,N)*gcd(b_j,N)=gcd(b_i*b_j,N)\)(因为\(gcd(b_i,b_j)=1\))
(ii)\(gcd(x,N)=k(k\ne1),\) 那么\(gcd(a_i,N)*gcd(a_j,N)\) 必定会出现两次\(k\),
即\(k^2\mid gcd(a_i,N)*gcd(a_j,N)\).
而\(gcd(a_i*a_j,N)\)可能只会出现一次\(k\)(如果\(N\)的\(k\)不够的话),
即\(k\mid gcd(a_i*a_j,N)\) 但 \(k^2 \nmid gcd(a_i*a_j,N)\) (当然,也有可能\(k^2 \mid gcd(a_i*a_j,N)\),如果\(N\)的\(k\)足够)
即\(gcd(a_i*a_j,N)\mid gcd(a_i,N)*gcd(a_j,N)\),不影响正确性
而我们分别对\(a_i,a_j\) 与N求公因数再相乘时反而会有重复.
至于累乘多少个差值,这就比较玄学了.这题取的是127(可能是面向数据编程?QAQ)
改进代码如下:
void Pollard_Rho(int x)
{
if(Test(x))
{
Ans=max(x,Ans);
return;
}
int t1=rand()%(x-1)+1;
int t2=t1,b=rand()%(x-1)+1;
t2=(qmul(t2,t2,x)+b)%x;
int p=1;
int i=0;
while(t1!=t2)
{
++i;
p=qmul(p,abs(t1-t2),x);
if(p==0) //①
{
int t=gcd(abs(t1-t2),x);
if(t!=1&&t!=x)
{
Pollard_Rho(t),Pollard_Rho(x/t);
}
return;
}
if(i%127==0)
{
p=gcd(p,x);
if(p!=1&&p!=x)
{
Pollard_Rho(p),Pollard_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)
{
Pollard_Rho(p),Pollard_Rho(x/p);
return;//②
}
}
还是有两个要注意的地方:
①:p为0,说明乘上当前差值后变为了x的倍数,那么当前差值与N的gcd一定为x的因子.
②:\(\rho\)环的长度可能小于127,所以需要在跳出循环时判断.
Code:
#include<bits/stdc++.h>
#define gc getchar
#define R register int
#define RG register
#define int long long
#define IL inline
using namespace std;
int rd()
{
int ans = 0,flag = 1;
char ch = gc();
while((ch>'9'||ch<'0')&&ch!='-') ch=gc();
if(ch == '-') flag=-1;
while(ch>='0'&&ch<='9') ans+=(ans<<2)+(ans<<1)+ch-48,ch=gc();
return flag*ans;
}
int Ans,n;
IL int qpow(int x,int a,int mod)
{
RG int ans = 1;
while(a)
{
if(a&1) ans=ans*x%mod;
a>>=1,x=x*x%mod;
}
return ans;
}
IL bool test(int x,int p)
{
return qpow(x,p-1,p)==1;
}
IL bool Miller_Rabin(int x,int p)
{
if(!test(x,p)) return false;
int k=p-1;
while(!(k&1))
{
k>>=1;
RG int t=qpow(x,k,p);
if(t!=1&&t!=p-1) return false;
if(t==p-1) return true;
}
return true;
}
IL bool Test(int p)
{
if(p==1) return false;
if(p==2||p==3||p==5||p==7||p==11) return true;
return Miller_Rabin(2,p)&&Miller_Rabin(3,p)&&Miller_Rabin(5,p)&&Miller_Rabin(7,p);
}
IL int qmul(int x,int y,int mod)
{
return (x*y-(long long)((long double)x/mod*y)*mod+mod)%mod;
}
inline int gcd(int x,int y)
{
return y?gcd(y,x%y):x;
}
void Pollard_Rho(int x)
{
if(Test(x))
{
Ans=max(x,Ans);
return;
}
int t1=rand()%(x-1)+1;
int t2=t1,b=rand()%(x-1)+1;
t2=(qmul(t2,t2,x)+b)%x;
int p=1;
int i=0;
while(t1!=t2)
{
++i;
p=qmul(p,abs(t1-t2),x);
if(p==0)
{
int t=gcd(abs(t1-t2),x);
if(t!=1&&t!=x)
{
Pollard_Rho(t),Pollard_Rho(x/t);
}
return;
}
if(i%127==0)
{
p=gcd(p,x);
if(p!=1&&p!=x)
{
Pollard_Rho(p),Pollard_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)
{
Pollard_Rho(p),Pollard_Rho(x/p);
return;
}
}
main()
{
//freopen("Divide.in","r",stdin);
//freopen("Divide.out","w",stdout);
ios::sync_with_stdio(false);
cin>>n;
for(R i=1;i<=n;i++)
{
RG int t;
cin>>t;
if(Test(t))
{
cout<<"Prime"<<endl;
continue;
}
Ans=0;
while(Ans==0)
Pollard_Rho(t);
cout<<Ans<<endl;
}
return 0;
}
基本没卡常.
复杂度
个人认为这一部分相对不是很重要,毕竟作为一名\(REIO\) ,知道如何用算法就可以了,感兴趣的同学可以看一下不保证正确
而且说实话,这种随机算法的复杂性我不会分析
为什么是\(O(n^\frac{1}{4})\) 呢?
之前有提到:
若一个数\(N=p*q\) (p,q均为质数)
那么满足\((N,x) \ne 1(x<N)\)的个数是\(p+q-2\).
所以找一次找出\(p,q\) 成功的概率为\(\frac{p+q-1}{N}\).
找i次找到的概率\(P_i=1-\prod_{j=0}^i \frac{(N-p-q+2)-j}{N}\)
我们用\(\sqrt{N}\)来近似$p,q $ (这个上界是比较松的\(p+q\ge2\sqrt{pq}=2\sqrt{N}\))
那么
\(P_i\approx1-\prod_{j=0}^i \frac{(N-2\sqrt{N}+2)-j}{N}\)
\(\approx 1-\prod_{j=0}^{i}\frac{\sqrt{N}-j}{\sqrt{N}}\)(这个还是上界)
令\(t=\sqrt{N}\).
那么\(P_i\approx1-\prod_{j=0}^{i}\frac{t-j}{t}\)
(这个式子很熟悉吧)
而生日悖论告诉我们,当\(i\)取到\(\sqrt{t}\)(也就是\(N^\frac{1}{4}\))时,\(P_i\)就会大于\(\frac{1}{2}\).
所以说复杂度为\(O(N^\frac{1}{4})\).
另一种说法是\(O(\sqrt{p})\)(p为N的较小因数),原因应该也类似.
