前言
$Miller-Robbin$ 與 $Pollard Rho$ 雖然都是隨機算法，不過用起來是真的爽。
$Miller Rabin$ 算法是一種高效的質數判斷方法。雖然是一種不確定的質數判斷法，但是在選擇多種底數的情況下，正確率是可以接受的。
$Pollard Rho$ 是一個非常玄學的方式，用於在 $O(n^{1/4})$ 的期望時間復雜度內計算合數$n$的某個非平凡因子。
事實上算法導論給出的是 $O(\sqrt p)$ ， $p$ 是 $n$ 的某個最小因子，滿足 $p$ 與 $\frac{n}{p}$ 互質。
但是這些都是期望，未必符合實際。但事實上 $Pollard Rho$ 算法在實際環境中運行的相當不錯。
注：以上摘自洛谷。
Miller-Robbin
前置芝士
1. 費馬小定理
-
內容：若 \(\varphi(p)=p-1,\,p>1\)，則\(a^{p}\equiv a\pmod{p}\) 或 \(a^{p-1}\equiv 1\pmod{p},\,(a<p)\)
-
證明：戳這里
2. 二次探測定理
-
內容：如果 \(\varphi(p)=p-1,\,p>1,\,p>X\) ，且\(X^2\equiv 1\pmod{p}\)，那么\(X=1\ or\ p-1\)
-
證明：
$\because$ $X^2\equiv 1\pmod{p}$
$\therefore$ $p|X^2-1$
$\therefore$ $p|(X+1)(X-1)$
$\because$ $p$是大於$X$的質數
$\therefore$ $p=X+1\ or\ p\equiv X-1\pmod{p}$，即$X=1\ or\ p-1$。
算法原理
由費馬小定理，我們可以有一個大膽的想法：滿足 $a^{p-1}\equiv 1\pmod{p}$ 的數字 $p$ 是一個質數。
可惜，這樣的猜想是錯誤的，可以舉出大量反例，如：$2^{340}\equiv 1\pmod{341}$，然鵝 $341=31*11$ 。
所以，我們可以取不同的 $a$ 多驗證幾次，不過，$\forall a<561,\,a^{560}\equiv 1\pmod{561}$，然鵝 $561=51*11$ 。
這時，二次探測就有很大的用途了。結合費馬小定理，正確率就相當高了。
這里推薦幾個 $a_i$ 的值： $2，3，5，7，11，61，13，17$。用了這幾個 $a_i$，就連那個被稱為強偽素數的 $46856248255981$ 都能被除去。
主要步驟
-
.將 \(p-1\) 提取出所有 \(2\) 的因數，假設有\(s\) 個。設剩下的部分為 \(d\)（這里提取所有\(2\)的因數，是為了下面應用二次探測定理） 。
-
枚舉一個底數 \(a_i\) 並計算 \(x=a_i^{d}\pmod p\)。
-
令 \(y=x^{2}\pmod p\)，如果沒有出現過 \(p-1\)，那么 \(p\) 未通過二次探測，不是質數。
-
否則，若底數已經足夠，則跳出；否則回到第二步。
簡易代碼
```cpp #define ll long long ll p,a[]={2,3,5,7,11,61,13,17}; inline ll mul(ll a,ll b,ll mod) { ll ans=0; for(ll i=b;i;i>>=1) { if(i&1) ans=(ans+a)%p; a=(a<<1)%p; } return ans%p; } inline ll Pow(ll a,ll b,ll p) { ll ans=1; for(ll i=b;i;i>>=1) { if(i&1) ans=mul(ans,a,p); a=mul(a,a,p); } return ans%p; } bool Miller_Robbin(ll p) { if(p==2) return true; if(p==1 || !(p%2)) return false; ll d=p-1;int s=0; while(!(d&1)) d>>=1,++s; for(int i=0;i<=8 && a[i]大致流程
-
先判斷當前數是否是素數（這里就可應用 \(Miller-Robbin\) ），如果是則直接返回
-
如果不是素數的話，試圖找到當前數的一個因子（可以不是質因子）
-
遞歸對該因子和約去這個因子的另一個因子進行分解
如何找因子
一個一個試肯定是不行的。而這個算法的發明者采取了一種清奇的思路。（即采取隨機化算法）
-
我們假設要找的因子為p
-
隨機取一個 \(x、y\)，不斷調整 \(x\) ，具體的辦法通常是 \(x=x*x+c\)（c是隨機的，也可以自己定）。
-
取 \(p=gcd(y-x,n)\) ，若\(p \in \left(1,n\right)\) ，則找到了一個因子，就返回。
-
如果出現 \(x=y\) 的循環，就說明出現了循環，並不斷在這個環上生成以前生成過一次的數，所以我們必須寫點東西來判環：我們可以用倍增的思想，讓\(y\)記住\(x\)的位置，然后\(x\)再跑當前跑過次數的一倍的次數。這樣不斷讓\(y\)記住\(x\)的位置，x再往下跑，因為倍增所以當\(x\)跑到\(y\)時，已經跑完一個圈。
-
同時最開始設定兩個執行次數\(i=1、k=2\)（即倍增的時候用） ，每次取 \(gcd\) 時 \(++i\) ；如果 \(i==k\) ，則令 \(y=x\) ，並將 \(k\) 翻倍。
完整代碼
#include<cstdio>
#include<ctime>
#include<map>
#include<algorithm>
using namespace std;
#define rg register int
#define V inline void
#define I inline int
#define db double
#define B inline bool
#define ll long long
#define F1(i,a,b) for(rg i=a;i<=b;++i)
#define F3(i,a,b,c) for(rg i=a;i<=b;i+=c)
#define ed putchar('\n')
#define bl putchar(' ')
template<typename TP>V read(TP &x)
{
TP f=1;x=0;register char c=getchar();
for(;c>'9'||c<'0';c=getchar()) if(c=='-') f=-1;
for(;c>='0'&&c<='9';c=getchar()) x=(x<<3)+(x<<1)+(c^48);
x*=f;
}
template<typename TP>V print(TP x)
{
if(x<0) x=-x,putchar('-');
if(x>9) print(x/10);
putchar(x%10+'0');
}
int T;
ll n,ans;
ll a[]={2,3,5,7,11,61,13,17,24251};
template<typename TP>inline ll Gcd(TP a,TP b) {return !b?a:Gcd(b,a%b);}
template<typename TP>inline ll mul(TP a,TP b,TP p)
{
ll ans=0;
for(TP i=b;i;i>>=1)
{
if(i&1) ans=(ans+a)%p;
a=(a<<1)%p;
}
return ans%p;
}
template<typename TP>inline ll Pow(TP a,TP b,TP p)
{
ll ans=1;
for(TP i=b;i;i>>=1)
{
if(i&1) ans=mul(ans,a,p);
a=mul(a,a,p);
}
return ans%p;
}
B Miller_Robbin(ll n)
{
if(n<2) return false;
if(n==2) return true;
if(n%2==0) return false;
ll d=n-1;int s=0;
while(!(d&1)) d>>=1,++s;
for(rg i=0;i<=8 && a[i]<n;++i)
{
ll x=Pow(a[i],d,n),y=0;
F1(j,0,s-1)
{
y=mul(x,x,n);
if(y==1 && x!=1 && x!=(n-1)) return false;
x=y;
}
if(y!=1) return false;
}
return true;
}
inline ll Pollard_Rho(ll n)
{
ll x,y,c,i,k;
while(true)
{
ll x=rand()%(n-2)+1;
ll y=rand()%(n-2)+1;
ll c=rand()%(n-2)+1;
i=0,k=1;
while(++i)
{
x=(mul(x,x,n)+c)%n;
if(x==y) break;
ll d=Gcd(abs(y-x),n);
if(d>1) return d;
if(i==k) {y=x;k<<=1;}
}
}
}
V Find(ll n)
{
if(n==1) return;
if(Miller_Robbin(n)) {ans=max(ans,n);return;}
ll p=n;
while(n<=p) p=Pollard_Rho(p);
Find(p);Find(n/p);
}
int main()
{
read(T);srand(time(0));
while(T--)
{
ans=0;
read(n);Find(n);
if(ans==n) puts("Prime");
else print(ans),ed;
}
return 0;
}
**emmmm \(\cdots\) **
這數據也太毒瘤了吧！！
看來要瘋狂卡常了
優化1（不如叫做卡常？）
蛋定的分析一波，我們發現除了 $Pollard-Rho$ 是 $O(n^{1/4})$ 的期望時間復雜度外， $gcd$ 和龜速乘都是 $O(\log N)$ 的。
雖然這種復雜度已經很優秀了，可對於本題的數據（\(T≤350\) 、 \(1≤n≤10^{18}\)），還是太 \(\cdots\)
所以我們要果斷摒棄這種很 $low$ 的龜速乘，改用一種暴力溢出的快速乘：
-
簡單原理： \(a \times b \ \ mod \ \ p=a \times b−\left \lfloor \frac{a \times b}{p} \right \rfloor\)
-
用
long double來處理這個 \(\left \lfloor \frac{a \times b}{p} \right \rfloor\) -
然后處理一下浮點誤差就可以了。
-
模數較大時可能會出鍋。
-
不過出鍋概率很小 \(\cdots\)
如下
```cpp inline ll mul(ll a,ll b,ll mod) { a%=mod,b%=mod; ll c=(long double)a*b/mod; ll ret=a*b-c*mod; if(ret<0) ret+=mod; else if(ret>=mod) ret-=mod; return ret; } ```實踐證明，戰果輝煌：$6pts -> 94pts$ ！！！
優化2（正解）
雖然關於龜速乘的 $O(\log n)$ 的惡劣影響得到了一定遏制，不過，我還是好想 $AC$ 啊！
通過辦法1 ：打表 \(\cdots\)
正確 $AC$ 姿勢如下：
- 我們發現在 $Pollard-Rho$ 中如果長時間隨機化而得不到結果，$gcd$帶來的 $O(\log n)$ 還是很傷腎的！！那有沒有辦法優化呢？答案是肯定的。-
在生成\(x\)的操作中，龜速乘所模的數就是\(n\)，而要求的就是\(n\)的某一個約數，即現在的模數並不是一個質數。
-
根據取模的性質：如果模數和被模的數都含有一個公約數，那么這次模運算的結果必然也會是這個公約數的倍數。所以如果我們將若干個\((y−x)\) 相乘，因為模數是 \(n\) ，所以如果若干個\((y−x)\)中有一個與\(n\)有公約數，最后的結果定然也會含有這個公約數。
-
所以可以多算幾次\((y−x)\)的乘積再來求\(gcd\) （一般連續算\(127\)次再求一次\(gcd\)）。
-
不過\(\cdots\)，如果在不斷嘗試\(x\)的值時碰上一個環，就可能會還沒算到\(127\)次就跳出這個環了，就無法得出答案；同時，可能\(127\)次計算之后，所有\((y−x)\)的乘積都變成了\(n\)的倍數（即\(\prod_{i=1}^{127} {(y-x)} \equiv 0 \pmod{n}\) ）
-
所以我們可以不僅在每計算\(127\)次之后求\(gcd\)、還要在倍增時（即判環時）求\(gcd\)，這樣既維護了其正確性，又判了環！！
-
完整\(AC\)代碼：
#include<cstdio>
#include<ctime>
#include<algorithm>
using namespace std;
#define rg register int
#define V inline void
#define I inline int
#define db double
#define B inline bool
#define ll long long
#define F1(i,a,b) for(rg i=a;i<=b;++i)
#define F3(i,a,b,c) for(rg i=a;i<=b;i+=c)
#define ed putchar('\n')
#define bl putchar(' ')
template<typename TP>V read(TP &x)
{
TP f=1;x=0;register char c=getchar();
for(;c>'9'||c<'0';c=getchar()) if(c=='-') f=-1;
for(;c>='0'&&c<='9';c=getchar()) x=(x<<3)+(x<<1)+(c^48);
x*=f;
}
template<typename TP>V print(TP x)
{
if(x<0) x=-x,putchar('-');
if(x>9) print(x/10);
putchar(x%10+'0');
}
int T;
ll n,ans;
ll a[]={2,3,5,7,11,61,13,17};
template<typename TP>inline ll Gcd(TP a,TP b) {return !b?a:Gcd(b,a%b);}
inline ll mul(ll a,ll b,ll mod)
{
a%=mod,b%=mod;
ll c=(long double)a*b/mod;
ll ret=a*b-c*mod;
if(ret<0) ret+=mod;
else if(ret>=mod) ret-=mod;
return ret;
}
template<typename TP>inline ll Pow(TP a,TP b,TP p)
{
ll ans=1;
for(TP i=b;i;i>>=1)
{
if(i&1) ans=mul(ans,a,p);
a=mul(a,a,p);
}
return ans%p;
}
B Miller_Robbin(ll n)
{
if(n<2) return false;
if(n==2) return true;
if(n%2==0) return false;
ll d=n-1;int s=0;
while(!(d&1)) d>>=1,++s;
for(rg i=0;i<=8 && a[i]<n;++i)
{
ll x=Pow(a[i],d,n),y=0;
F1(j,0,s-1)
{
y=mul(x,x,n);
if(y==1 && x!=1 && x!=(n-1)) return false;
x=y;
}
if(y!=1) return false;
}
return true;
}
inline ll Pollard_Rho(ll n)
{
while(true)
{
ll x=rand()%(n-2)+1;
ll y=rand()%(n-2)+1;
ll c=rand()%(n-2)+1;
ll i=0,k=1,b=1;
while(++i)
{
x=(mul(x,x,n)+c)%n;
b=mul(b,abs(y-x),n);
if(x==y || !b) break;
if(!(i%127) || i==k)
{
ll d=Gcd(b,n);
if(d>1) return d;
if(i==k) y=x,k<<=1;
}
}
}
}
V Find(ll n)
{
if(n<=ans) return;
if(Miller_Robbin(n)) {ans=max(ans,n);return;}
ll p=Pollard_Rho(n);
while(n%p==0) n/=p;
Find(p),Find(n);
}
int main()
{
read(T);srand(time(0));
while(T--)
{
ans=0;
read(n);Find(n);
if(ans==n) puts("Prime");
else print(ans),ed;
}
return 0;
}