Miler-Rabbin素数判定

前言

素数判定？
小学生都可以打的出来！
直接暴力$O(\sqrt n)$……
然后就会发现，慢死了……
于是我们想到了筛法，比如前几天说到的詹欧筛法。
但是线性的时间和空间成了硬伤……如果是long long范围内的数，就搞不出来了。
那么素数判定就止步于此了吗？
不可能的，优化永无止境。我们可以用正确率来换时间啊！
Miller-Rabbin素数判定法就是这样的一个水法好方法。

---

前置数论知识

“费马小定理的逆定理”

首先是人人皆知的费马小定理：如果$p$为素数，对于任意非零整数$a$，满足$a^{p-1}\equiv 1(\mod p)$
怎么证明呢？抱歉，我初三了也不会证……（WHH五年级时就会证了……）
接下来有个猜想，如果有个正整数$p$，满足对于任意$a$，有$a^{p-1}\equiv 1( \mod p)$，那么$p$是素数。
这个猜想是正确的吗？
当然不是……
有种可恶的数叫Carmichael数，你可以将它称为伪素数。
这种恶心的数啊，它不是素数，但是它满足这条性质。
它的个数是无限的，但是分布很稀疏，在$10^8$内，只有$255$个Carmichael数。
所以正确率还是很高的……
可是我们要精益求精，并且，要知道这个世界上有很多毒瘤出题人……有时候数据怎样不取决于你的运气，全在于出题人……

二次探测定理

为了提高上面那东西的正确率，$Miller$和$Rabbin$决定引入二次探测定理。
定理内容很简单：若$p$是一个素数，方程$x^2\equiv 1(\mod p)$的唯二解为$p=\pm1(\mod p)$
证明？
移项得$x^2-1\equiv 0(\mod p)$，
所以$(x-1)(x+1)\equiv 0 (\mod p)$
所以$p \mid (x+1)(x-1)$
因为$p$是个素数，所以$p \mid (x+1)$或$p \mid (x-1)$
所以$x=\pm1(\mod p)$

我们反过来看，如果存在非零整数$x$使得这个定理不成立，那么$p$一定是合数。
所以这个东西可以大大地加强正确率。

---

Miller-Rabbin素数判定

这个算法就是将上面两个数论知识结合起来。
对于待测数$p$，先将$p-1$分解成$2^kt$的形式，其中$t$为奇数。
然后取几个基数$a$，分别做以下操作：
首先算出$a^t$，判断它是否为$\pm1$，如果是，则暂时判定为素数，退出。（这样它在经过自乘后会保持$1$，满足性质）
如果不是，就枚举$t-1$次，不断自乘，判断他是否为$-1$，如果是，则暂时判定为素数，退出。（原因类似，不用判断$1$，因为如果它是$1$，在之前必定是$\pm 1$，早就退出了）
搞完之后，如果之前没有退出，那么就将它判定为合数（再自乘一次就变成$a^{p-1}$，由于前面不是$\pm1$，所以一定不成立），退出。
所有基数计算完毕，如果判定为素数，就返回素数。

可以想象一下，从一个杂乱的数，变成$-1$，再变成$1$，后面都是$1$。
具体可以见代码。超短，特别好理解。

---

如何拥有更高的正确率？

前面说过Miller_Rabbin是牺牲正确性的算法。
所以基数的取值会和程序的效果有很大关系。
一般来说，可以取一堆随机数，这样就可以达到很高的正确率了。
但是我们要精益求精。
有一种方法是取前几个素数：

图片截自64位以内Rabin-Miller_强伪素数测试和Pollard_rho_因数分解算法的实现.doc

假如选取素数2 3 5 7，在$2.5*10^{13}$内唯一一个强伪素数为$3,215,031,751$。
假如选取素数2 3 7 61 24251，在$10^{14}$内唯一一个强伪素数为$46,856,248,255,981$。
背下来就可以了……（以防毒瘤出题人）

我脑子不好，所以枚举前面的几个素数。
虽然说Miller_Rabbin牺牲了正确率，但是在一定范围内，只要你的基数取得好，那么正确率是可以达到$100\%$的。

还有在判定之前，先拿几个素数来判一下，判完再走。前面的$7$个素数可以判掉$81.9\%$的数。
当然不要判太多，判越多次，比起先前贡献的增长就越小。

---

代码

代码比较简单，主体部分很短。
记得配上快速幂和龟速乘，不然会爆炸。

long long mul(long long a,long long b,long long mo){
	long long res=0;
	for (;b;b>>=1,a=(a<<1)%mo)
		if (b&1)
			res=(res+a)%mo;
	return res; 
}
inline long long mpow(long long x,long long y,long long mo){
	long long res=1;
	for (;y;y>>=1,x=mul(x,x,mo))
		if (y&1)
			res=mul(res,x,mo);
	return res;
}
const int p[7]={2,3,5,7,11,13,17};
inline bool mr(long long x){
	if (x==0 || x==1)
		return 0;
	for (int i=0;i<7;++i){
		if (x==p[i])
			return 1;
		if (!(x%p[i]))
			return 0;
	}
	int k=0;
	long long y=x-1;
	while (!(y&1))
		y>>=1,++k;
	for (int i=0;i<7;++i){
		long long s=mpow(p[i],y,x);
		if (s==1 || s==x-1)
			continue;
		for (int j=1;j<k;++j){
			s=mul(s,s,x);
			if (s==x-1)
				break;
		}
		if (s!=x-1)
			return 0;
	}
	return 1;
}