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