Miller-Rabin素性测试

可以去看本文的大部分来源：数论部分第一节：素数与素性测试

在判断一个数是否为质数时，我们通常选择根号 n 的试除法来进行判断

但这是在询问比较少/数据比较小的前提下的

Miller-Rabin 素性测试为我们提供了一种 logn 的做法

我们知道 Fermat 小定理 a^(p - 1) mod p = 1 是对于 p 为质数且 a，p 互质的情况成立的

那么其逆命题成立吗？

过去人们是这样认为的

但由于伪素数的发现，人们终于认识到它是错误的

不过这样的数还是占据少数的

所以我们至少可以用它来先验证出一部分数，尽管这个集合中有伪素数

可以通过多取一些 a 的值来降低错误概率

以上是Fermat素性测试。

那么我们自然会这样想，如果我们考虑了所有的小于 n 的底数 a ，是否就足够区分素数和其他数呢？

很可惜，这样的合数依然存在

他们被称为“Carmichael数”

你一定会认为这样极端的伪素数一定很大

然而第一个 Carmichael数仅仅是一个三位数：561

而且前 10 亿个自然数中 Carmichael数有600个之多

Carmichael数的存在说明，我们还需要继续加强素性判断的算法。

Miller和Rabin两个人的工作让Fermat素性测试迈出了革命性的一步，建立了传说中的Miller-Rabin素性测试算法。

新的测试基于下面的定理：

　　如果 p 是素数，x 是小于 p 的整数，且 x^2 mod p = 1 ，那么要么 x = 1 ，要么 x = p - 1，这是显然的，因为 x^2 mod p = 1 相当于 p 能整除 x^2 - 1，也即p能整除 (x + 1)(x - 1) 。由于 p 是素数，那么只可能是 x - 1 = 0(mod p) (此时 x = 1 )或 x + 1 = 0(mod p) 整除(此时 x = p - 1 )。

下面用 341 来演示一下上面的定理如何运用到 Fermat定理上：

 　　2^340 mod 341 = 1  => (2^170)^2 mod 341 = 1

　　若 341 为素数，则上面方程的解可能是 1 或 340，若等于 1，我们继续上面的变形得：

　　(2^85)^2 mod 341 = 1

　　然而计算得 2^85 mod 341 = 32

　　那么 341 不是素数

这就是 Miller-Rabin 素性测试的方法，通过不断提取指数 p - 1中的因子 2 ，将方程化为 x^2 mod p = 1 的形式， 根据解 只能为 1 或 p - 1 ，若解为 1，继续重复上述过程，解为 p - 1，停止，这个数目前为止我们认为它可能是一个素数，解为其他数，则这个数不是素数

时间复杂度 O(logn)

Miller-Rabin素性测试同样为不确定算法，我们把可以通过以a为底的Miller-Rabin测试的合数称作以a为底的强伪素数(strong pseudoprime)。第一个以2为底的强伪素数为2047。第一个以2和3为底的强伪素数则大到1 373 653。

以上是 Miller-Rabin 素性测试算法。

通常我们使用该算法时，都是要来判断 long long 级别的数字，由于可能在快速幂里乘法的步骤乘爆了，需要快速乘， O(logn)，写起来跟快速幂差不多

附上一道例题：hihocoder 1287 数论一·Miller-Rabin质数测试

代码：

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cctype>
#include<cstdio>
using namespace std;
 
typedef long long ll;
 
int m, n;
ll ps[12] = {2ll, 3ll, 5ll, 7ll, 11ll, 13ll, 17ll, 19ll, 23ll, 29ll, 31ll, 37ll};
 
inline ll turtle_mul(ll a, ll b, ll mod) {
    ll ans = 0ll;
    while(b) {
        if(b & 1) ans = (ans + a) % mod;
        a = (a << 1) % mod;
        b >>= 1;
    }
    return ans;
}
inline ll fast_pow(ll low, ll up, ll mod) {
    ll ans = 1ll;
    low %= mod;
    while(up) {
        if(up & 1ll) ans = turtle_mul(ans, low, mod);
        up >>= 1;
        low = turtle_mul(low, low, mod);
    }
    return ans;
}
inline bool dvd_chk(ll bot, ll top, ll mod) {
	register ll tmp;
	while (!(top & 1)) {
		tmp = fast_pow(bot, top, mod);
		if (tmp == 1ll) top >>= 1;
		else if (tmp == mod - 1ll) return true;
		else return false;
	}
	return true;
}
inline bool Test(ll a) {
    if(a <= 1) return false;
    if(a == 2ll) return true;
    if(!(a & 1ll)) return false;
    for(int i = 0; i < 12; ++i) {
        if(ps[i] == a)
            return true;   
        if(fast_pow(ps[i], a - 1ll, a) != 1ll)
            return false;  
        else if(!dvd_chk(ps[i], a - 1ll, a))
            return false;
    }
    return true;
}
 
int main() {
    scanf("%d", &m);
    ll num;
    while(m--) {
        scanf("%lld", &num);
        if(Test(num)) puts("Yes");
        else puts("No");
    }
    return 0;
}