Pollard Rho 算法及其優化詳解

$ PollardRho $ 算法總結：

$ Pollard~Rho $ 是一個很神奇的算法，用於在 $ O(n^{1/4}) $ 的期望時間復雜度內計算合數n的某個非平凡因子（除了1和它本身以外能整除它的數）。事書上給出的復雜度是 $ O(\sqrt p) $ ， $ p $ 是 $ n $ 的某個最小因子，滿足 $ p $ 與 $ n/p $ 互質。雖然是隨機的，但 $ Pollard~Rho $ 算法在實際環境中運行的相當不錯，不會被卡。

先說一些與 $ PollardRho $ 算法經常一起考的東西：

1. 快速乘
2. Miller Rabin判質數

$ PollardRho $ 算法：

問題模型：

給一個數 $ n $ ，你需要快速求出它的一個非平凡因子。

對於一個比較小的數（ $ n\leq10^9 $ ），我們直接暴力枚舉質因子就行但太大了我們就必須考慮一下隨機算法了。

對於一個非常大的合數 $ n\leq 10^{18} $ （如果可能是質數，我們可以用Miller Rabin判一下）我們要求 $ n $ 的某一個非平凡因子，如果 $ n $ 的質因子很多（就是約數很多）我們也可以暴力隨機弄一弄，但如果是一些（像由兩個差不多的而且很大的質數乘得的 $ n $ ）它的非平凡因子就只有兩個而且 $ n $ 本身還很大，此時如果我們隨機的話復雜度 $ O(n) $ ，這個太難接受了，所以我們想辦法優化一下。

1.我們求 $ gcd $

如果我們直接隨機求 $ n $ 的某一個約數復雜度很高，我們可以考慮求一下 $ gcd $ 。因為我們可以保證一個合數它絕對存在一個質因子小於 $ \sqrt n $ ，所以在 $ n $ 就存在至少 $ \sqrt n $ 個數與 $ n $ 有大於一的公約數。於是我們隨機的可能性就提高到了 $ O(\sqrt n) $ 。

2. 玄學隨機法： $ x^2+c \quad _{x=rand(),c=rand()} $

其實網上大多都叫這種方法為生日悖論（悖論是不符合經驗但符合事實的結論），有興趣的可以去看看。但博主也不懂覺得這個方法跟我們的 $ PollardRho $ 的關聯。在1中我們用gcd的方案，現在我們考慮有沒有辦法能夠快速的找到這些與 $ n $ 有公約數的數（或者說找到一種隨機生成方法，使隨機生成這類我們需要的數的概率變大一點）。

這個隨機生成方法最終被 $ Pollard $ 研發出來了：就是標題那個式子！ 我們定義 $ (y=x\quad x=x\times x+c \quad g=gcd(y-x,p)) $ 為一次操作，我們發現 $ g>1 $ （就是找到我們需要的數）的幾率出奇的高，根據國際上大佬的瘋狂測試，在 $ n $ 中找到一對 $ y-x $ 使兩者有公約數的概率接近 $ n^{\frac{1}{4}} $ ，不得不說太玄學又太優秀了！

判環

這種隨機生成方法雖然優秀，但也有一些需要注意之處，比如有可能會生成一個環，並不斷在這個環上生成以前生成過一次的數，所以我們必須寫點東西來判環：

1. 我們可以讓y根據生成公式以比x快兩倍的速度向后生成，這樣當y再次與x相等時，x一定跑完了一個圈且沒重復跑多少！
2. 我們可以用倍增的思想，讓y記住x的位置，然后x再跑當前跑過次數的一倍的次數。這樣不斷讓y記住x的位置，x再往下跑，因為倍增所以當x跑到y時，已經跑完一個圈，並且也不會多跑太多（這個必須好好想一想，也可以看看代碼如何實現的）（因為這種方法會比第一種用的更多，因為它在 $ PollardRho $ 二次優化時還可以起到第二個作用）（博主這里就給第二種的代碼吧）

inline ll rho(ll n){
    ll x,y,c,g; rg i,j;
    x=y=rand(); c=rand();//初始化
    i=0,j=1;//倍增初始化
    while(++i){
		x=(ksc(x,x,n)+c)%n;//x只跑一次
		if(x==y)break;//跑完了一個環
        g=gcd(Abs(y-x),n);//求gcd（注意絕對值）
   		if(g>1)return g;
		if(i==j)y=x,j<<=1;//倍增的實現
	}
}

$ PollardRho $ 算法的二次優化：

我們發現其實 $ PollardRho $ 算法其實復雜度不止 $ O(n^{\frac{1}{4}}) $ ，它每一次操作都會求一次 $ gcd $ ，所以復雜度會變成 $ O(n^{\frac{1}{4}}\times log) $ 我們發現這個 $ log $ 實在有些拖慢速度（雖然實際上也很快了）。於是一個很奇妙的優化誕生了！

