Miller-Rabin质数测试

这种质数算法是基于费马小定理的一个扩展。

费马小定理：对于质数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)

伪代码：

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

写成 C++ 代码：

// 快速判断是否为质数 
#include <iostream>
#include <time.h>
#include <algorithm>
#define SS 100
#define ll long long 
using namespace std;
 
ll mod_mul(ll a, ll b, ll n){
    ll res = 0;
    while (b){
        if(b & 1) res = (res + a) % n;
        
        a = (a + a) % n;
        b >>= 1;
    }
    return res;
}
// 快速幂 (a^b)%n 
ll mod_exp(ll a, ll b, ll n){
    ll res = 1;
    while(b){
        if(b & 1) res = mod_mul(res, a, n);
        
        a = mod_mul(a, a, n);
        b >>= 1;
    }
    return res;
}

bool MillerRabin(ll x){
    if(x == 2) return true;
    if(x % 2 == 0 || x < 2) return false;
    
    ll u = x-1;
    while(u%2 == 0) u = u/2;
    
    for(int i = 1; i <= SS; i++){
        srand((unsigned)time(NULL));
        ll a = (rand()%(x - 2)) + 2;
        ll xx = mod_exp(a, u, x);
        while(u< x){
            ll yy = mod_exp(xx, 2, x);
            if(yy == 1 && xx != 1 && xx != x-1) return false;
            
            xx = yy,u = u*2;
        }
        if(xx!=1) return false; 
    }
    return true;
}

int main(){
    int n;
    cin>>n;
    while(n--){
        ll number;
        cin >> number;
        cout << (MillerRabin(number)?"Yes":"No") << endl;
    } 
    return 0;
}