Miller_rabin 素数测试
一种用来判断素数的算法。
前置芝士
威尔逊定理
若 (p) 为素数，((p-1)! equiv -1 (mod p))。
证明：
充分性证明：
如果 (p) 不是素数，那么他的因数必定存在于$ 1,2,3,dots,p−1$ 之中，所以 (gcd((p-1)!,p))，那么 ((p-1)! notequiv -1)。
必要性证明：
首先，我们知道 $$p-1equiv -1 (mod p)$$
那么我们只需要证明 ((p-2)! equiv 1 (mod p)) 就可以了。
设 (A={2,3dots,p-2})
对于 (xin A)，肯定存在一个 (a in A)，使得 (ijequiv 1(mod p))
当 (x=a) 时，考虑二次剩余 (x^2 equiv 1(mod p))。
移项就可以得到 ((x+1)(x-1) equiv 0 (mod p))，
那么只有两个解，(x equiv 1(mod p)),(x equiv p-1(mod p))，不成立。
所以其他情况都可以找到对应的 (a)。
所以 ((p-2)!equiv 1(mod p))，((p-1)!equiv p-1(mod p))。
费马小定理
若 (pinmathbb{P},gcd(a,p)=1)，则 (a^{p-1}equiv 1(mod p))，
证明：
因为 (1,2,3,dots ,p-1) 是 (p) 的完全剩余系，(a,p) 互质，
所以 (a,2* a,3* a,4* a, dots (p-1)* a) 也是 (mod p) 的完全剩余系。
所以 (1 * 2 * 3 * dots * (p-1) equiv a * 2a * 3a * dots * (p-1)a (mod p))。
就是 ((p-1)! equiv (p-1)!a^{p-1} (mod p))。
由威尔逊定理得出，((p-1)! equiv -1(mod p))，
两边同时约去 ((p-1)!)，
所以可以得出 (a^{p-1} equiv 1(mod p))。
二次探测定理
若 (p) 为素数，(x^2 equiv 1(mod p))，那么(x equiv pm 1(mod p))。
证明：
原式移项就可以得到 ((x+1)(x-1) equiv 0 (mod p))，
那么只有两个解，(x equiv 1(mod p)),(x equiv p-1(mod p))。
但是这些又和 Miller_rabin 有什么关系呢？
我们把 (p-1) 分解为 (2^k* t)，当 (p) 是素数时，则有 (a^{2^k* t} equiv 1(mod p))。
然后随机出一个数 (a)，可以算出 (a^t)，然后再自乘，就可以得到 (a^{2^k* t}) ，也就是 (a^{p-1})。
但如果我们在自乘之后 (equiv 1(mod p))，而之前 (notequiv 1(mod p))，那么就违背了二次探测定理，则 (p) 不是素数。
else
(Zwad) 告诉我，若 (p) 通过一次测试，则p不是素数的概率仅为25%。
那么用 (2,325,9375,28178,450775,9780504,1795265022) 这几个数进行判断,在 $long long $ 范围内能保证正确。
Code
#include<bits/stdc++.h>
#define int __int128
#define N_4 10004
#define N_5 100005
#define N_6 1000006
#define Mod 1000000007
#define FOR(i,j,k) for(long long i=j;i<=k;++i)
#define ROF(i,j,k) for(long long i=j;i>=k;--i)
using namespace std;
inline int read(){
int x=0,f=1;
char ch;
ch=getchar();
while(ch<'0'||ch>'9'){
if(ch=='-') f=-f;
ch=getchar();
}
while(ch>='0'&&ch<='9'){
x=x*10+(ch-'0');
ch=getchar();
}
return x*f;
}
int T,n,tot;
int gg[6]={0,2,7,61,63,97};
bool isprime[500005];
int prime[500005],an[500005];
inline void Euler(int n){
isprime[1]=true;isprime[0]=true;
for(register int i=2;i<=n;i++){
if(!isprime[i])
prime[++tot]=i;
an[i]=tot;
for(register int j=1;j<=tot&&prime[j]*i<=n;j++){
isprime[i*prime[j]]=true;
if(i%prime[j]==0)
break;
}
}
return;
}
inline int ksm(int a,int b,int mod){
int ans=1,p=a,g=b;
while(g){
if(g&1) ans=(ans*p)%mod;
p=(p*p)%mod;
g>>=1;
}
return ans;
}
int mul(int a,int b,int p){return (a*b-(int)((__float128)a/p*b)*p+p)%p;}
inline bool miller_rabin(int P){
if(P==1) return 0;
int t=P-1,k=0;
while(!(t&1)) ++k,t>>=1;
for(register int i=1;i<=5;++i){
if(P==gg[i]) return true;
int a=ksm(gg[i],t,P),nxt=a;
for(register int j=1;j<=k;++j){
nxt=mul(nxt,nxt,P);
if(nxt==1&&a!=1&&a!=P-1) return false;
a=nxt;
}
if(a!=1) return false;
}
return true;
}
signed main(){
T=read();
Euler(500000);
while(T--){
n=read();
if(n<500000){
if(isprime[n]) cout<<"NO"<<endl;
else cout<<"YES"<<endl;
}
else{
if(!miller_rabin(n)) cout<<"NO"<<endl;
else cout<<"YES"<<endl;
}
}
return 0;
}
文章来源: 博客园
- 还没有人评论，欢迎说说您的想法！
客服
