Miller-Rabin学习笔记

$Miller-Rabin$

$Miller-Rabin$ 用于判定一个大整数是不是素数，且速度非常快
应该是 $O(klog^3n)$，其中 $k$ 为测试的次数，$n$ 为要判定的数
算法本质上是一种概率算法，存在误判的可能性，但是出错的概率非常小。出错的概率到底是多少，存在严格的理论推导。（网上copy的）

两个重要的东西

$Farmet$测试

根据费马小定理，如果 $p \in \mathbb P$ 且 $\gcd(a , p) = 1$ 那么 $a^{p-1} \equiv 1 \pmod p$。它给我们提供了一种检验素数的思路：对于数n，它的基本思想是不断地选取在 $[2,n-1]$ 的数让它为 $a$ 并检验是否每次都有 $a^{n-1} \equiv 1 \pmod n$，如果没有，那么 $n$ 肯定不是素数

其实有它我们就可以大概率判定素数了
为什么是大概率呢？
谁告诉你费马小定理的逆定理是成立的？
也就是说费马小定理的逆定理并不成立，换言之，满足了 $a^{n-1} \equiv 1 \pmod n$，$n$ 也不一定是素数。

例子：1819年有人发现了费马小定理逆命题的第一个反例：虽然2的340次方除以341余1，但341=11*31，后来，人们又发现了561, 645, 1105等数都表明a=2时Fermat小定理的逆命题不成立。虽然这样的数不多，但不能忽视它们的存在。统计表明，在前10亿个自然数中共有50847534个素数，而满足 $2^{n-1} ≡ 1 \pmod n$ 的合数n有5597个
如果用费马小定理的逆命题来判断一个正整数n是不是素数，在前10亿个自然数中出错的可能性为 0.011%
这个出错的可能性还是很高的，但是仍然可以用这个技巧来排除大量的合数，这种方法就是费马检测

上面的做法中随机地选择 $a$ ，很大程度地降低了犯错的概率。但是仍有一类数，上面的做法并不能准确地判断。
对于合数 $n$ ，如果对于所有正整数 $a$ ， $a$ 和 $n$ 互素，都有同余式 $a^{n-1} ≡ 1 \pmod n$ 成立，则合数 $n$ 为卡迈克尔数（$\texttt{Carmichael Number}$），又称为费马伪素数。
比如， $561 = 3 \times 11 \times 17$ 就是一个卡迈克尔数。
而且我们知道，若 $n$ 为卡迈克尔数，则 $m = 2^n-1$ 也是一个卡迈克尔数，从而卡迈克尔数的个数是无穷的。

于是 $Farmet$ 测试变得很玄学

当然，因为卡迈克尔数的数量在一定范围内是很少的，$1 \sim 10^8$ 内只有255个，如果提前打好一个卡迈克尔数的表来特判，还是可以的
但不是主流
于是我们可以愉快请出第二位 $boss$

二次探测定理

根据有限域上的平方根定理：$p$为奇素数（$p \ne 2$ 的素数），有 $x^2 \equiv 1 \pmod {p^e}$， $x$ 仅有两根 $x=1$ 或 $x=-1$，在模 $p$ 的意义下 $x=-1$ 相当于 $x=p-1$，我们称 $\pm 1$ 为1的平凡平方根。很容易有一个推论，如果模n存在1的非平凡平方根，n一定是合数

当然，为什么仅有两根 $\pm 1$ ？一个因式分解就好了

$$\begin{aligned} x^2 \equiv 1 \pmod {p^e} \\ x^2 -1 \equiv 0 \pmod {p^e} \\ (x+1)(x-1) \equiv 0 \pmod {p^e} \end{aligned}$$

可解出 $x = 1$ 或 $x=-1$

根据卡迈克尔数的性质，可知其一定不是 $p^e$
不妨将费马小定理和二次探测定理结合起来使用：

先使用二次探测定理，将 $a^{n-1}$ 分成 $a^{2^b*d}$，可以从 $x = a^d$ 开始，依次平方 $b$ 次，每次平方的时候模上 $n$，按照之前的平方根定理，如果模上 $n$ 的结果为 $1$ 的话，那么 $x$ 一定是 $1$，或者是 $n-1$，如果不满足则不是素数，$x=x^2$ ，再次循环。最后，我们其实算出了 $a^{n-1}$ 模 $n$ 意义的值，再用费马小定理判断一下就好了

当然，因为是大整数，所以在算两数相乘时极有可能爆 $\texttt{long long}$，所以我们用上快速乘
具体细节看代码

$Code$

#include<cstdio>
#include<ctime>
#include<algorithm>
using namespace std;
typedef long long LL;
 
inline LL fmul(LL x , LL y , LL p) //快速乘
{
	LL res = 0;
	while (y)
	{
		if (y & 1) res = (res + x) % p;
		y >>= 1 , x = (x + x) % p;
	}
	return res;
}
 
inline LL fpow(LL x , LL y , LL p) //快速幂
{
	LL res = 1;
	while (y)
	{
		if (y & 1) res = fmul(res , x , p);
		y >>= 1 , x = fmul(x , x , p);
	}
	return res;
}
 
 
inline int Miller_Rabin(LL m , int test_time) //test_time 测试的次数，每次选取随机的在[2..n-1]范围内的 a
{
	if (m < 3) return (m == 2);
	if (!(m & 1)) return 0;
	
	LL d = m - 1;
	int b = 0;
	while (!(d & 1)) d = d >> 1 , b++; //将 a^{n-1} 分成 a^{2^b*d}
	
	srand((unsigned)time(NULL));
	for(register int i = 0; i < test_time; i++)
	{
		LL a = rand() % (m - 2) + 2 , v = fpow(a , d , m) , u = 0;
		for(register int j = 0; j < b; j++)
		{
			u = fmul(v , v , m); //不断平方
			if (u == 1 && v != 1 && v != (m - 1)) return 0; //平方根定理推论判定
			v = u;
		}
		if (v != 1) return 0;
	}
	return 1;
}
 
int main()
{
	int n;
	scanf("%d" , &n);
	LL m;
	while(n--)
	{
		scanf("%lld" , &m);
		if (Miller_Rabin(m , 6)) printf("Prime\n");
		else printf("Not prime\n");
	}
}