Miller-Rabin素数测试

蒟蒻的第一篇博客，自然要写一些玄学的东西

 

前置知识：                                                                                                 

快速幂，费马小定理，龟速乘/O(1)快速乘 

 

问题探究：                                                                                                 

判断素数的方法有哪些呢？下面给大家几个方法，看看你会哪一种

1、循环判断法：时间复杂度O(n)

#include<bits/stdc++.h>
using namespace std;
int main()
{
    long long n;
    cin>>n;
    for(long long i=2;i<n;i++)
    if(n%i==0) {cout<<"No"; return 0;}
    cout<<"Yes";
    return 0;
}

      这种暴力算法应该没人用吧……

2、循环判断法（优化）：时间复杂度O()

      循环只用判断到就可以了，因为如果 i|n，则(n/i)|n

#include<bits/stdc++.h>
using namespace std;
int main()
{
    long long n;
    cin>>n;
    for(long long i=2;i<=sqrt(n);i++)
    if(n%i==0) {cout<<"No"; return 0;}
    cout<<"Yes";
    return 0;
}

     这也是通常使用的判断素数的方法，适用于大多数题目

 

但如果是一个较大的数，即使是的算法也会超时，例如下面这题

（模板）输入T，表示有T组数据，每组数据有一个整数n，判断其是否是素数，
输出“Yes”或“No”。
范围规定  T≤1000，n≤1×1019

这题中的大小也已经达到了10⁹，再加上多组数据，更是严重超时，那么我们就需要寻找一种效率更高的方法判断一个数是否是素数。

 

正文引入：                                                                                                 

在此之前需要了解费马小定理，即

对于质数p，任意整数a，均有
       ap-1≡1(mod p)  或   ap≡a(mod p)

详细可以见这里，或者当公式记下来

 

这样便给我们提供了一个思路，如果对于给定的一个数p,存在一个数a，使a^p-1≡1(mod p)不成立，那么p一定不是素数，我们可以任取一些数作为a，依次判断是否满足a^p-1≡1(mod p)，若对于所取的所有情况下式子都成立，那么p为素数的可能性就越大。通常情况下，我们一般取10次不同的a对p进行测试，若全部通过，则判定p素数。

若p通过一次测试，则p不是素数的概率为25%，若p通过了N次测试，则p不是素数的概率是1/4^N，若N取10时，p不是素数的概率已经远小于99.99%，在不是特别多的素数判断下可以大大提升效率，其中计算a的p-1次方时可以用快速幂实现，时间复杂度为O(k·logn)。

 

深入研究：                                                                                                 

1、假如模数p很大时，在快速幂中的两个数相乘会超过long long 范围时，通常会采用龟速乘或O(1)快速乘，防止在相乘中爆精度。

当然为了防止某些奇异错误，在数据范围不大时尽量少用。

2、如何进一步提升效率？既然都拿分数换时间了，还不想着怎么更快更好吗。。

在实际应用中，可以先预先筛出500个小素数（第500个素数为3571）对其先进行测试，进一步提升效率。

 

现在再回到一开始的题，下面是代码：

 

#include<bits/stdc++.h>
#define ll long long
#define ld long double
#define count 10
using namespace std;
ll t,n;
ll pr[600],npr[4000];
ll multi(ll x,ll y,ll mod)
{
    ll tmp=(x*y-(ll)((ld)x/mod*y+1.0e-8)*mod);
    return tmp<0?tmp+mod:tmp;
}
ll ksm(ll a,ll b,ll mod)
{
    ll ans=1;
    a%=mod;
    while(b)
    {
        if(b&1) ans=multi(ans,a,mod);
        b>>=1; a=multi(a,a,mod);
    }
    return ans;
}
void pre()
{
    int num=1;
    npr[0]=npr[1]=1;
    for(int i=2;num<501;i++)
    {
        if(!npr[i]) pr[num++]=i;
        for(int j=1;j<num&&i*pr[j]<4000;j++)
        {
            npr[i*pr[j]]=1;
            if(i%pr[j]==0) break;
        }
    }
    return;
}
bool miller(ll x)
{
    for(int i=1;i<=500&&pr[i]<=x;i++)
    if(pr[i]==x) return 1;
    else if(x%pr[i]==0) return 0;
    for(int i=1;i<=count;i++)
    {
        ll a=rand()%(n-2)+2;
        if(ksm(a,n-1,n)!=1) return 0;
    }
    return 1;
}
int main()
{
    srand(time(NULL));
    pre();
    cin>>t;
    while(t--)
    {
        cin>>n;
        cout<<(miller(n)?"Yes":"No")<<endl;
    }
    return 0;
}

 

 

其实这个算法不难理解，也比较好打，但为了避免某些毒瘤出题人，学一学总是有好处的