Miller-Rabin素性测试

有时候我们想快速的知道一个数是不是素数，而这个数又特别的大导致 $O(\sqrt n)$ 的算法也难以通过，这时候我们可以对其进行 Miller-Rabin 素数测试，可以很大概率测出其是否为素数。

两个理论基础

（1）费马小定理：当 $p$ 为质数，有 $a^{p-1}\equiv 1(mod \ p)$，反过来不一定成立，也就是说，如果 $a, \ p$ 互质，且 $a^{p-1}\equiv 1(mod \ p)$，不能推出 $p$ 是质数，比如 $Carmichael$ 数

（2）二次探测：如果 $p$ 是素数，$0 < x < p$，则方程 $x^2 \equiv 1(mod \ p)$ 的解为 $x=1$ 或 $x=p-1$.

证一下二次探测定理：

易知 $ x^2-1 \equiv 0(mod \ p)$

$\because(x+1)(x-1) \equiv 0\ (mod \ p)$

$\because$ $p$ 为素数，只能分解为

$\because x+1 \equiv 0\ (mod \ p)$ 或 $x-1 \equiv 0\ (mod \ p)$

根据 $x$ 的取值范围，$x=1$ 或 $x=p-1$.

Fermat 素性测试

只根据费马小定理，易得一种检测质数的思路：

它的基本思想就是多次从 $\left [ 2, n-1 \right ]$ 中选取基 $a$，并检验是否每次都有 $a^{n-1} \equiv 1\ (mod \ n)$

bool millerRabin(int n) {
  if (n < 3) return n == 2;
  // test_time 为测试次数,建议设为不小于 8
  // 的整数以保证正确率,但也不宜过大,否则会影响效率
  int  test_time = 20;
  for (int i = 1; i <= test_time; ++i) {
    int a = rand() % (n - 2) + 2;
    if (qpow(a, n - 1, n) != 1)
    {
        printf("a:%d\n", a);
        return false;
    }
  }
  return true;
}

遗憾的是，由于费马小定理的逆定理并不成立，即满足 $a^{n-1} \equiv 1\ (mod \ n)$，$n$ 也不一定是素数。

卡迈克尔数

上面的做法中随机地选择 $a$，很大程度上降低了犯错概率，但是仍有一类数，上面的做法并不能准确地判断。

对于合数，如果对于所有与 $n$ 互素的正整数 $a$，都有同余式 $a^{n-1} \equiv 1\ (mod \ n)$ 成立，则合数 $n$ 为卡迈克尔数（Carmichael Number），又称费马伪素数。

而且我们知道，若 $n$ 为卡迈克尔数，则 $2^n-1$ 也是卡迈克尔数，从而卡迈克尔数的个数是无穷的。

在1~100000000范围内的整数中，只有255个Carmichael数。前3个Carmichael数是561,1105,1729。

Miller-Rabin测试

根据卡迈克尔数的性质，可知其一定不是 $p^e$.

不妨将费马小定理和二次探测定理结合起来使用：

对偶数和0、1、2直接判断
设要测试的数为 $n$，设$m$、$k$，满足 $n-1 = m2^k$（使其中 $m$ 为奇数）
我们先算出 $a^t$，然后不断地平方且进行二次探测（进行 $k$ 次）
最后根据费马小定理探测一次
多次取不同的 $a$ 进行上面3步

可以证明Miller-Rabbin算法单次出错的概率小于等于 $\frac{1}{4}$，若重复 $k$ 次，则出错概率降低至 ${(\frac{1}{4})}^k$。这是一个很保守的估计，实际效果要好得多。

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

typedef long long int ll;

ll mod_mul(ll a, ll b, ll mod)
{
    ll res = 0;
    while (b)
    {
        if (b & 1)
            res = (res + a) % mod;
        a = (a + a) % mod;
        b >>= 1;
    }
    return res;
}

ll mod_pow(ll a, ll n, ll mod)
{
    ll res = 1;
    while (n)
    {
        if (n & 1)
            res = mod_mul(res, a, mod);
        a = mod_mul(a, a, mod);
        n >>= 1;
    }
    return res;
}

// Miller-Rabin随机算法检测n是否为素数
bool Miller_Rabin(ll n)
{
    if (n == 2)
        return true;
    if (n < 2 || !(n & 1))
        return false;
    ll m = n - 1, k = 0;
    while (!(m & 1))
    {
        k++;
        m >>= 1;
    }
    for (int i = 1; i <= 10; i++)  // 20为Miller-Rabin测试的迭代次数
    {
        ll a = rand() % (n - 1) + 1;
        ll x = mod_pow(a, m, n);
        ll y;
        for (int j = 1; j <= k; j++)
        {
            y = mod_mul(x, x, n);
            if (y == 1 && x != 1 && x != n - 1)
                return false;
            x = y;
        }
        if (y != 1)
            return false;
    }
    return true;
}

int main()
{
    ll n;
    while(scanf("%lld", &n) == 1)
    {
        printf("%d\n", Miller_Rabin(n));
    }
}