hihoCoder 1287 Miller-Rabin质数测试

Miller-Rabin质数测试

摘自hihoCoder：

这种质数算法是基于费马小定理的一个扩展，首先我们要知道什么是费马小定理：

费马小定理：对于质数p和任意整数a，有a^p ≡ a(mod p)(同余)。反之，若满足a^p ≡ a(mod p)，p也有很大概率为质数。
将两边同时约去一个a，则有a^(p-1) ≡ 1(mod p)

也即是说：假设我们要测试n是否为质数。我们可以随机选取一个数a，然后计算a^(n-1) mod n，如果结果不为1，我们可以100%断定n不是质数。

否则我们再随机选取一个新的数a进行测试。如此反复多次，如果每次结果都是1，我们就假定n是质数。

该测试被称为Fermat测试。需要注意的是：Fermat测试不一定是准确的，有可能出现把合数误判为质数的情况。

Miller和Rabin在Fermat测试上，建立了Miller-Rabin质数测试算法。

与Fermat测试相比，增加了一个二次探测定理：

如果p是奇素数，则 x^2 ≡ 1(mod p)的解为 x ≡ 1 或 x ≡ p - 1(mod p)

如果a^(n-1) ≡ 1 (mod n)成立，Miller-Rabin算法不是立即找另一个a进行测试，而是看n-1是不是偶数。如果n-1是偶数，另u=(n-1)/2，并检查是否满足二次探测定理即a^u ≡ 1 或 a^u ≡ n - 1(mod n)。

举个Matrix67 Blog上的例子，假设n=341，我们选取的a=2。则第一次测试时，2^340 mod 341=1。由于340是偶数，因此我们检查2^170，得到2^170 mod 341=1，满足二次探测定理。同时由于170还是偶数，因此我们进一步检查2^85 mod 341=32。此时不满足二次探测定理，因此可以判定341不为质数。

将这两条定理合起来，也就是最常见的Miller-Rabin测试。

但一次MR测试仍然有一定的错误率。为了使我们的结果尽可能的正确，我们需要进行多次MR测试，这样可以把错误率降低。

写成伪代码为：

Miller-Rabin(n):
	If (n <= 2) Then
		If (n == 2) Then
			Return True
		End If
		Return False
	End If
	
	If (n mod 2 == 0) Then
		// n为非2的偶数，直接返回合数
		Return False
	End If
	
	// 我们先找到的最小的a^u，再逐步扩大到a^(n-1)
	
	u = n - 1; // u表示指数
	while (u % 2 == 0) 
		u = u / 2
	End While // 提取因子2
	
	For i = 1 .. S // S为设定的测试次数
		a = rand_Number(2, n - 1) // 随机获取一个2~n-1的数a
		x = a^u % n
		While (u < n) 
			// 依次次检查每一个相邻的 a^u, a^2u, a^4u, ... a^(2^k*u)是否满足二次探测定理
			y = x^2 % n 
			If (y == 1 and x != 1 and x != n - 1)	// 二次探测定理
				// 若y = x^2 ≡ 1(mod n)
				// 但是 x != 1 且 x != n-1
				Return False
			End If
			x = y
			u = u * 2 
		End While
		If (x != 1) Then	// Fermat测试
			Return False
		End If
	End For
	Return True

QAQ：（本渣太弱了，并没有实现这种伪代码（实现后，wrong了），哪位聚聚看到后，帮助弱一下。）

值得一提的是，Miller-Rabin每次测试失误的概率是1/4；进行S次后，失误的概率是4^(-S)。

时间复杂度：

每一次单独的MR测试，需要O(log n)的时间。一共要进行S次MR测试，也就是O(Slog n)。

AC代码：

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;

//****************************************************************
// Miller_Rabin 算法进行素数测试
//速度快，而且可以判断 <2^63的数
//****************************************************************
const int S = 20;//随机算法判定次数，S越大，判错概率越小

//计算(a*b)%c.   a,b都是ll的数，直接相乘可能溢出的
//a,b,c <2^63
ll mult_mod(ll a,ll b,ll c){
    a %= c;
    b %= c;
    ll res = 0;
    while(b){
        if(b & 1){
            res += a;
            res %= c;
        }
        a <<= 1;
        if(a >= c)
            a %= c;
        b >>= 1;
    }
    return res;
}

//计算x^n %c
ll pow_mod(ll x,ll n,ll mod){//x^n%c
    if(n == 1)
        return x%mod;
    x %= mod;
    ll tmp = x;
    ll res = 1;
    while(n){
        if(n & 1)
            res = mult_mod(res,tmp,mod);
            tmp = mult_mod(tmp,tmp,mod);
        n >>= 1;
    }
    return res;
}

//以a为基,n-1=x*2^t      a^(n-1)=1(mod n)  验证n是不是合数
//一定是合数返回true,不一定返回false
bool check(ll a,ll n,ll x,ll t){
    ll res = pow_mod(a,x,n);
    ll last = res;
    for(int i = 1; i <= t; ++i){
        res = mult_mod(res,res,n);
        if(res==1 && last !=1 && last != n-1)
            return true;//合数
        last = res;
    }
    if(res != 1)
        return true;
    return false;
}

// Miller_Rabin()算法素数判定
//是素数返回true.(可能是伪素数，但概率极小)
//合数返回false;
bool Miller_Rabin(ll n){
    if(n < 2) return false;
    if(n == 2) return true;
    if((n&1) == 0) return false;//偶数
    ll x = n-1;
    ll t = 0;
    while((x&1)==0){
        x >>= 1;
        ++t;
    }
    for(int i = 0; i < S; ++i){
        ll a = rand()%(n-1)+1;
        if(check(a,n,x,t))
            return false;//合数
    }
    return true;
}

int main(){
    int T;
    scanf("%d",&T);
    while(T--){
        ll n;
        scanf("%lld",&n);
        if(Miller_Rabin(n))
            puts("Yes");
        else
            puts("No");
    }
    return 0;
}