素性测试+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;

}

巴特西