[转载]Hastin's 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\)为质数,则对于\(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;
}