Miller-Rabin算法

转自：http://blog.csdn.net/Sunshine_cfbsl/article/details/52425798

判定一个整数n(n>1)是否为素数

算法1:

直接根据素数的定义枚举i从2到(n-1)，如果n%i==0，n为合数。
时间复杂度：O($n$)

bool is_prime(int n) {
    int i;
    for(i = 2; i < n; i++)
        if(n % i == 0) return false;
    return true;
}

---

算法2:

发现若存在$i<n$使得n%i==0，则必有n%(n/i)==0。
所以只需枚举i从2到$\sqrt n$即可。
时间复杂度：O($\sqrt n$)

bool is_prime(int n) {
    int i;
    for(i = 2; i * i <= n; i++)
        if(n % i == 0) return false;
    return true;
}

---

Miller-Rabin算法:

这是一种随机性素数判定算法，也就是说，答案可能出错，但是可能性极小。
先是讲两个定理：

• 费马小定理：

对于一个质数$p$，取任意整数$a$，满足$gcd(p,a)=1$，则有

$a^{p-1}≡1(mod p)$

• 二次探测定理：

对于$x∈(0,p)$，若$p$是素数，则方程：

$x^{2}≡1(mod p)$
的解为：

$x_1=1,x_2=p-1$

因为费马小定理的逆命题不成立，而逆否命题成立，所以我们可以利用一下一点：

对于任意整数$p>a$，不满足

$a^{p-1}≡1(mod p)$
则p为合数。

所以我们可以不断在区间$[2,p-1]$范围内随机取$a$，并进行判定。在$s$次判定不为合数之后，我们就可以说这个数是质数。
但是这还不够精确，我们可以先把$p-1$分解成

$2^{t}×u(u∈{x|x=2k+1,k∈N})$ 的形式
然后令 $x[0]=a^u mod p$
那么将 $x[0]$平方 $t$次就是(a^u)^2^t mod p的值
我们设$x[i]$为 $x[0]$平方i次的值，根据二次探测定理:
若 $x[i]$等于 $1$，则$x[i-1]$等于 $1$或$p-1$，不满足则 $p$为合数。
注意以上操作中所有的形如$a^bmodp$的式子都要用快速幂运算，当 $n$比较大时，形如$a*bmodp$的式子也要使用分治的思想来计算。

这就是Miller-Rabin算法的主要内容。
时间复杂度：考虑常数后为O($s(log_2n)^{3}$)

代码如下：

const int MAXN = 65;
long long n, x[MAXN];

long long multi(long long a, long long b, long long p) {
    long long ans = 0;
    while(b) {
        if(b&1LL) ans = (ans+a)%p;
        a = (a+a)%p;
        b >>= 1;
    }
    return ans;
}

long long qpow(long long a, long long b, long long p) {
    long long ans = 1;
    while(b) {
        if(b&1LL) ans = multi(ans, a, p);
        a = multi(a, a, p);
        b >>= 1;
    }
    return ans;
}

bool Miller_Rabin(long long n) {
    if(n == 2) return true;
    int s = 20, i, t = 0;
    long long u = n-1;
    while(!(u & 1)) {
        t++;
        u >>= 1;
    }
    while(s--) {
        long long a = rand()%(n-2)+2;
        x[0] = qpow(a, u, n);
        for(i = 1; i <= t; i++) {
            x[i] = multi(x[i-1], x[i-1], n);
            if(x[i] == 1 && x[i-1] != 1 && x[i-1] != n-1)
                return false;
        }
        if(x[t] != 1) return false;
    }
    return true;
}