[學習筆記]Pollard-rho算法

一.什么是Pollard-rho

​ 這是一個用來尋找一個合數的因子的算法。很顯然的,我們可以使用試除法，1~\(\sqrt{n}\)之間一個一個試。很顯然他很慢。

二.朴素的代碼

​ 我們來看一個沙雕代碼。

	n=read();
	int a=rand()%(n-2)+2;
	if(n%a==0)cout<<"I found "<<a;
	else cout<<"I'm such a stupid.";

​ 很顯然的，這段代碼很沙雕，完全就是在碰運氣抽簽。在這 𝑁 − 1 個數中，我們恰好只有兩個因子 𝑝 和 𝑞。因此概率就是 \(\frac{n-1}{2}\), 如果𝑁<=\(10^{10}\)，那么成功的機會約為 \(\frac{1}{𝟓𝟎𝟎𝟎𝟎𝟎𝟎𝟎𝟎𝟎}\),這比購買一張彩票中獎的概率還要小，換句話說，我們不得不使用不同的隨機數重復將近 𝑁 次來找到一個因子，顯然，這比試除法還要糟。但是我們可以優化他。

三.Birthday-trick

​ Birthday-trick悖論被稱為生日悖論。先不討論他跟生日的關系是什么。但是他可以對上面的沙雕進行優化。

​ 很顯然的，我們在\(10^{3}\)中，選一個28的概率就很小。但是如果選兩個數的差是28，概率大大滴增加，如果是三個數任意兩個差為28呢？k個數呢？

​ 概率大概如下（盜的圖）

​ 這張圖很直觀的告訴我們在\(10^{3}\)下，k=30時，概率就超過一半，k=100時，美滋滋。

四.優化沙雕

​ 對於沙雕程序，我們我們不考慮選取一個單獨的a,而是考慮選擇兩個數的差看看是否滿足條件。又Birthday-trick我們看出當我們選的足夠多，這個概率就很大了。但是直接除，再加上要選很多的數，這種操作是不行的。再加上答案只有p,q。我們很難找出答案。

​ 但是幸運的是，我們不止這一種方法。對於差為a，我們轉向考慮gcd(a,n)>1.對於這個，出答案的可能性會大大增加。至於有多少種答案我們暫不分析。

五.Pollard-rho

​ 很遺憾，我們去要選取大約\(N^{\frac{1}{4}}\)個數字，當N足夠大的時候，我們的內存扛不住。於是我們采用Pollard-rho算法。

​ 我們不一口氣生成k個數，而是依次生成並進行比較。我們使用一個函數來生成（看起來是的）隨機數。

\[𝑓(𝑥) = (𝑥^2 + 𝑎)\mod\;𝑁 \]

​ ( 我們可以自己指定𝑎，也可以用 rand()隨即函數來生成，這不是重點 )
​ 我們從𝑥1 = 2 或者其他數開始，讓𝑥2 = 𝑓(𝑥1), 𝑥3 = 𝑓(𝑥2), … … 生成規則為\(x_{i+1}=f(x_i)\) .

​ 假設 𝑁 = 55 , 𝑓(𝑥) = ( 𝑥2 + 2 ) 𝑚𝑜𝑑 55

​ 你可以發現對於大部分的數據這個算法能夠正常運行，但是對於某些數據，它將會進入無限循環。為什么呢？這是因為存在 𝑓 環的原因。當它發生的時候，我們會在一個有限數集中進行無限循環
​ 例如，我們可以構造一個偽隨機函數並生成如下偽隨機數：
​ 2, 10, 16, 23, 29, 13, 16, 23, 29, 13 … …
​ 在這個例子中，我們最終將會在16, 23, 29, 13這個圈中無限循環，永遠找不到因子

​ 然后我們大致可以將這個過程畫為一個\(\rho\),這就是算法名稱的來歷。

六.關於環

​ 在上面我們發現了算法是可能進入死循環的，為了判定是否出現了循環，我們使用Floyd的方法（當然可以全部記錄下來一個個找，不過內存不夠）

​ 對於兩個由f產生的數a,b來說，我們假設運行一次\(f()\)是走一步的話，我們讓a一次走一步，b一次走兩步，當b第一次趕上a的時候，就必定出現了環。

七.后記

​ 首先對於這個算法來講，要特判一下因子為2的數，以及使用Miller-rabin可以提前判斷一下是不是素數再進行算法。

八.一道例題

​ luogu的模板題。（超時）

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#include<fstream>
#include<ctime>
using namespace std;
#define m_for(i,a,b) for(int i=(a);i<=(b);++i)
#define ll long long
const int times=10;
ll multi(ll a,ll b,ll m)
{
    ll ans=0%m;
    while(b)
    {
        if(b&1)ans=(ans+a)%m,b--;
        b >>= 1;
        a=(a+a)%m;
    }
    return ans;
}
ll KSM(ll a,ll b,ll m)
{
    ll ans=1%m;
    while(b)
    {
        if(b&1)ans=multi(ans, a, m),b--;
        b>>=1;
        a=multi(a,a,m);
    }
    return ans;
}
inline bool Miller_Rabin(ll p){
	if(p<=2)return p==2;
	if(!(p&1))return 0;
	ll u=p-1;
	ll power=0;
	while(!(u&1))u>>=1,power++;
	m_for(i,1,times){
		ll a=rand()%(p-2)+2,x=KSM(a,u,p),y;
		for(int i=1;i<=power;++i,x=y){
			if((y=multi(x,x,p))==1&&x!=1&&x!=p-1)return 0;
		}
		if(x!=1)return 0;
	}
	return 1;
}
inline ll m_max(ll a,ll b){return a>b?a:b;}
ll ans=0;
ll mod,c;
ll f(ll x){return (multi(x,x,mod)+c)%mod;}
ll gcd(ll a, ll b){return b ? gcd(b, a % b) : a;}
inline ll m_abs(ll x){return x>0?x:-x;}
ll find(ll p){
	if(!(p&1))return 2;
	if(Miller_Rabin(p))return p;
	mod=p;
	ll a,b;
	a=b=rand()%(p-2)+2,c=rand()%(p-2)+2;
	do{
		a=f(a);b=f(f(b));
		ll x=gcd(m_abs(a-b),p);
		if(x>1&&x<p)return x;
	}while(a!=b);
	return find(p);
}
void Pollard_rho(ll x){
	if(x<=1)return;
	if(Miller_Rabin(x)){
		ans=m_max(x,ans);return ;
	}
	ll u=find(x);
	Pollard_rho(u);
	Pollard_rho(x/u);
}
int main(){
srand(time(NULL));
int cas;
ll p;
scanf("%d",&cas);
while(cas--){
	scanf("%lld",&p);
	if(Miller_Rabin(p))printf("Prime\n");
	else{
		ans=0;
		Pollard_rho(p);
		printf("%lld\n",ans);
	}
}
	return 0;
}