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的較小因數),原因應該也類似.