[转载] miller rabin

Miller-Rabin

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

step 0

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

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

然而显然会炸.

于是我们就有了Miller-Rabin.

step 1

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

若$p$为质数,则对于$a\in[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},i\in[1,p-1]$,所以都能被$p$整除.

而$n^p\equiv n(Mod$ $p)$,

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

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

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

于是就有了

step 2

二次探测!

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

证明:

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

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

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

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

step 3

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

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

对于要判定的数$p$,

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

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

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

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

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

并且一次的失误率是$\frac{1}{4}$,

那么$T$次测试的失误率就是$4^{-T}$,时间复杂度为$O(Tlog_{\tiny 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;
}