假装学习Miller_Rabin

Miller_Rabin素性测试

---

不是，说真的，数学真难学。

别人都是把参考资料放在结尾，我就放开头。

参考1

参考2

感谢两篇对本蒟蒻的帮助。

1. 什么是素性测试？

正常情况下，判断一个数是否为素数可以使用一个 $ O( \sqrt{n}) $ 的算法。
但是作为一个 $ OIer $ ，怎么可能跟正常扯上关系。于是ta来了。

这个数据量显然不是$ O( \sqrt{n}) $ 能拿下的，于是就得看我们的 $ Miller_Rabin $ 了。

算法逻辑

温故而知新：

1.费马小定理：如果 $ p $ 为素数，则满足 $ a^{p-1} \equiv 1 ( mod \ p ) $
2.对于一个素数 $ p $，存在一个 $ x，0<x<p $ ，使得方程 $ x^{2} \equiv 1 ( mod \ p) $，解有且只有 $ x=1 $ 或 $ x=p-1 $ 。

---

对于第二点的证明

由 $ x^{2} \equiv 1 ( mod \ p ) $ ，可以推得 $ (x-1) \cdot (x+1) mod \ p = 0 $。由于此处的 $ p $ 是一个素数，又因为 $ x<p $ ，所以 $ (x-1) $ 和 $ (x+1) $ 就都不是 $ p $ 的因子，那么 $ (x-1) \cdot (x+1) %p $ 就不会为 $ 0 $ ，因此要让 $ x^{2} \equiv 1 ( mod \ p ) $ 成立，只有让 $ (x-1) \cdot (x+1) == 0 $ ，于是就可以得到两个解 $ x=1 $ 和 $ x=p-1 $ 。

---

二次探测定理

而这也就是二次探测定理，但是难不成让我一个个 $ x $ 去判断吗？那不 $ T $ 飞到下辈子？从反面考虑，咱去把结果为 $ 1 $ 的 $ x $ 拎出来，看下他们是不是等于 $ p-1 $ 或 $ 1 $ 。

这个时候就要关注一下我们的费马小定理了，一个数要是个素数就必须得通过费马素性测试，即满足 $ a^{p-1} \equiv 1 ( mod \ p ) $ 。但反过来则不一定成立，而反例就是 $ Carmichael $ 数。

注意到 $ p-1 $ 是一个偶数，于是他就可以写成 $ d \cdot 2^{t} $ 的形式，那么原式就可以表述为 $ a^{d \cdot 2^{t}} \equiv 1 ( mod \ p ) $，将 $ a^{d} $ 看成 $ k $ ，那么就得到了一个简洁的 $ k^{2{t}} \equiv 1 ( mod \ p ) $ 。同时 $ k^{2{t}} $ 也可以看成 $ (k \times k…… \times k)^{2} $ 。

再来看一下我们的目标，此时这就发现所有要找的 $ x $ 存在于 $ k^{2{t}} $ 的变化过程中，所以我们只要在过程中判断，只要其等于 $ 1 $ 或 $ p-1 $即可，但是要想严格确定一个数为素数，就得遍历所有的 $ k $ ，那么必然叒要 $ T $ 。然后我们可以牺牲一点点正确率，从而在时间上有巨大收获，也就是仅尝试少数个 $ k $ ，这样我们用可以忽略不计的错误概率换取一个跑得飞快的算法，赢麻了。

补充：在一定范围内，概率可以达到 $ 100% $ ,

1. $ if(n<1373653) a=2,3 $ ;
2. $ if(n<9080191) a=31,73 $ ;
3. $ if(n<4759123141) a=2,7,61 $ ;
4. $ if(n<2152302898747) a=2,3,4,5,11 $ 。

总结

核心部分：

如果n为一个素数，则 $ \forall a<n $ ，设 $ n-1=d \cdot 2^{t} $ ， $ d $ 为奇数，则以下两命题至少一个成立：

