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 10typedef 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;}