【算法杂谈】Miller-Rabin素性测试算法

额，我们今天来讲一讲Miller-Rabin素性测试算法。

读者：怎么又是随机算法！！！(⊙o⊙)…

【好了，言归正传】

【费马小定理】

费马小定理只是个必要条件，符合费马小定理而非素数的数叫做Carmichael

Carmichael数是非常少的。

在1~100000000范围内的整数中，只有255个Carmichael数。

为此又有二次探测定理，以确保该数为素数。

这就构成了Miller-Rabin的基本原理 ╰(￣▽￣)╭

【二次探测定理】

二次探测定理 如果p是一个素数，0<x<p,则方程x^2≡1(mod p)的解为x=1,p-1

【Miller-Rabin】

让我们来模拟一下：

第一步：先计算出m、j，使得n-1=m*2^j，其中m是正奇数，j是非负整数

第二步：随机取一个b，2<=b

第三步：计算v=b^m mod n

第四步：如果v==1，通过测试，返回

第五步：令i=1

第六步：如果v=n-1，通过测试，返回

第七步：如果i==j，非素数，结束

第八步：v=v^2 mod n，i=i+1

第九步：循环到S（S根据数据进行改变）

【算法分析】

既然是随机算法，那么我们就来计算一下这种算法的正确率。

Miller-Rabin的正确率为75%。

怎么这么低(⊙o⊙)…，你可以多调用几次就好了！这样可以使正确概率提高为1-(1/4)^S

【代码】

#include<iostream>
#include<cstdio>
#include<cstring>
#include<windows.h>
#include<time.h>
using namespace std;
#define ll long long
const int S=5; //可以根据数据进行更改！！！ 
ll mod_mul(ll a, ll b, ll n) 
{
    ll res=0;
    while(b) 
    {
        if(b&1) res=(res + a)%n;
        a=(a+a)%n;
        b>>=1;
    }
    return res;
}

ll mod_exp(ll a, ll b, ll n) 
{
    ll res=1;
    while(b) 
    {
        if(b&1) res=mod_mul(res,a,n);
        a=mod_mul(a,a,n);
        b>>=1;
    }
    return res;
}
bool miller_rabin(ll n) 
{
    if(n==2 || n==3 || n==5 || n==7 || n==11) return true;
    if(n == 1 || !(n%2) || !(n%3) || !(n%5) || !(n%7) || !(n%11)) return false;
    ll x,pre,u;
    int i,j,k=0;
    u=n-1;    
    //要求x^u % n
    while(!(u&1))
    {   //如果u为偶数则u右移，用k记录移位数
        k++;
        u >>= 1;
    }
    srand((ll)time(0));
    for(i=0;i<S;++i) 
    {   //进行S次测试
        x=rand()%(n-2) + 2;    
        //在[2, n)中取随机数
        if((x%n)==0) continue;
        x=mod_exp(x,u,n);    
        //先计算(x^u) % n，
        pre=x;
        for(j=0;j<k;++j) 
        {   //把移位减掉的量补上，并在这地方加上二次探测
            x=mod_mul(x,x,n);
            if(x==1 && pre!=1 && pre!=n-1) return false;    
            //二次探测定理，这里如果x = 1则pre 必须等于 1，或则 n-1否则可以判断不是素数
            pre=x;
        }
        if(x!=1) return false;    
        //费马小定理
    }
    return true;
}
int main()
{
    long long n;
    scanf("%lld",&n);
    cout<<miller_rabin(n);
    system("pause");
    return 0;
}

P.S. 感谢吉林大学的模板，当然这已经改过了，进行了优化。

【你以为这就完了？】

下面我们来研究一下Miller-Rabin的更高效算法

【代码】

#include<windows.h>
#include<iostream>
#include<time.h>
#include<cstdio>
#include<cmath>
using namespace std;
bool witness(long long a,long long n)
{
    long long t,d,x;
    d=1;
    int i=ceil(log(n-1.0)/log(2.0)) - 1;
    for(;i>=0;i--)
    {
        x=d;  d=(d*d)%n;
        if(d==1 && x!=1 && x!=n-1) return true;
        if( ((n-1) & (1<<i)) > 0)
            d=(d*a)%n;
    }
    return d==1? false : true;
}
bool miller_rabin(long long n)
{
    int s[]={2,7,61};
    if(n==2) return true;
    if(n==1 || ((n&1)==0)) return false;
    for(int i=0;i<3;i++)//注意这里的i<3不可以改，因为每次S都只为2,7,61！ 
        if(witness(s[i], n)) return false;
    return true;
}
int main()
{
    long long n;
    scanf("%lld",&n);
    cout<<miller_rabin(n);
    //system("pause");
    return 0;
}

我们发现，它只调用witness函数3次，每次都a为2,7,61。时间很快啊。

但正确率……

:-D   //迷之微笑