二次優化：我們發現我們在每一次生成操作中，乘法之后所模的數就是我們的 $ n $ ，而我們要求的就是 $ n $ 的某一個約數！也就是說我們現在的模數並不是一個質數！而根據取模的性質：如果模數和被模的數都含有一個公約數，那么這次模運算的結果必然也會是這個公約數的倍數！！！所以如果我們將若干個 $ (y-x) $ 乘到一起，因為中途模的是 $ n $ ，所以如果我的若干個 $ (y-x) $ 中有一個與n有公約數，最后的結果定然也會含有這個公約數！所以我們完全可以多算幾次 $ (y-x) $ 的乘積在來求 $ gcd $ （一般連續算127次再求一次gcd（這個應該有大佬測試過了）），這樣我們的復雜度就能降一個 $ log $ 級別，跑的飛快！

對原算法的一些影響：任何一個優化都是要考慮其對原算法的影響的。這個優化也會有一些影響：首先，因為我們會等大概127次后再去gcd，然而我們有可能在生成時碰上一個環，我們有可能還沒生成127次就跳出這個環了，這樣就無法得出答案；其次，我們的 $ PollardRho $ 算法實在太玄學優秀了，可能我們跑127次之后，所有 $ (y-x) $ 的乘積就變成了n的倍數（模 $ n $ 意義下得到 $ 0 $ ），所以我們不能完全就只呆板的等127次在取模。

那怎么辦呢？（博主在上文就已經提到過了：倍增時我們的 $ i $ 和 $ j $ 會在二次優化時起到另一種作用。）沒錯！就是倍增！我們可以用倍增，分別在生成（1次，2次，4次，8次，16次，32次......）之后進行 $ gcd $ ！這樣就可以完全避免上述的兩個影響，而且我們是倍增來gcd的這對我們的復雜度影響不大！！！！！然后我們看代碼吧！

inline ll rho(ll p){//求出p的非平凡因子
    ll x,y,z,c,g; rg i,j;//先擺出來（z用來存（y-x）的乘積）
    while(1){//保證一定求出一個因子來
        y=x=rand()%p;//隨機初始化
        z=1; c=rand()%p;//初始化
        i=0,j=1;//倍增初始化
        while(++i){//開始玄學生成
            x=(ksc(x,x,p)+c)%p;//可能要用快速乘
            z=ksc(z,Abs(y-x),p);//我們將每一次的（y-x）都累乘起來
            if(x==y||!z)break;//如果跑完了環就再換一組試試（注意當z=0時，繼續下去是沒意義的）
            if(!(i%127)||i==j){//我們不僅在等127次之后gcd我們還會倍增的來gcd
                g=gcd(z,p);
                if(g>1)return g;
                if(i==j)y=x,j<<=1;//維護倍增正確性，並判環（一箭雙雕）
            }
        }
    }
}

然后是兩道道例題，都有點難度的。

$ code1： $

博主給出的樣例代碼是對應多組數據的，對於每個數字檢驗是否是質數，是質數就輸出Prime；如果不是質數，輸出它最大的質因子是哪個。對應洛谷P4718 【模板】Pollard-Rho算法

#include<iostream>
#include<cstdio>
#include<iomanip>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<ctime>
#include<queue>
#include<vector>
#include<stack>
#include<map>
#include<set>

#define ull unsigned long long
#define lb long double
#define ll long long
#define rg register int

using namespace std;

int n;
ll ans,m;

inline ll qr(){//快讀
    char ch;
    while((ch=getchar())<'0'||ch>'9');
    ll res=ch^48;
    while((ch=getchar())>='0'&&ch<='9')
        res=res*10+(ch^48);
    return res;
}

inline ll Abs(ll x){return x<0?-x:x;}//取絕對值
inline ll gcd(ll x,ll y){//非遞歸求gcd
    ll z;
    while(y){z=x;x=y;y=z%y;}
    return x;
}

inline ll ksc(ull x,ull y,ll p){//O(1)快速乘（防爆long long）
    return (x*y-(ull)((lb)x/p*y)*p+p)%p;
}

inline ll ksm(ll x,ll y,ll p){//快速冪
    ll res=1;
    while(y){
        if(y&1)res=ksc(res,x,p);
        x=ksc(x,x,p); y>>=1;
    }return res;
}

inline bool mr(ll x,ll p){//mille rabin判質數
    if(ksm(x,p-1,p)!=1)return 0;//費馬小定理
    ll y=p-1,z;
    while(!(y&1)){//二次探測
        y>>=1; z=ksm(x,y,p);
        if(z!=1&&z!=p-1)return 0;
        if(z==p-1)return 1;
    }return 1;
}

inline bool prime(ll x){ if(x<2)return 0;//mille rabin判質數
    if(x==2||x==3||x==5||x==7||x==43) return 1;
    return mr(2,x)&&mr(3,x)&&mr(5,x)&&mr(7,x)&&mr(43,x);
}

