Miller Rabin 大素数测试

　　PS:本人第一次写随笔，写的不好请见谅。

　　接触MillerRabin算法大概是一年前，看到这个算法首先得为它的神奇之处大为赞叹，竟然可以通过几次随机数据的猜测就能判断出这数是否是素数，虽然说是有误差率，但是相对于他比其他素数判断的高效，真的是可以说是完美级。那时候忙于找工作，所以也没有细究，现在空下来终于对这个算法有了一定的理解。

　　先说两个定理:

　　　　(1) 当x<p时，满足x^(p-1) % p = 1，说明x与p互质；

　　　　(2) 当x<p时，满足x^2 % p = 1;　　x的解为 x = 1 或者 x = p -1；

　　根据第一个定理，我们可以在[2,p-1)中间随机选择一个数，然后根据定理1判断是否是素数（测试准确率跟测试次数有关，多次检测）。定理2可以做二次探索。

　　

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<ctime>
using namespace std;
#define S 10

long mod_mul(long a,long b,long n) {
    long res = 0;
    while(b) {
        if(b&1)
            res = (res + a) % n;
        a = (a + a) % n;
        b >>= 1;
    }
    return res;
}

long mod_exp(long a, long b, long n) {
    long 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(long 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;

    long x, pre, u;
    int i, j, k = 0;
    u = n - 1;    //要求x^u % n

    while(!(u&1)) {    //如果u为偶数则u右移，用k记录移位数
        k++; u >>= 1;
    }

    srand(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()
{
    for(int i=2;i<100;++i)
        if(miller_rabin(i))
            cout<<i<<" ";
    return 0;
}

　　mod_mul与mod_exp是关于求模。不要看这两个函数很高大上，其实换种形式会更加理解。

long mod_mul(long a,long b,long n) {

    //求 (a*b)%n
    //在不考虑数据上移时，先求a*b的值
    long res = 0;  // res表示a*b的积
   //举例 5（十进制） * 1101(2进制) = 5 * 2^0 + 5*0*2^1 + 5*2^2+5*2^3;
    while(b) {
        if(b&1)
            res = res + a;
        a = 2*a;
        b >>= 1;
    }
    return res%n;     //最后求模
}

long mod_exp(long a, long b, long n) {
    long res = 1; 
  // res = a^b,把b分解成二进制
    while(b) {
        if(b&1)
            res = res * a;
        a = a*a;
        b >>= 1;
    }
    return res%n;
}    

现在讲米勒算法：

　　1、n为待判定数字，先取一个数u = n - 1;

     2、通过循环除二，分解成 u（原） = u(新) * 2^k;

     3、随机出x，x为[2,n),判断 x ^(u*2^k)%n == 1, 式子分解  [(x^u)%n ]^(2^k) %n == 1; 先计算 x= x ^ u % n;然后计算 x^(2^k) % n == 1;

     4、如何计算x^(2^k) = (x^2)^(2^(k-1));  并且可知，当x^2 % n == 1时，无需计算2^(k-1);  当x^2 % n == 1时， x ! = 1 或者 x != n-1则可判断该数为合数

     5、循环测试，如果最后都没有return false，则可认为该数时质数。

例题有hdu2138：

　　注意点：虽然题目中数字都在32位整型之内，但是再做加法的时候，大小会超。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<ctime>
using namespace std;
//a容易越界,res容易越界。因为有加号。
int mod_mul(__int64 a,int b,int n)
{
    __int64 res = 0;
    while(b)
    {
        if(b&1) res = (res+a)%n;
        a = (2*a)%n;
        b >>= 1;
    }
    return (int)res;
}

int mod_exp(int a,int b,int n)
{
    int res = 1;
    while(b)
    {
        if(b&1) res = mod_mul(res,a,n);
        a = mod_mul(a,a,n);
        b >>= 1;
    }
    return res;
}

bool rabin_miller(int n)
{
    if(n==2 || n==3 || n==5 || n==7 || n==11 || n==13) return true;
    if(n<=1 || n%2==0 || n%3==0 || n%5==0 || n%7==0 || n%11==0 || n%13==0) return false;
    int u,k,x,pre;
    u = n-1;
    k = 0;
    while(u&1 ==0)
    {
        ++k;
        u >>= 1;
    }
    srand(time(0));
    int s = 3;
    while(s--)
    {
        x = rand()%(n-2)+2;
        pre = x = mod_exp(x,u,n);
        for(int i=0;i<k;++i)
        {
            x = mod_mul(x,x,n);
            if(x ==1 && pre != 1 && pre != n-1) return false;
            pre = x;
        }
        if(x != 1) return false;
    }
    return true;
}

int main()
{
    int n;
    while(~scanf("%d",&n))
    {
       int ans = 0;
       for(int i=0;i<n;++i)
       {
           int a;
           scanf("%d",&a);
           if(rabin_miller(a))
               ++ans;
       }
       printf("%d
",ans);
    }

    return 0;
}