Miller-Rabin大素数判定

[原文](http://www.cnblogs.com/vongang/archive/2012/03/15/2398626.html)
#include <algorithm>
#include <iostream>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
using namespace std;
const int MAX = 1010;
const int INF = 0x3f3f3f3f;

#define S 10
typedef long long ll;
ll mod_mul(ll a, ll b, ll n);
ll mod_exp(ll a, ll b, ll n);
bool Miller_Rabin(ll n);
bool witness(ll a, ll b, ll n,int t);

int main()
{
    int N, ans, num;
    while(~scanf("%d", &N))
    {
        ans=0;
        for(int i=1; i<=N; i++)
        {
            scanf("%d", &num);
            if(Miller_Rabin(num))ans++;
        }
        printf("%d\n", ans);
    }
    return 0;
}

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, u=n-1;
    int k=0;

    while(!(u&1))
    {
        k++;
        u>>=1;
    }

    for(int i=1; i<=S; i++)
    {
        x=rand()%(n-2)+2;  
        if(x%n==0)continue;

        if(!witness(x,u,n,k))return false;
    }
    return true;
}

ll mod_mul(ll a, ll b, ll n)  \\计算a*b%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)  \\计算a^b%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;
}
\\令n-1=2^t*u, u是奇数，即n-1的二进制表示数是奇数u的二进制表示数后

面加上t个零，a^(n-1)=(a^(u*2^t))(mod n),所以可以先求a^u(mod n),
然后把结果平方后取模t次，就是a^(n-1)(mod n)。其中a^u(mod n),结果

平方过程中，会有公式可以减少判断次数。x*x(mod n)=1,解为x=1||x=n-1.


bool witness(ll a, ll b, ll n,int t)
{
    ll pre, cur;

    cur=mod_exp(a,b,n);
    pre=cur;
    for(int i=1; i<=t; i++)
    {
        cur=mod_mul(cur,cur,n);
        if(cur==1&&pre!=1&&pre!=n-1)return false;
        pre=cur;
    }

    if(cur!=1)return false;
    else return true;
}