inline ll rho(ll p){//求出p的非平凡因子
    ll x,y,z,c,g; rg i,j;//先擺出來（z用來存（y-x）的乘積）
    while(1){//保證一定求出一個因子來
        y=x=rand()%p;//隨機初始化
        z=1; c=rand()%p;//初始化
        i=0,j=1;//倍增初始化
        while(++i){//開始玄學生成
            x=(ksc(x,x,p)+c)%p;//可能要用快速乘
            z=ksc(z,Abs(y-x),p);//我們將每一次的（y-x）都累乘起來
            if(x==y||!z)break;//如果跑完了環就再換一組試試（注意當z=0時，繼續下去是沒意義的）
            if(!(i%127)||i==j){//我們不僅在等127次之后gcd我們還會倍增的來gcd
                g=gcd(z,p);
                if(g>1)return g;
                if(i==j)y=x,j<<=1;//維護倍增正確性，並判環（一箭雙雕）
            }
        }
    }
}

inline void prho(ll p){//不斷的找他的質因子
    if(p<=ans)return ;//最優性剪紙
    if(prime(p)){ans=p;return;}
    ll pi=rho(p);//我們一次必定能求的出一個因子，所以不用while
    while(p%pi==0)p/=pi;//記得要除盡
    return prho(pi),prho(p);//分開繼續分解質因數
}

int main(){
    n=qr(); srand(time(0));//隨機數生成必備！！！
    for(rg i=1;i<=n;++i){
        ans=1; prho((m=qr()));
        if(ans==m)puts("Prime");
        else printf("%lld\n",ans);
    }
    return 0;
}

$ code2： $

對應：杭電3864 $ D\_num $

根據這道題的答案性質，我們知道這個數在除了1以外還有三個約數，仔細想一下不難發現這類數最多只能有兩個質因子，而這兩個質因子分別就是它的兩個約數，而第三個約數自然就是它自己了！（當然還有可能出現質數的三次方這類特殊的數，所以我們要特判一下）

所以我們的思路就是，先求出 $ n $ 的一個小於它的約數，然后我們可以用這個約數算出 $ n $ 的另一個約數，然后我們判一下它是不是質數的三次方這類特殊的數（是就可以直接輸出答案了），不是的話我們判斷這兩個約數是否分別是不同的質因子，是就能輸出了，不是就說明這不是 $ D\_num $

#include<iostream>
#include<cstdio>
#include<iomanip>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<ctime>
#include<queue>
#include<stack>
#include<vector>
#include<map>
#include<set>

#define ll long long
#define db double
#define rg register int
#define cut {puts("is not a D_num");continue;}//這個可以看一看

using namespace std;

ll n;

inline ll qr(){
    char ch;
    while((ch=getchar())<'0'||ch>'9');
    ll res=ch^48;
    while((ch=getchar())>='0'&&ch<='9')
        res=res*10+(ch^48);
    return res;
}

inline ll Abs(ll x){return x<0?-x:x;}
inline ll gcd(ll x,ll y){
    ll z;
    while(y){z=x,x=y,y=z%y;}
    return x;
}

inline ll ksm(ll x,ll y,ll p){
    ll res=1;
    while(y){
        if(y&1)res=(__int128)res*x%p;
        x=(__int128)x*x%p; y>>=1;
    }return res;
}

inline bool mr(ll x,ll p){
    if(ksm(x,p-1,p)!=1)return 0;
    ll y=p-1,z;
    while(!(y&1)){
        y>>=1; z=ksm(x,y,p);
        if(z!=1&&z!=p-1)return 0;
        if(z==p-1)return 1;
    }return 1;
}

inline bool prime(ll p){ if(p<2)return 0;
    if(p==2||p==3||p==5||p==7||p==43)return 1;
    return mr(2,p)&&mr(3,p)&&mr(5,p)&&mr(7,p)&&mr(43,p);    
}

inline ll rho(ll p){
    ll x,y,z,c,g; rg i,j;
    while(1){
        y=x=rand()%p;
        z=1,c=rand()%p;
        i=0,j=1;
        while(++i){
            x=((__int128)x*x+c)%p;
            z=(__int128)z*Abs(y-x)%p;
            if(x==y)break;
            if(i%72||i==j){
                g=gcd(z,p);
                if(g>1&&g!=p)return g;
                if(i==j)y=x,j<<=1;
            }
        }
    }
}

int main(){
//freopen(".in","r",stdin);
//freopen(".out","w",stdout);
    srand(time(0)); ll p,q;
    while(scanf("%lld",&n)!=EOF){
        if(prime(n)||n<=2)cut;
        p=rho(n); q=n/p;
        if((__int128)p*p==q||(__int128)q*q==p)//這題很卡細節的
        {printf("%lld %lld %lld\n",p,q,n);continue;}
        if(!prime(p)||!prime(q)||p==q)cut;//注意兩個相等也不行
        printf("%lld %lld %lld\n",min(p,q),max(p,q),n);
    }
    return 0;
}