Miller-Rabin素数测试

素数的测试:

 

费尔马小定理:如果p是一个素数,且0<a<p,则a^(p-1)%p=1.
            利用费尔马小定理,对于给定的整数n,可以设计素数判定算法,通过 计算d=a^(n-1)%n来判断n的素性,当d!=1时,n肯定不是素数,当d=1时,n   很可能是素数.

 

二次探测定理:如果p是一个素数,且0<x<p,则方程x^2%p=1的解为:x=1或    x=p-1.
            利用二次探测定理,可以再利用费尔马小定理计算a^(n-1)%n的过程 中增加对整数n的二次探测,一旦发现违背二次探测条件,即得出n不是素数的结论.
   
    如果n是素数,则(n-1)必是偶数,因此可令(n-1)=m*(2^q),其中m是正奇数( 若n是偶数,则上面的m*(2^q)一定可以分解成一个正奇数乘以2的k次方的形式 ),q是非负整数,考察下面的测试:
    序列:
         a^m%n; a^(2m)%n; a^(4m)%n; …… ;a^(m*2^q)%n
    把上述测试序列叫做Miller测试,关于Miller测试,有下面的定理:

定理:若n是素数,a是小于n的正整数,则n对以a为基的Miller测试,结果为真.
Miller测试进行k次,将合数当成素数处理的错误概率最多不会超过4^(-k).

 

Miller_Rabin测试T次时，它产生一个假的素数所花费的时间不超过1/4^T
一些Carmichael数:
3
561 1105 1729
6
294409 56052361 118901521 172947529 216821881 228842209

 

15MS:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

 

typedef long long LL;

LL exp_mod(LL a,LL b,LL m)
{
    bool bit[64];
    int i=0;
    while(b)
    {
        bit[i++] = b&1 ? 1 : 0;
        b >>= 1;
    }
    LL s=1;
    for(i--; i>=0; i--)
    {
        s = s * s % m;
        if(bit[i]) s = s * a % m;
    }
    return s;
}

bool Witness(LL a,LL n)
{
    LL m = n-1;
    int j=0;
    while( !(m&1) )
    {
        j++;
        m >>= 1;
    }
    LL x = exp_mod(a,m,n);
    if(x == 1 || x == n-1) return false; // n may be prime
    while(j--)
    {
    //    x = exp_mod(x,2,n);
        x = x * x % n;
        if(x == n-1) return false;
    }
    return true; //n must be composite
}

/*
bool Witness(LL a,LL n)
{
    LL u=n-1;
    int t=0;
    while(!(u&1))
    {
        t++;
        u>>=1;
    }
    LL y,x=mod_exp(a,u,n);
    while(t--)
    {
        y = x * x % n;
        if(y==1 && x!=1 && x!=n-1) return true;
        x = y;
    }
    if(y!=1) return true; // composite
    return false; //may be prime
}
*/

bool Miller_Rabin(LL n,int T=1)
{
    if(n<2) return false;
    if(n==2) return true;
    if( !(n&1) ) return false;
    while(T--)
    {
        LL a = rand() * (n-2) / RAND_MAX + 1; //[1,n+1]
        if( Witness(a,n) ) return false;
    }
    return true;
}

int main()
{
#ifndef ONLINE_JUDGE
    freopen("tdata.txt","r",stdin);
#endif
    srand(unsigned(time(0)));
    int n,cnt;
    LL p;
    while(scanf("%d",&n)!=EOF)
    {
        cnt=0;
        while(n--)
        {
            scanf("%I64d",&p);
            if(Miller_Rabin(p)) cnt++;
        }
        printf("%d\n",cnt);
    }
    return 0;
}