Miller-Rabin算法

算法背景

米勒-拉宾素性检验（Miller Rabin算法），是一种素数判定法则，利用随机化算法判断一个数是合数还是可能是素数。卡内基梅隆大学的计算机系教授Gary Lee Miller首先提出了基于广义黎曼猜想的确定性算法，由于广义黎曼猜想并没有被证明，其后由以色列耶路撒冷希伯来大学的Michael O. Rabin教授作出修改，提出了不依赖于该假设的随机化算法。

涉及知识点

费马小定理
费马小定理就是说如果p是一个质数，而整数a不是p的倍数，则有a^（p-1）≡1（mod p）。也可以写做a^p ≡ p (mod p).

二次探测
如果p是一个素数，那么使得x² ≡ 1 (mod p)的 x的解只有两种可能，就是x = 1 或者 x = p-1。
证明如下：
x² ≡ 1 (mod p)可以转换为 x²-1 = 0（mod p），然后就化为了(x-1)(x+1) = 0 (mod p), 因此p是整除(x-1)(x+1), 用数学语言也就是说p|(x-1)(x+1)。因此x = 1 或者 x = p-1。

Miller-Rabin证明过程

假设n为奇素数，那么满足方程x²≡1（mod n）的解为x = 1 或者 x = n-1。
又已知费马小定理 a^n-1≡1（mod n）

已知n为奇数，那么n-1一定为偶数，我们可以设n-1 = m * 2^t。其中m一定为奇数，否则就可以并到2^t里面。其实原本是n-1 = m * q^t。为了计算就给q特殊化为2了。那么如何求m和t呢，上面其实也说了，m一定是奇数，所以只要给n-1不断地除2，直到n-1变为奇数，那么m就是这个奇数，设初始t为0，每次(n-1)除以2后，t都加一。这样t和m都求出来了。

接下来我们从[1, n-1]里面随机取一个a，容易得到a^n-1 = a^(m*(2 ^ t)), 又由二次探索可知

不断的递归下去，不断的用二次探索定理，直到最后a^m = 1(mod n) 或者a^m = n-1(mod n)，中间只要有一个不满足二次测试定理，那么这个数就一定是合数，否则就有3/4的概率为质数。这样其实为合数的几率也不小，有25%的大小。因此我们要多取几次a进行计算，这样是合数的概率就变为了（1/4）^k。最后其实也就是相当与肯定是质数了。

算法流程

上面说了整体的证明，可能有一些人看到代码还是蒙，那我们就说一下算法的流程。

1. 首先对于要判断的数n，先判断是不是2，是就直接返回true。
2. 判断是不是小于2或者为偶数，是就返回false。
3. 令n-1 = m*2^t，求出u和t。
4. 利用rand随机找一个a， a属于[1， n-1]
5. 令x = (a^m)%n
6. 然后进行t次循环，为什么要循环t次呢，上面证明是从2^t倒着推到2⁰，我们写程序完全可以正着来，上面的x其实就是
   
7. 那么在循环中每次都进行 y =（x*x）%n， x = y 的操作，最终会变为a^n-1，也就是
   
   又因为x在循环时是一定小于n的，所以可以用二次探测定理
   运算过程中如果(x²)%n = 1, 也就是y = 1，假如n是素数，那么x == 1 或者x == n-1，否则n就一定不是素数，直接返回false。
8. 整个循环结束后，程序运行到最后x = a^n-1, 根据费马小定理，如果x还是不等于1，那么肯定不是素数，直接返回false。
9. Miller-Rabin算法正确率为75%，所以要进行大约4~8次来提高正确率。

下面具体操作看代码

#include <iostream>
#include <cstdlib> 

using namespace std;

typedef long long ll;

ll qul_mul(ll a, ll b, ll n)  //快速积运算，快速求出a^m^*a^m^并且能防止爆掉
{
	ll num = 0;
 
 	while (b)
	{
		if (b&1)
		{
			num=(num+a)%n;	
		}
	  
		a = (a+a)%n;
		b>>=1;
	}
 
	return num;
}

ll qul_pow(ll a, ll b, ll n) //快速幂计算快速求出a^m^的值
{
	ll sum = 1;
 
	while (b)
	{
		  if (b&1)
		  {
		  	sum = sum * a % n;
		  }
		  
		  a=a*a%n;
		  b>>=1;
	}
		 
	 return sum;
}

bool Miller_Rabin(ll n)
{
	int m = n-1;
	int t = 0;
 
	if (n==2)
	{
 		return true;
	}
	else if (n<2 || !(n&1))
	{
 		return false;
	}
 
 	//求出n-1 = m*2^t的m和t。
	while (!(m&1))
	{
		m>>=1;
		t++;
	}
 
	for (int i=0; i<20; i++)
	{    
		//随机数a
		ll a = rand()%(n-1) + 1;
		// 计算a^m
		ll x = qul_pow(a, m, n), y;
  		
		for (int j=0; j<t; j++)
		{
			y = qul_mul(x, x, n);  //进行(x*x)%n操作。
   			
   			//不满足二次探测定理，也就是y得1了但是x并不等于1或者n-1，那么n就一定不是质数。
			if (y==1 && x!=1 && x!=n-1)
			{
 				return false;
			}
   
			x = y;
		}
  		//上面循环跑完了，如果最后x都不等于1，那么一定是一个合数了。
		if (x!=1)
		{
			return false;
		}
	}
 
	return true;
}

int main()
{
	ll n;
 
	while (cin >> n)
	{
		cout << Miller_Rabin(n) << endl;
	}
 
	return 0;
}

例题推荐：
hdu2138