与数论的爱恨情仇--01：判断大素数的Miller-Rabin

　　在我们需要判断一个数是否是素数的时候，最容易想到的就是那个熟悉的o(√n)的算法。那个算法非常的简单易懂，但如果我们仔细想想，当n这个数字很大的时候，这个算法其实是不够用的，时间复杂度会相对比较高。

　　怎么解决呢？我们先来了解一下“费马小定理”。假设我们有一个素数p，且另一个数a和p互素，就可以得到a^p-1≡1(mod p)。这个定理很巧妙啊，有人就想了，能不能通过费马小定理来判断一个数是否是素数呢？也就是说，当我们判断一个数p是否是素数时，只需要判断a^p-1≡1(mod p)是否成立即可。这里的a因为是任意数，干脆就让它等于2，那么判断一个数p是否是素数就转化成了判断2^p-1≡1(mod p)是否成立。乍一看，这好像没什么问题。当这个式子不成立时，p一定是一个合数，这没毛病；但是当式子成立的时候p就一定是素数吗？我们举个反例。当p = 341时，2³⁴⁰≡1(mod 341)成立，然而很不巧，341 = 11 * 31是一个纯正的合数。这就是数学中所说的，对于所有的a都存在对应的伪素数。(ps:这个问题并不能完全通过改变基数解决)

　　那我们该怎么办呢？其实很简单，只需要进行“二次探测”把伪素数揪出来就ok了。这可不是乱改，是有依据的：当p为素数时，方程x²≡1(mod p)有两个根 x = 1 和 x = p - 1。这两个根被我们赋予了一个奇怪的名字：平凡平方根。那么，判断一个数p是否是素数这个问题就又被我们转化，变成了判断在模p意义下是否存在1的非平凡平方根，若存在则p为合数，反之则为素数。这一测试被我们“亲切地”称为“miller - rabin测试”。

　　具体操作步骤如下：

　　①选取多个基数a进行测试；

　　②寻找模p为1的非平凡平方根；令p - 1 = 2^t*u(t >= 1, u为奇数)，a^p-1 = a^2t*u = a^2t  ，先算出x=a^u (mod p)，再把 x 平方 t 次，每次模上 p，这样我们就得到了一个长度为 t + 1 的序列。我们希望这个序列以 1 结尾。若中间某一项为 1，则这一项的前一项必须为 1 或 n - 1，否则p就是合数。

　　在miller - rabin测试中进行s次测试，这并不代表这项测试是简单地验证费马小定理，它大大降低了出错的概率，研究表明，miller - rabin测试的出错概率至多为 2^-s 这可以说是非常小了。所以不用担心它的准确度和严谨性。

　　代码如下：(你没看错就是rand())

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

int pow(long long a, long long b, long long n) {
    long long ans = 1;
    while (b) {
        if (b & 1) {
            ans = ans * a % n;
        }
        a = a * a % n;
        b >>= 1;
    }
    return ans;
}

bool judge(long long n, long long a) {
    long long u = 0, t = n - 1;
    while (t % 2 == 0) {u++; t /= 2;}
    long long x = pow(a, t, n);
    
    for (int i = 1; i <= u; i++) {
        long long nx = x * x % n;
        if ((nx == 1) && (x != 1) && (x != n - 1))
            return true;
        x = nx;
    }
    if (x != 1) return true;
    return false;
}

bool miller(long long n, int s) {
    if (n == 2) return true;
    if (n < 2 || n % 2 == 0) return false;
    for (int i = 1; i <= s; i++) {
        long long a = rand() % (n - 2) + 2;
        if (judge(n, a)) return false;
    }
    return true;
}

int main() {
    int t;
    cin >> t;
    long long a;
    while (t--) {
        cin >> a;
        if (miller(a, 10)) {
            cout << "yes" << endl;
        }
        else cout << "no" << endl;
    }
    return 0;
}