Miller Rabin质数判定模板（2）

 

首先介绍一个相关的引理。

  

和

  

总是得到1，称这两个数为1的“平凡平方根”。当p是素数且p>2时，不存在

  

的“非平凡平方根”。为了证明该引理，我们假设x是

  

的平方根，于是有 [1] 

推出p|(x+1)(x-1)，也就是说，p|x+1或p|x-1。

现在假设n是一个奇素数，且 n>2。于是n-1是一个偶数，可以被表示为

  

的形式，s和d都是正整数且d是奇数。对任意在

  

范围内的a 和

  

，必满足以下两种形式的一种：

因为由于费马小定理，对于一个素数n，有

又由于上面的引理，如果我们不断对

  

取平方根，我们总会得到 1 或 -1。如果我们得到了 -1，意味着②式成立，n是一个素数。如果我们从未得到 -1，那么通过这个过程我们已经取遍了所有2的幂次，即①式成立。

米勒-拉宾素性测试就是基于上述原理的逆否，也就是说，如果如果我们能找到这样一个a，使得对任意

  

以下两个式子均满足：

那么n就不是一个素数。这样的a称为n是合数的一个凭证（witness）。否则a可能是是一个证明n是素数的“强伪证”（strong liar），即当n确实是一个合数，但是对当前选取的a来说上述两个式子均不满足，这时我们认为n是基于a的大概率素数。

每个奇合数n都有很多个对应的凭证a，但是要生成这些a并不容易。当前解决的办法是使用概率性的测试。我们从集合

  

中随机选择非零数a，然后检测当前的a是否是n为合数的一个凭证。如果n本身确实是一个合数，那么大部分被选择的a都应该是n的凭证，也即通过这个测试能检测出n是合数的可能性很大。然而，仍有极小概率的情况下我们找到的a是一个“强伪证”（反而表明了n可能是一个素数）。通过多次独立测试不同的a，我们能减少这种出错的概率。

对于测试一个大数是否是素数，常见的做法随机选取基数a，毕竟我们并不知道凭证和伪证在 [1,n-1]这个区间如何分布。典型的例子是 Arnault 曾经给出了一个397位的合数n，但是所有小于307的基数a都是n是素数的“强伪证”。不出所料，这个大合数通过了 Maple 程序的isprime()函数（被判定为素数）。这个函数通过检测 a=2,3,5,7,11这几种情况来进行素性检验。

//#include<pch.h>
#include <iostream>
#include <cstdio>
#include <bits/stdc++.h>
#include <map>
#include <algorithm>
#include <stack>
#include <iomanip>
#include <cstring>
#include <cmath>
#define DETERMINATION main
#define lldin(a) scanf_s("%lld", &a)
#define println(a) printf("%lld\n", a)
#define reset(a, b) memset(a, b, sizeof(a))
const int INF = 0x3f3f3f3f;
using namespace std;
const double PI = acos(-1);
typedef long long ll;
typedef unsigned long long ull;
typedef long double ld;
const int mod = 1000000007;
const int tool_const = 19991126;
const int tool_const2 = 2000;
inline ll lldcin()
{
	ll tmp = 0, si = 1;
	char c;
	c = getchar();
	while (c > '9' || c < '0')
	{
		if (c == '-')
			si = -1;
		c = getchar();
	}
	while (c >= '0' && c <= '9')
	{
		tmp = tmp * 10 + c - '0';
		c = getchar();
	}
	return si * tmp;
}
///Untersee Boot IXD2(1942)
/**Although there will be many obstructs ahead,
the desire for victory still fills you with determination..**/
/**Last Remote**/
const int primes[7] = { 2,3,5,31,37,211,317 };
ll quickm(ll a, ll b, ll c)
{
	return ((a*b - ll((ld)a / c * b + 0.005)*c) % c + c) % c;//有风险的快速乘,加上一个数使得它误差减小
}
ll quickp(ll a, ll b, ll c)
{
	ll ans = 1;
	while (b)
	{
		if (b & 1)
			ans = quickm(ans, a, c);
		b >>= 1;
		a = quickm(a, a, c);
	}
	return ans % c;
}
bool Miller_Robin(ll num)
{
	if (num == 2)
		return true;
	if ((num &1)==0)
		return false;
	ll tmp = num - 1,tmp5=0;
	while (!(tmp & 1))
	{
		tmp >>= 1;
		++tmp5;
	}
	for (int i = 0; i < 7; i++)
	{
		if (primes[i] == num)
			return true;
		ll solution = quickp(primes[i], tmp, num);//a^tmp=solution（mod num）
		ll tmpsolution = solution;
		for (int j = 1; j <= tmp5; j++)
		{
			solution = quickm(solution, solution, num);
			if (solution == 1 && tmpsolution != 1 && tmpsolution != num - 1)//如果二次检测原则两个式子都不能满足
				return false;
			tmpsolution = solution;//迭代
		}
		if (solution != 1)//回到了费马小定理逆定理，如果不是1，直接判否
			return false;
	}
	return true;
}
int DETERMINATION()
{
	ios::sync_with_stdio(false);
	cin.tie(0), cout.tie(0);
	ll x;
	while (cin >> x)
	{
		if (x == 1)
			cout << "N" << endl;
		else
		{
			if (Miller_Robin(x))
				cout << "Y" << endl;
			else
				cout << "N" << endl;
		}
	}
	return 0;
}