1. $ a^{d} \equiv 1 ( mod \ n ) $
2. $ \exists 0 \leq i < t ， (a^{d}){2^{i}} \equiv -1 ( mod \ n ) $

注意这里是存在命题！！！

复杂度：$ O(k \cdot log (n) ) $ ， $ k $ 为所取测试数的个数。

---

完整代码

#ifdef ONLINE_JUDGE 
#else 
#define Qiu_Cheng 
#endif 
#include <bits/stdc++.h> 
#define i8  __int128
#define int long long 
#define fuck inline
#define lb long double 
using namespace std; 
// typedef longlong ll; 
const int N=1e6+5,M=1e6+520,mod=1e9+7;
const int inf=INT_MAX,INF=1e9+7; 
// const int mod1=469762049,mod2=998244353,mod3=1004535809;
// const int G=3,Gi=332748118; 
// const int M=mod1*mod2;
fuck int read()
{
    int x=0,f=1;
    char c=getchar();
    while(c<'0'||c>'9'){if(c=='-'){f=-1;}c=getchar();}
    while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+(c-'0');c=getchar();}
    return x*f;
}
fuck void write(int x)
{
    if(x<0){putchar('-');x=-x;}
    if(x>9) write(x/10);
    putchar(x%10+'0');
}
//线性筛
int prime[N],is_prime[N]={1,1},tot=0;
fuck void xxs()
{
    for(int i=1;i<=N;i++)is_prime[i]=0;
    for(int i=2;i<=N;i++)
    {
        if(!is_prime[i])prime[++tot]=i;
        for(int j=1;j<=tot&&i*prime[j]<=N;j++)
        {
            is_prime[i*prime[j]]=1;
            if(i%prime[j]==0)break;
        }
    }
}
// 龟速乘
fuck int gsc(int a,int b,int Mod)
{
    a%=Mod,b%=Mod;
    int res=0;
    while(b)
    {
        if(b&1)res=(res+a+Mod)%Mod;
        b>>=1;a=(a+a+Mod)%Mod;
    }
    return (res+Mod)%Mod;
}
//快速幂
fuck int ksm(int a,int b,int Mod)
{
    int res=1;
    while(b)
    {
        if(b&1)res=gsc(res,a,Mod);
        b>>=1;a=gsc(a,a,Mod);
    }
    return (res+Mod)%Mod;
}
fuck int Miller_Rabin(int n,int a)
{
    int m=n-1,t=0;
    while((m&1)==0)m>>=1,t++;
    int k=ksm(a,m,n);
    if(k==1)return 1;//特判一下
    for(int i=1;i<=t;i++,k=gsc(k,k,n))
        if(k==n-1)return 1;//如果存在n-1的解，说明其为素数
    return 0;
}
fuck int isprime(int n)
{
    if(n<=1)return 0;
    for(int i=1;i<=20;i++)
    {
        if(n==prime[i])return 1;
        if(n%prime[i]==0||!Miller_Rabin(n,prime[i]))return 0;
    }
    return 1;
}
fuck void solve()
{
    int n=read();
    if(isprime(n))cout<<"YES"<<"\n";
    else cout<<"NO"<<"\n";
}
signed main() 
{ 
#ifdef Qiu_Cheng 
    freopen("1.in","r",stdin); 
    freopen("1.out","w",stdout); 
#endif 
    // ios::sync_with_stdio(false); 
    // cin.tie(0); cout.tie(0); 
    // int fuckccf=read();
    xxs();
    int QwQ=read();
    while(QwQ--)solve(); 
    // solve(); 
    return 0; 
}
//  6666   66666  666666 
// 6    6  6   6      6 
// 6    6  6666      6 
// 6    6  6  6    6 
//  6666   6   6  6666666

---

完结收工！！！！！

个人主页

看完点赞，养成习惯

\(\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\Downarrow\)