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 = ;
while(b) {
if(b&)
res = (res + a) % n;
a = (a + a) % n;
b >>= ;
}
return res;
} long mod_exp(long a, long b, long n) {
long res = ;
while(b) {
if(b&)
res = mod_mul(res, a, n);
a = mod_mul(a, a, n);
b >>= ;
}
return res;
} bool miller_rabin(long n) { //过滤器
if(n == || n == || n == || n == || n == ) return true;
if(n == || !(n%) || !(n%) || !(n%) || !(n%) || !(n%)) return false; long x, pre, u;
int i, j, k = ;
u = n - ; //要求x^u % n while(!(u&)) { //如果u为偶数则u右移，用k记录移位数
k++; u >>= ;
} srand(time());
for(i = ; i < S; ++i) { //进行S次测试
x = rand()%(n-) + ; //在[2, n)中取随机数
if((x%n) == ) continue; x = mod_exp(x, u, n); //先计算(x^u) % n，
pre = x;
for(j = ; j < k; ++j) { //把移位减掉的量补上，并在这地方加上二次探测
x = mod_mul(x, x, n);
if(x == && pre != && pre != n-) return false; //二次探测定理，这里如果x = 1则pre 必须等于 1，或则 n-1否则可以判断不是素数
pre = x;
}
if(x != ) return false; //费马小定理
}
return true;
} int main()
{
for(int i=;i<;++i)
if(miller_rabin(i))
cout<<i<<" ";
return ;
}
mod_mul与mod_exp是关于求模。不要看这两个函数很高大上，其实换种形式会更加理解。
long mod_mul(long a,long b,long n) {
//求 (a*b)%n
//在不考虑数据上移时，先求a*b的值
long res = ; // res表示a*b的积
//举例 5（十进制） * 1101(2进制) = 5 * 2^0 + 5*0*2^1 + 5*2^2+5*2^3;
while(b) {
if(b&)
res = res + a;
a = *a;
b >>= ;
}
return res%n; //最后求模
}
long mod_exp(long a, long b, long n) {
long res = ;
// res = a^b,把b分解成二进制
while(b) {
if(b&)
res = res * a;
a = a*a;
b >>= ;
}
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 = ;
while(b)
{
if(b&) res = (res+a)%n;
a = (*a)%n;
b >>= ;
}
return (int)res;
} int mod_exp(int a,int b,int n)
{
int res = ;
while(b)
{
if(b&) res = mod_mul(res,a,n);
a = mod_mul(a,a,n);
b >>= ;
}
return res;
} bool rabin_miller(int n)
{
if(n== || n== || n== || n== || n== || n==) return true;
if(n<= || n%== || n%== || n%== || n%== || n%== || n%==) return false;
int u,k,x,pre;
u = n-;
k = ;
while(u& ==)
{
++k;
u >>= ;
}
srand(time());
int s = ;
while(s--)
{
x = rand()%(n-)+;
pre = x = mod_exp(x,u,n);
for(int i=;i<k;++i)
{
x = mod_mul(x,x,n);
if(x == && pre != && pre != n-) return false;
pre = x;
}
if(x != ) return false;
}
return true;
} int main()
{
int n;
while(~scanf("%d",&n))
{
int ans = ;
for(int i=;i<n;++i)
{
int a;
scanf("%d",&a);
if(rabin_miller(a))
++ans;
}
printf("%d\n",ans);
} return ;
}