MillerRabin学习笔记

\(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");
	}
}