Miller_Rabin素数测试

Miller_Rabin素数测试

算法用处

Miller_Rabin为随机算法，有一定的误判率，算法中会有常数S，S越大算法耗时越大，误判率越小。S这里选20

误判率 = 1 / 2^S，大概等于9*10^-7

由于之前的素数判断通常经过埃氏筛法，线性筛选出素数再在O(1)的时间内判断，但是素数很大的时候打表非常耗时，非常耗空间，所以用Miller_Rabin素数测试，虽然有很小的误判率，但是效率很高，O(logN)就能出结果，有很高的实用性。但是对于那些极其注重安全的地方不能使用。应为例如核电站控制，核武器控制。一单出错就容易发生严重后果。这里扯远了。。。。

算法原理

1.费马小定理

对于素数p来说任意一个整数a，总是满足a^p≡a(mod p)，特别的当gcd(a, p) = 1(最大公约数)的时候又有a^(p - 1)≡1(mod p)。所以a^k≡a^(k % (p - 1)) (mod p)

2.二次探测定理

当存在X，取值范围为0 < x < p时，x^2≡1(mod p)的解为x = 1或x = p - 1。因为除了1和p-1在mod p的情况下的逆元等于自己，其他的逆元都不是自身。   (有兴趣深入的话请看威尔逊定理，这里只是教怎么使用和编码)

算法过程

首先我们根据原理1，找出这个数a，那么求的就是a^(p - 1)，这里将p - 1里面的所有的2提取出来变成形如p - 1 = d * 2^t，原式就成为了(a ^ d)^(2^t)那么就成了将a^d平方t次，那么每次平方的时候就可以使用二次探测加强判定准确性，最后结果看是不是和1对于p同余，如果不满足肯定是合数。

这里就涉及到快速幂以及为了防止溢出的类似快速幂的乘法，可以看multi_mod()和pow_mod()。

如果Miller_Rabin测试出来某个数是合数的时候它一定是合数，如果是素数的话有1 / 2^S的概率为合数

HDU 2138关于素数的可以试试自己写的Miller_rabin测试

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
#include<ctime>
#include<cstdlib>
#define ll __int64

const int S = 20;//素数测试误判率=1 / 2^S   S越大误差率越小

using namespace std;


ll multi_mod(ll a, ll b, ll p){
    ll ret = 0, q = a;
    while(b){
        if(b & 1){
            ret = (ret + q) % p;
        }
        q = (q + q) % p;
        b >>= 1;
    }
    return ret % p;
}

ll pow_mod(ll a, ll b, ll p){
    ll ret = 1, q = a;
    while(b){
        if(b & 1){
            ret = ret * q % p;
        }
        q = q * q % p;
        b >>= 1;
    }
    return ret % p;
}

bool Miller_Rabin(ll p){
    int i, j, k;
    ll u, t, x, y;
    if(p == 2)  return true;
    if(p % 2 == 0 || p == 1) return false;
    u = p - 1, t = 0;
    while(u % 2 == 0){
        t++;
        u >>= 1;
    }
    for(int i = 0; i < S; i++){
        x = rand() % (p - 1) + 1;
        x = pow_mod(x, u, p);
        for(j = 0; j < t; j++){
            y = x;
            x = multi_mod(x, x, p);
            if(x == 1 && y != 1 && y != p - 1){
                return false;
            }
        }
        if(x != 1)  return false;
    }
    return true;
}



int main(){
    int i, j, k;
    ll a, b, p;
    int n;
    srand(time(NULL));
    while(cin >> n){
        int ans = 0;
        for(i = 0; i < n; i++){
            cin >> p;
            if(Miller_Rabin(p)) ans++;
        }
        cout << ans << endl;
    }
    return 0;
}