Miller-Rabin 素数测试

相关定理——费马小定理：假设P是素数，且(a,p)=1，那么
由此我们知道这样一个事实： p是素数，(a,p)=1 -> ； -> p不是素数
定义：a是正整数，p是合数，且 ， 那么称p是以a为基的伪素数。
Miller-Rabin算法原理： 取多个a（底）进行试验，次数越多，p是素数的概率越大。

相关例题：
POJ 3641 Pseudoprime numbers

http://poj.org/problem?id=3641

求解p是否是以a为基的伪素数

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
typedef long long LL;
bool ispri(LL x){
    bool yes=1;
    for(int i=2;i<=(LL)sqrt(double(x));i++){
        if(x%i==0){
            yes=0;
            break;
        }
    }
    return yes;
}
LL quick_mod(LL p,LL a){
    LL ans=1,pow=p;
    while(pow){
        if(pow&1) ans=ans*a%p;
        a=a*a%p;
        pow>>=1;
    }
    return ans;
}
int main()
{
    LL p,a;
    while(cin>>p>>a&&(p+a)){
        if(ispri(p)==false&&quick_mod(p,a)==a) puts("yes");
        else puts("no");
    }
    return 0;
}

然而，事情没有这样简单。。有一类神奇的数字，它叫Carmichael数。
它具有这样的性质：对于Carmichael数p，它是合数，如果对于所有正整数a，(a，p)=1，都有同余式 成立。
所以，在进行Miller-Rabin素数测试时，需要进行Carmichael数的排除操作。
二次探测定理：如果P是一个素数，而且0<x<P，则方程 的解是x=1或者x=p-1.
所谓的Carmichael数排除操作也即是检测P是否满足这样的性质。
即将产生的比P小的随机数字，进行幂取模操作 带入x,随后不断验证。
详见代码：

hdu 2138  How many prime numbers

   
   

    
    http://acm.hdu.edu.cn/showproblem.php?pid=2138
   
   
   
   

    
    求解给定的一堆数字中有多少素数？
   
   

#include <iostream>
#include <stdlib.h>
#include <cstdio>
#include <time.h>
using namespace std;
typedef long long LL;
LL multi(LL a,LL b,LL p){
    LL ans=0;
    a=a%p;
    while(b){
        if(b&1) ans=(ans+a)%p;
        a=(a+a)%p;
        b>>=1;
    }
    return ans;
}
LL quick_mod(LL a,LL m,LL p){
    LL ans=1;
    a=a%p;
    while(m){
        if(m&1) ans=multi(ans,a,p);
        a=multi(a,a,p);
        m>>=1;
    }
    return ans;
}
bool Miller_Rabin(LL p){
    if(p==2) return 1;
    if(p<2 || (p&1)==0) return 0;
    LL m=p-1;
    int sum=0;
    while((m&1)==0){
        m>>=1;
        sum++;
    }
    //srand(time(NULL)); 用不着
    for(int i=0;i<10;i++){
        LL a=rand()%(p-1)+1;  // 产生1-P-1的随机数
        LL x=quick_mod(a,m,p);
        LL g=0;
        for(int j=0;j<sum;j++){
            g=multi(x,x,p);
            if(g==1&&x!=1&&x!=p-1) return 0;
            x=g;
        }
        if(g!=1) return 0;
    }
    return 1;
}
int main()
{
    int t;
    LL p;
    while(cin>>t){
        int ans=0;
        for(int i=0;i<t;i++){
            scanf("%lld",&p);
            if(Miller_Rabin(p)) ans++;
        }
        printf("%d\n",ans);
    }
    return 0;
}