翻了好多部落格和題解，感覺都講得不是很清晰qwq，很多地方就一個顯然輕飄飄地帶過，自己想了好久才想通。

目錄

• $Miller\ Rabin$素性測試
  • 前置定理
    • 費馬小定理
    • 二次探測定理
  • 演算法分析
  • Code View
• $Pollard\ Rho$演算法
  • 演算法分析
  • Code View

\(Miller\ Rabin\)素性測試

\(MillerRabin\)演算法是一種高效的單個質數判定方法。雖然是一種不確定的質數判斷法，但是在選擇多種底數的情況下，正確率是可以接受的。它可以判定的數字範圍較大，速度也比較優秀，所以是一種比較實用的演算法。

前置定理

費馬小定理

若\(p\)是質數，則\(a^p\equiv a\mod p\)

二次探測定理

若\(p\)是素數且\(x^2\equiv 1\mod p\),則滿足\(x\equiv 1\mod p\)或\(x\equiv p-1\mod p\)

證明：

因為\(x^2\equiv 1\mod p\)，所以\(x^2-1\equiv 0\mod p\)，則\(p|(x+1)(x-1)\)

又因為\(p\)是質數，所以\((x+1)\)或\((x-1)\)有因子\(p\)，則\(p|(x+1)\)或\(p|(x-1)\)

\(x+1\equiv p \mod p\)

\(x\equiv p-1 \mod p\)或\(x\equiv 1 \mod p\)

\(Q.E.D.\)

演算法分析

設我們要判定的數為\(x\)，我們用一個素數\(p\)來進行判定。

首先，如果\(x==p\)，那麼\(x\)是素數；如果\(p|x\)，那麼\(x\)不是素數。可以特判掉。

然後，先用費馬小定理進行測試（這一步也叫做費馬測試），如果\(p^{x-1}\%x!=1\)，那麼\(x\)不是質數。

否則，我們用二次探測定理進行測試。

令\(k=x-1\),如果\(2|k\)，顯然有\((p^{\frac k 2})^2\%x==1\)

令\(t=p^{\frac k 2}\%x\)，根據二次探測定理，如果\(t!=1||t!=x-1\)，那麼\(x\)不是素數。

如果\(t==1\)，那麼把\(p^{\frac k 2}\%x==1\)看作是新的一個條件，如果\(2 | k\),將\(k/2\),繼續重複剛才的內容，判定\(p^{\frac k 4}\)；當然，如果\(k\)已經是奇數，那麼無法繼續判定，所以認定\(x\)是素數。

如果\(t==x-1\)，不符合二次探測定理的那個條件式，那麼就沒有辦法繼續判定，所以認定\(x\)是素數。

以上就是\(MillerRabin\)的演算法流程了。

事實上存在很少一部分強偽素數是沒有辦法被\(MillerRabin\)演算法篩掉，所以可以多選幾個底數\(p\)進行判定，它能逃脫所有底數的篩選的概率很小，正確率是在可接受範圍內的。

經過大佬的經驗傳授，如果\(x<=10^{12}\)，\(p\)取\(\{2,13,23,1662803\}\)就可以。

如果\(x<=10^{18}\)\(p\)取\(\{2,3,5,7,11,13,17,19,23\}\)就可以。

Code View

const int P[]={2,3,5,7,11,13,17,19,23},pn=9;
int ksm(int a,int b,int MOD)
{
	int res=1;
	while(b)
	{
		if(b&1) res=1ll*res*a%MOD;
		a=1ll*a*a%MOD;
		b>>=1;
	}
	return res;
}
bool check(int x,int p)
{
	if(ksm(p%x,x-1,x)!=1) return 0;
	int k=x-1,t;
	while(!(k&1))
	{
		k>>=1;
		t=ksm(p%x,k,x);
		if(t!=1&&t!=x-1) return 0;
		if(t==x-1) return 1;//不符合二次探測的條件式 沒有辦法繼續判定 
	}
	return 1;//k變成了奇數 仍然沒有篩出來 
}
bool Miller(int x)
{
	for(int i=1;i<=pn;i++)
	{
		if(x==P[i]) return 1;
		if(x%P[i]==0) return 0;
		if(!check(x,P[i])) return 0;
	}
	return 1;
} 

在快速冪會\(T\)的情況下，可以先把\(2\)提出來，然後逆著做，倒著乘過去。

\(Pollard\ Rho\)演算法

\(Pollard\ Rho\)演算法是一個大數質因數分解演算法，它的實現是基於\(Miller\ Rabin\)素性測試。它是一種比較玄學的隨機化演算法，《演算法導論》給出的時間複雜度是\(O(\sqrt p)\)的，\(p\)是\(n\)的一個最小因子。

演算法分析

設我們要分解的數是\(x\)。

首先，我們用\(Miller\ Rabin\)判斷一下\(x\)是否為質數，如果是，那麼就可以統計資訊然後返回。

接下來，我們考慮如果可以找到一個數\(y\)使得\(1<gcd(x,y)<x\)，\(gcd(x,y)\)就是\(x\)的一個非平凡因子，然後可以把\(x\)分成\(gcd(x,y)\)和\(\frac x{gcd(x,y)}\)兩部分進行遞迴計算。

然後考慮怎麼求這個\(y\)。首先隨機化一個數\(v_0∈[0,x-1]\)，然後生成一個序列\(v_i=f(v_{i-1})%x\)，其中\(f(n)\)是一個偽隨機數函式，例如\(f(n)=n^2+t\),\(t\)是常數。

