Miller_Rabin 素数判定算法

我们先看一道模板题
LOJ #143 素数判定
也许你会想，这不就是素数判定吗，直接上试除法，就这！
然而，这道题目有\(10^5\)次询问，每次询问的\(x \le 10^{18}\)
普通的\(O(\sqrt x)\)的试除法根本过不去，这是就轮到我们的主角：Miller_Rabin算法登场了

2个重要定理

1.费马小定理
对于任意的质数\(p\)，满足\(a^{p-1} \equiv 1(mod\ p)\)

2.二次探测定理
对于任意的质数\(p\)及整数\(x\)，若满足\(x^2 \equiv 1(mod\ p)\)，则一定有\(x \equiv 1或p-1(mod\ p)\)

算法流程

想必你一定会想用以上2个定理的逆命题来判定质数吧。
遗憾的是，二次探测定理逆命题的反例非常多，而对于费马小定理，也已被证明存在无数多个合数满足费马小定理逆命题，叫做\(Carmichael\)数。
但好消息是，\(Carmichael\)数的密度极小，\(1e8\)范围内只有\(255\)个，并且在进行次数较多的二次探测定理逆否命题判定后，\(Carmichael\)数都会被判定出来。

于是我们考虑通过结合二者来进行判定：
令我们需要判定的\(x=2^s*t+1\)
然后通过多次寻找一个较小的数\(a\)，先计算出\(a^t\ mod\ x\)，在进行\(s\)次平方，每次平方的同时使用二次探测定理进行判定
平方\(s\)次后，最终得到的数\(=(a^t)^{2^s}=2^{x-1}\)，于是我们用费马小定理的逆否命题来进行判定

容易看出，当\(a\)是较小的质数时，判定出错的概率会更小
因为有时间复杂度的限制，所以我们选择使用最小的10个质数作为\(a\)依次进行上述判定，顺带的，我们可以先将这些质数的倍数特判以降低一些常数
这样一来，我们进行了10次费马小定理的判定与数十次二次探测定理，这些结合起来之后，算法的错误概率几乎已经可以忽略不计了，因此\(Miller\_Rubin\)是一个可靠的算法

需要注意的是，因为\(x\)是\(long\ long\)范围的数，因此在进行快速幂时，直接进行乘法可能会爆\(long\ long\)
所以我们选择将乘法用类似快速幂的方式，使用\(log\)次加法来实现
于是我们就得到了一个十分优秀，足以通过本题的素数判定算法

code

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll x;
ll pr[20]={2,3,5,7,11,13,17,19,23,29};
inline ll add(ll x,ll y,ll mod){return (x+y>=mod)?x+y-mod:x+y;}
inline ll mul(ll a,ll b,ll mod){
	ll ret=0;
	for(;b;a=add(a,a,mod),b>>=1) if(b&1) ret=add(ret,a,mod);
	return ret;
}
inline ll fsp(ll a,ll b,ll mod){
	ll ret=1;
	for(;b;a=mul(a,a,mod),b>>=1) if(b&1) ret=mul(ret,a,mod);
	return ret;
}
inline bool Miller_Rabin(ll x){
	if(x<2) return false;
	for(int i=0;i<10;++i){
		if(x==pr[i]) return true;
		if(x%pr[i]==0) return false;
	}
	int s=0;ll t=x-1;
	while(!(t&1)){s++;t>>=1;}
	for(int i=0;i<10;++i){
		ll a=pr[i];
		ll cur=fsp(a,t,x),nxt;
		for(int j=0;j<s;++j,cur=nxt){
			nxt=mul(cur,cur,x);
			if(cur!=1&&cur!=x-1&&nxt==1) return false;
		}
		if(cur!=1) return false;
	}
	return true;
}
int main(){
	while(scanf("%lld",&x)!=EOF){
		if(Miller_Rabin(x)) puts("Y");
		else puts("N");
	}
	return 0;
}