[转载] miller rabin

本文转载自https://www.cnblogs.com/zsq259/p/11602175.html

Miller-Rabin

事先声明,因为菜鸡Hastin知识水平有限就是菜,因此语言可能不是特别规范,仅供理解.

step 0

问一个数(p)是否为质数,(p<=10^{18}).

一个简单暴力的办法是(O( sqrt{n}))枚举约数.

然而显然会炸.

于是我们就有了Miller-Rabin.

讲了好多废话...

step 1

首先了解一下费马小定理:

若(p)为质数,则对于(ain[1,p-1])有(a^{p}equiv a(Mod) (p))

那么就有(a^{p-1} equiv 1( Mod) (p))

下面我们用数学归纳法来证明一下费马小定理:

显然(a=1)时结论成立,

若(a=n)时结论成立,

当(a=n+1)时,有

((n+1)^p)(=sum_{i=0}^{p}C_{p}^{i}n^{p-i})(二项式定理)

那么除了(n^{p})和(1)这两项外,

其它的都有一个系数(C_{p}^{i},iin[1,p-1]),所以都能被(p)整除.

而(n^pequiv n(Mod) (p)),

所以((n+1)^{p}equiv n+1(Mod) (p)),结论成立.

所以回到正题,如果对于一个数(p),存在(a in [1,p-1]),(a^{p-1} otequiv1(Mod) (p)),则(p)一定不是质数.

然而仍然有一些数能够逃掉这种检测,

于是就有了

step 2

二次探测!

对于一个质数(p),若(a in[1,p-1]),且(a^2equiv1(Mod) (p)),则(a=1)或(a=p-1).

证明:

若(a^2equiv1(Mod) (p)),则(a^2-1equiv 0(Mod) (p))

即(p|(a+1)(a-1)).

因为(p)为质数,且(ain [1,p-1]),所以只有当(a=1)或(a=p-1)时上式才成立.

所以反过来想当(a)等于其它数时,(p)就不是质数了.

step 3

下面再来讲一下具体的实现.

首先大的思路是费马小定理,在中间用二次探测判定.

对于要判定的数(p),

设(u=p-1),则根据费马小定理有(a^uequiv 1(Mod) (p)),(ain [1,p-1]).

然后把(u)写成(d*2^n)的形式,也就是((((d^2)^2)^2)^{2...})(n个2)

那么从(d^2)开始用二次探测判定,

最后再用费马小定理判定就行了.

而时间复杂度是(O(log_{ iny 2}) (p)).

并且一次的失误率是(frac{1}{4}),

那么(T)次测试的失误率就是(4^{-T}),时间复杂度为(O(Tlog_{ iny 2}p)),能够接受.

并且如果(n<2^{64}),只用选取(a=2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)测试即可.

code:(似乎要开int128)

#include <iostream>
#include <cstdio>
#include <cstring>
#define int __int128
#define fre(x) freopen(x".in","r",stdin),freopen(x".out","w",stdout)
using namespace std;

inline int read(){
	int sum=0,f=1;char ch=getchar();
	while(ch>'9'||ch<'0'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){sum=sum*10+ch-'0';ch=getchar();}
	return f*sum;
}

int p[101]={2,3,5,7,11,13,17,19,23,29,31,37};
int T,n,tot=12;

inline int fpow(int a,int b,int p){
	int ret=1;
	for(;b;a=a*a%p,b>>=1) if(b&1) ret=ret*a%p;
	return ret;
}

inline bool Miller_Rabin(int n){
	if(n<=2){
		if(n==2) return 1;
		return 0;		
	}
	int u=n-1;
	while(!(u%2)) u>>=1;
	int d=u;
	for(int i=1;i<tot&&p[i]<n;i++){
		int x=p[i];u=d;x=fpow(x,u,n);		
		while(u<n){
			int y=fpow(x,2,n);
			if(y==1&&x!=1&&x!=n-1) return 0;
			x=y;u<<=1;
		}
		if(x!=1) return 0;
	}
	return 1;
}

signed main(){
	T=read();
	while(T--){
		n=read();
		if(Miller_Rabin(n)) puts("Yes");
		else puts("No");
	}
	return 0;
}