\(v[]\)是會形成一個環的，環的最長長度是\(x\)，因為最長到第\(x+1\)個數時就會重複，而這個對映關係是唯一的，所以就會成環。

根據生日悖論，環長的期望是\(\sqrt x\)，所以複雜度可以保證。不過這個結論我不知道怎麼證它，所以大概可以用一個期望\(dp\)來驗證的方法來說明它是對的？

定義\(f[i]\)表示序列已經排到了第\(i\)個數，\([1,i-1]\)中的數都互不相同形成的環長的期望。

那麼\(f[x+1]=x,f[i]=\frac {i-1} x*\frac{(1+(i-1))*i}{2*i}+\frac{x-i+1}{x}*f[i+1]\)

前面是和\([1,i-1]\)的數一樣產生的期望，和\(v[i-1]\)一樣環長是\(1\),和\(v[i-2]\)一樣環長是\(2\)...概率是一樣的，所以求一個平均值就是期望；後面是和前面的數都不一樣的產生的期望。

我們可以寫出程式碼：

double f[N];
void work()
{
	int x=rd();
	f[x+1]=x;
	for(int i=x;i>=1;i--)
		f[i]=1.0*(i-1.0)/x*i/2+1.0*(x-i+1)/x*f[i+1];
	printf("%.9f\n%.9f",f[1],f[1]*f[1]);
}

實際結果的話，比\(\sqrt x\)小，\(\sqrt x\)大概是結果的\(2\)到\(3\)倍。

對於求出來的\(v_i\)，算出\(d=gcd(|v_i-v_0|,x)\)，並判斷\(1<d<x\)，如果滿足就記錄\(d\)並繼續遞迴計算。

如果已經成環了，就沒用了，就分解失敗，可以調整\(t\)的值並重新分解。

關於如何探測環的出現，一個稍微有點暴力的方法是每次都存下來，判斷是否有出現過，但這樣做在資料範圍比較大的時候會\(MLE\)。一種比較有意思的做法是\(Floyd\)提出來的（怎麼又是他）。在一個很長的圓形軌道上行走，要判斷什麼時候至少走了一圈，我們可以讓\(A,B\)同時從同一起點開始走，\(B\)的速度是\(A\)的兩倍，當\(B\)第一次趕上\(A\)的時候，\(B\)就走了至少一圈了（準確地來說，應該是2圈？\(t=\frac L v,S_B=2v*t=2L\)）。

優化：我們每求出一個\(v_i\)就求了一個\(gcd\)，\(gcd\)求得太過頻繁，我們完全可以把很多個數累在一起求\(gcd\)，不會影響正確性。因為如果\(gcd(p,x)>1\)，那麼\(gcd(p*q,x)>1\)

Code View

洛谷板題：P4718 【模板】Pollard-Rho演算法

#include<cstdio>
#include<algorithm>
#include<queue>
#include<cstring>
#include<iostream>
#include<cstdio>
#include<cmath>
#include<map>
using namespace std;
#define INF 0x3f3f3f3f
#define LL long long
LL rd()
{
	LL x=0,f=1;char c=getchar();
	while(c<'0'||c>'9'){if(c=='-')f=-1; c=getchar();}
	while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+(c^48); c=getchar();}
	return f*x;
}
LL ans;
LL Abs(LL x)
{
	if(x>=0) return x;
	return -x;
}
LL gcd(LL a,LL b)
{
	if(b==0) return a;
	return gcd(b,a%b);
}
LL qmul(LL x,LL y,LL MOD)
{//快速乘 
	return (x*y-(long long) ((long double) x/MOD*y)*MOD+MOD)%MOD;
}
const int P[]={0,2,3,5,7,11,13,17,19,23},pn=9;
int ksm(int a,int b,int MOD)
{
	int res=1;
	while(b)
	{
		if(b&1) res=1ll*res*a%MOD;
		a=1ll*a*a%MOD;
		b>>=1;
	}
	return res;
}
bool check(LL x,LL p)
{
	if(ksm(p%x,x-1,x)!=1) return 0;
	LL k=x-1,t;
	while(!(k&1))
	{
		k>>=1;
		t=ksm(p%x,k,x);
		if(t!=1&&t!=x-1) return 0;
		if(t==x-1) return 1;
	}
	return 1;
}
bool Miller(LL x)
{
	for(int i=1;i<=pn;i++)
	{
		if(x==P[i]) return 1;
		if(x%P[i]==0) return 0;
		if(!check(x,P[i])) return 0;
	}
	return 1;
}
void Rho(LL x)
{//
	if(Miller(x))
	{
		ans=max(x,ans);
		return ;
	}
	LL t1=rand()%(x-1)+1;
	LL t2=t1,b=rand()%(x-1)+1;
	t2=(qmul(t2,t2,x)+b)%x;
	LL p=1,i=0;
	while(t1!=t2)
	{
		i++;
		p=qmul(p,Abs(t1-t2),x);
		if(p==0)
		{
			LL t=gcd(Abs(t1-t2),x);
			if(t!=1&&t!=x)
			{
				Rho(t);
				Rho(x/t);
			}
			return ;
		}
		if(i%127==0)//為什麼是127...玄學
		{
			p=gcd(p,x); 
			if(p!=1&&p!=x)
			{
				Rho(p);
				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)
	{
		Rho(p);
		Rho(x/p);
		return ;
	}
}
int main()
{
	int T=rd();
	while(T--)
	{
		LL x=rd();
		if(Miller(x))
		{
			puts("Prime");
			continue;
		}
		ans=0;
		while(ans==0)
			Rho(x);
		printf("%lld\n",ans); 
	}
    return 0;
}

（MR還好,PR就是真的腦殼大