质数检验有不少算法，一般使用的质数检验复杂度是 $O(\sqrt{n})$ ；

又如线性筛可以在 $O(n)$ 的时间内求出所有1~n的质数

但是，当n非常大，连 $O(\sqrt{n})$ 的复杂度也难以接受时，上述算法便不能满足要求

这篇blog记录了一些关于Miller_Rabin算法的内容

大家都知道的著名的费马小定理：

$a^{p-1}\equiv1\pmod p$

其中 $a,p$ 互质

我们猜想，任意选取 $a$ ，如果一个数 $p$ 满足以上式子，那么它就很有可能是一个质数

但是这个猜想很容易找到反例： $a=2$ ， $p=341,561, 645, 1105$ 时，费马小定理逆命题不成立，人们把 $a^{p-1}\equiv1\pmod p$ 的合数称为以 $a$ 为底的“伪素数”

在 $[1,10^9]$ 中，质数一共有 $50847534$ 个，而满足 $2^{p-1}\equiv1\pmod p$ 的合数 $p$ 有 $5597$ 个，算法出错的概率太高了

一个想法是，同时使用多个a来进行判断，例如同时检验 $a=2,3$

在 $[1,10^9]$ 中，同时以 $2,3$ 为底的伪素数只有 $1272$ 个，不到以2为底的 $\frac{1}{4}$

选取的 $a$ 越多时，算法越准确，这便是费马素性检验

（以下内容引自Matrix67的博客）

人们自然会想，如果考虑了所有小于n的底数a，出错的概率是否就可以降到0呢？没想到的是，居然就有这样的合数，它可以通过所有a的测试（这个说法不准确，详见我在地核楼层的回复）。Carmichael第一个发现这样极端的伪素数，他把它们称作Carmichael数。你一定会以为这样的数一定很大。错。第一个Carmichael数小得惊人，仅仅是一个三位数，561。前10亿个自然数中Carmichael数也有600个之多。Carmichael数的存在说明，我们还需要继续加强素性判断的算法。

由此观之，费马素性检验仍不能满足我们的要求，我们需要改进素性检验算法

引理：

对于 $\forall a<p$ ， $p$ 是质数

若 $a^2\equiv1\pmod p$

那么有 $a=1$ 或 $a=p-1$

证明：

$a^2\equiv1\pmod p$

$\Rightarrow a^2-1=(a+1)*(a-1)\equiv0\pmod p$

Miller_Rabin素数测试的方法是，对于待检验数 $x$ ，不断地提取 $x-1$ 中的因子 $2$ ，把 $x-1$ 表示成 $2^r*d$ 的形式，那么我们需要计算的东西变成了 $a^{2^r*d}\mod x$

于是，如果 $x$ 是质数， $a^{2^r*d}\mod x$ 要么等于 $1$ ，要么等于 $x-1$

如果 $a^{2^{r-1}*d}\mod x =1$ ，定理同样适用于 $a^{2^{r-2}*d}\mod x$

不断地开方，直到存在一个 $i\in[0,r)$ ，使得 $a^{2^i*d}\mod x=x-1$

至此，费马素性检验被强化为如下形式：

尽可能提取因子 $2$ ，将待检验数 $x$ 表示为 $2^r*d$ 的形式，其中 $d$ 是一个奇数
如果 $a^d\mod n=1$ ，或者找到一个 $i\in[0,r)$ ，使得 $a^{2^i*d}\mod x=x-1$ ，那么 $x$ 极有可能是一个质数，否则 $x$ 就不是一个质数

上文提到极有可能，是因为Miller_Rabin同样是不确定算法，我们称能够通过以 $a$ 为底的Miller_Rabin检验的合数称为以 $a$ 为底的“强伪素数”，第一个以 $2$ 为底的强伪素数为 $2047$ ，而第一个以 $2$ 和 $3$ 为底的强伪素数则为1373653

如果选用 $2,3,7,61$ 和 $24251$ 作为 $a$ ，那么 $10^{16}$ 内唯一的强伪素数为 $46\space856\space248\space255\space981$ ，正确率可以接受

Miller_Rabin核心代码

这里，我采用先将n-1的所有因子2提取出来，再不断地平方回到原来的n-1的枚举方法

typedef long long ll;
const ll a[10]={0,2,3,7,61,24251};

ll fastpow(ll a,ll b,ll p)
{
	ll re=1,base=a;
	while(b)
	{
		if(b&1)
			re=re*base%p;
		base=base*base%p;
		b>>=1;
	}
	return re;
}

bool Miller_Rabin(ll x,ll a)
{
	if(x==a)
		return true;
	if(x%a==0 || x<2)
		return false;
	ll d=x-1,r=0;
	while(!(d&1))
	{
		d>>=1;
		++r;
	}
	ll k=fastpow(a,d,x);
	if(k==1 || k==x-1)
		return true;
	for(ll i=1;i<r;++i)
	{
		k=k*k%x;
		if(k==x-1)
			return true;
	}
	return false;
}

bool test(ll x)
{
	for(int i=1;i<=5;++i)
		if(!Miller_Rabin(x,a[i]))
			return false;
	return true;
}