素性測試+PollardRho
素數判定
暴力
本質上是檢查其是否能夠不用其本身進行質因數分解。
直接列舉從區間 \([2,n)\) 的數看其是否整除 \(n\) 即可。但是其實發現我們只要列舉到 \(\sqrt n\) 即可，複雜度 \(O(\sqrt n)\)。
inline bool prime(ll n){
for(int i=2;i*i<=n;++i){
if(n%i==0) return 0;
}
return 1;
}
素性測試
素性測試是能夠在不對給定數進行質因數分解的情況下，測試一個數是否為素數。
而素性測試分為兩種，可以絕對的判斷一個數是否為素數的叫 確定性測試。以一定概率期望判斷其為素數的叫 概率性測試
Fermat 素性測試
使用費馬小定理驗證素數。我們知道對於所有素數 \(n\) 有 \(a^{n-1}=1(mod\ n)\)，那麼只有知道 \(a^{n-1}\not=1(mod\ n)\) ，那麼 \(n\) 必定非質數。
那麼我們不斷選取 \([2,n-1]\) 之間的數 \(a\)，如果一直滿足 \(a^{n-1}=1\)，那麼我們暫且可以認為其為素數。
inline bool Fermat(ll p){ if(p<3) return p==2; for(int i=1;i<=10;++i){ int a=rd()%(n-2)+2; if(qpow(a,p-1,p)!=1) return 0; } return 1; }
但是畢竟是概率性測試，其對於某些偽素數會檢驗通過，所以我們一般不單獨使用它。
這些能通過 Fermat 測試的我們稱之為 費馬偽素數。
MillerRabin 素性測試
首先我們得知道一個定理：
二次探測定理
對於一個素數 \(p\) ，一定有 \(x^2=1(mod\ p)\) 的解為 \(x=1(mod\ p)\) 或 \(x=p-1(mod\ p)\) .
證明簡單，只要開方之後取模意義即可。
因此，我們可以將指數 \(p-1\) 分解為 \(t\times 2^k\) ，我們先算出 \(a^t\) ,然後不停平方。沒平方一次就用二次探測檢驗一次，最後再用費馬小定理檢驗一下，就可以大概率完成一個優秀的素性測試。
其實一般來說我們不選取隨機的 \(a\) ，而是選取一些固定的質數。這裡推薦一組：\(\{2,3,5,7,13,19,61,233\}\)
int ts[8]={2,3,5,7,13,19,61,233};
inline bool MillerRabin(ll p){
if(p<3) return p==2;
if(!(p&1) ) return 0;
ll t=p-1,k=0;
while(!(t&1) ) {t>>=1;++k;}
for(int i=0;i<8;++i){
ll a=qpow(ts[i],t,p);
for(int j=0;j<k;++j){
ll b=mul(a,a,p);
if(b==1 and a!=1 and a!=p-1) return 0;
a=b;
}
if(a!=1) return 0;
}
return 1;
}
其實還可以資料小的直接跑暴力，但是沒必要。
分解質因數
暴力
還是列舉區間 \([2,\sqrt n]\) 內的數，然後隨意弄一下就好了。
std::vector<int> CutbyPrime(ll n){
static std::vector<int>ans;
ans.clear();
for(int i=2;i<=n;++i){
if(n%i==0){
while(n%i==0) n/=i;
ans.push_back(i);
}
}
if(n!=1){
ans.push_back(n);
}
return ans;
}
生日悖論
具體就是說，在一年有 \(365\) 天的情況下，房間中有 \(23\) 個人，至少兩個人生日相同的概率達到及以上 \(50\%\) 。
證明可以看 OI wiki ，我就不證啦。
因此大約就是在一年有 \(n\) 天的情況下，房間中有 \(\sqrt n\) 個人，至少兩個人生日相同的概率達到及以上 \(50\%\) 。
那麼我們選數 \(\sqrt n\) 次得到重複的概率就達到了 \(50\%\) 。
PollardRho
我們發現對於一個非素數的數 \(p\) ，如果有兩個數 \(x1,x2\) 和 \(n\) 的一個非平凡因子 \(q\) 滿足 \(x1=x2(mod\ q)\) 且 \(x1\not=x2(mod\ p)\)，那麼我們可以通過 \(gcd(\lvert x1-x2\rvert,p)\) 求得 \(p\) 的一個非平凡因子。
有了上面的方法，看起來很強，但其實我們根本用不了，所以考慮怎麼弄到這個東西。
我們先用一個偽隨機函式 \(F(x)=x^2+t\) （\(t\) 為常數） 生成序列 \(f(x)\)，並定義 \(g(x)=f(x)\%q\)。我們由生日悖論可知，\(g(x)\) 期望在 \(O(\sqrt q)\le O(n^{\frac{1}{4} })\) 步內就會出現重複的數，從而出現環（所以形如 \(\rho\) 以此得名）。我們讓兩個指標錯位的在 \(f(x)\) 上跳，然後逐個檢驗，只要出現一個 gcd 不為 \(1\) 的就返回。
發現 \(gcd\) 帶了一個 \(\log\) ，特別難受。但是我們知道 \(gcd(a,p)>1\ or\ gcd(b,p)>1\Leftrightarrow gcd(ab\%p,p)>1\)
，所以我們把一連串的要測試的值乘起來並對 \(p\) 取模， 達到一個上界（一般取 \(128\)）就算一次 gcd，如果滿足大於 \(1\) 的條件，那麼就對裡面元素掃一遍求一下 gcd 找到其中一個可行解即可。
我們有兩種錯位方法：
1.追及法
設定一個指標步長為 \(1\)，另一個為 \(2\) 即可。
2.倍增法
設定一個指標只取到 \(f_{2^i}\) ，另一個在 \(f_{(2^i,2^{i+1}]}\) 上跳，一個一個判斷。
我們給一個倍增法的模板吧。
inline ll PollardRho(ll p){
ll c=rd()%(p-1)+1;
ll x,y=0,buf=1;
int top=0;
for(int st=1;;st<<=1){
x=y;
for(int i=0;i<st;++i){
y=inc(mul(y,y,p),c,p);
buf=mul(dec(y,x,p),buf,p);
sav[++top]=dec(y,x,p);
if(top==lim){
if(gcd(buf,n)>1) break;
top=0;
}
}
if(top==lim) break;
}
for(int i=1;i<=top;++i){
buf=gcd(sav[i],p);
if(buf>1) return buf;
}
return n;
}
我們每找到一個非平凡因子，就可以繼續分治 \(q,n/q\) 。總之就是這樣了。
給一個參考程式碼：
#include<bits/stdc++.h>
#define ll long long
#define db double
#define filein(a) freopen(#a".in","r",stdin)
#define fileot(a) freopen(#a".out","w",stdout)
#define sky fflush(stdout);
#define gc getchar
#define pc putchar
namespace IO{
inline bool blank(const char &c){
return c==' ' or c=='\n' or c=='\t' or c=='\r' or c==EOF;
}
inline void gs(char *s){
char ch=gc();
while(blank(ch) ) {ch=gc();}
while(!blank(ch) ) {*s++=ch;ch=gc();}
*s=0;
}
inline void gs(std::string &s){
char ch=gc();s+='#';
while(blank(ch) ) {ch=gc();}
while(!blank(ch) ) {s+=ch;ch=gc();}
}
inline void ps(char *s){
while(*s!=0) pc(*s++);
}
inline void ps(const std::string &s){
for(auto it:s)
if(it!='#') pc(it);
}
template<class T>
inline void read(T &s){
s=0;char ch=gc();bool f=0;
while(ch<'0'||'9'<ch) {if(ch=='-') f=1;ch=gc();}
while('0'<=ch&&ch<='9') {s=s*10+(ch^48);ch=gc();}
if(ch=='.'){
db p=0.1;ch=gc();
while('0'<=ch&&ch<='9') {s=s+p*(ch^48);p*=0.1;ch=gc();}
}
s=f?-s:s;
}
template<class T,class ...A>
inline void read(T &s,A &...a){
read(s);read(a...);
}
};
using IO::read;
using IO::gs;
using IO::ps;
int ts[8]={2,3,5,7,13,19,61,233};
inline ll inc(ll a,ll b,ll c){
return (a+=b)>=c?a-c:a;
}
inline ll dec(ll a,ll b,ll c){
return (a-=b)<0?a+c:a;
}
inline ll mul(ll a,ll b,ll c){
return (__int128)a*b%c;
/*
//有時可以考慮龜速乘，但是巨慢，慎用
ll res=0;
while(b){
if(b&1) res=inc(res,a,c);
a=inc(a,a,c);
b>>=1;
}
return res;
*/
}
inline ll qpow(ll a,ll b,ll c){
ll res=1;
while(b){
if(b&1) res=mul(res,a,c);
a=mul(a,a,c);
b>>=1;
}
return res;
}
ll gcd(ll x,ll y){
if(!x) return y;
return gcd(y%x,x);
}
inline bool MillerRabin(ll p){
if(p<=2) return p==2;
if(!(p&1) ) return 0;
ll t=p-1;int k=0;
while(!(t&1) ) {t>>=1;++k;}
for(int i=0;i<8;++i){
if(p==ts[i]) return 1;
ll a=qpow(ts[i],t,p);
for(int j=0;j<k;++j){
ll b=mul(a,a,p);
if(b==1 and a!=1 and a!=p-1) return 0;
a=b;
}
if(a!=1) return 0;
}
return 1;
}
std::mt19937_64 rd(time(0)*1000+clock() );
inline ll PollardRho(ll p){
static ll sav[128+3];
int top=0;
ll c=rd()%(p-1)+1;
ll x,y=0,buf=1;
for(int st=1;;st<<=1){
x=y;
for(int i=0;i<st;++i){
y=inc(mul(y,y,p),c,p);
buf=mul(buf,dec(y,x,p),p);
sav[++top]=dec(y,x,p);
if(top==128){
if(gcd(buf,p)>1) break;
top=0;buf=1;
}
}
if(top==128) break;
}
for(int i=1;i<=top;++i){
buf=gcd(sav[i],p);
if(buf>1) return buf;
}
return p;
}
ll mx;
inline void divide(ll p){
if(p<=mx or p<2) return;
if(MillerRabin(p) ){
mx=std::max(mx,p);
return;
}
ll q=PollardRho(p);
while(p==q) q=PollardRho(p);
while(p%q==0) p/=q;
divide(p);divide(q);
}
int main(){
filein(a);fileot(a);
int T;read(T);
while(T--){
ll n;read(n);
mx=0;divide(n);
if(mx==n){
printf("Prime\n");
}else printf("%lld\n",mx);
}
return 0;
}