【技能点】Miller_Rabin算法(大质数判定)

1. 一个数$p$如果是质数，那么一定满足费马小定理，$x^{p-1}\equiv 1(modp)$；
2. 对于能满足$条件1$的数$p$，显然不一定都是质数。所以需要进行所谓的二次探测。
   即若$p$是质数且$x^2\equiv 1(modp)$，则$x\equiv 1(modp)$或$x\equiv p-1(modp)$。

二次探测证明如下：
$$\begin{array}{rl}{x}^2\equiv 1(& modp)\\ x^2-1\equiv & 0(modp)\\ (x-1)(x+1& )\equiv 0(modp)\\ p|(x-1& )(x+1)\end{array}$$
$\because$质数$p$只有$1$和他本身两个因数
$\therefore$所以必定有$p|x-1$或$p|x+1$
即$x-1\equiv 0(modp)$或$x+1\equiv 0(modp)$
解得$x\equiv 1(modp)$或$x\equiv p-1(modp)$。

所以算法流程如下：

1. 先用费马小定理进行判断
2. 如果$p-1$是偶数，则有$x^{p-1}\equiv 1(modp)\Rightarrow (x^{\frac{p-1}{2}}{)}^2\equiv 1(modp)$，所以可以不断循环用二次探测判定。

由于Miller_Rabin算法再判断的时候容易出现失误，所以我们最好选择多个$x$进行测试判定。
我们在选择$x$的时候，最好选择质数。显然当我们选择的$x$质数越多，第一判断失误的数$T$也会越大。
有人发现，选择 { 2 , 3 , 7 , 61 , 24251 } \{2,3,7,61,24251\} {2,3,7,61,24251}这$5$个数作为底数的时候，在$10^{16}$的范围内，只有判断$46856248255981$时会出现失误，这个数会被判定为素数，而他是个合数$46856248255981=4840261\ast 9680521$。

也有一篇文章《64位以内Rabin-Miller 强伪素数测试和Pollard rho 因数分解算法的实现》提到，
选取前$k$个素数为基进行测试判定，并用${\psi }_k$表示以前$k$个素数为基进行素数判定第一个出现失误的合数，Riesel在1994年给出下表：
$$\begin{array}{rl}{\psi }_1& =2,047\\ {\psi }_2& =1,373,653\\ {\psi }_3& =25,326,001\\ {\psi }_4& =3,215,031,751\\ {\psi }_5& =2,152,302,898,747\\ {\psi }_6& =3,474,749,660,383\\ {\psi }_7& =341,550,071,728,321\\ {\psi }_8& =341,550,071,728,321\\ {\psi }_9& \le 41,234,316,135,705,689,041\\ {\psi }_{10}& \le 1,553,360,566,073,143,205,541,002,401\\ {\psi }_{11}& \le 56,897,193,526,942,024,370,326,972,321\end{array}$$

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<queue>
#define ll long long
#define ull unsigned long long
#define ld long double
#define inl inline
#define re register
#define MAXN 10100
using namespace std;
inl ll ksc(ull x,ull y, ll p){return (x*y-(ull)((ld)x/p*y)*p+p)%p;}
inl ll ksm(ll x,ll y,ll p)
{
	ll res=1;
	for(;y;y>>=1,x=ksc(x,x,p))if(y&1)res=ksc(res,x,p);
	return res;
}
inl bool mr(ll x,ll p)
{
	if(ksm(x,p-1,p)!=1)return 0;
	ll y=p-1,z;
	while(!(y&1))
	{
		y>>=1;z=ksm(x,y,p);
		if(z!=1&&z!=p-1)return 0;
		if(z==p-1)return 1;
	}
	return 1;
}
int te_pri[20]={0,2,3,5,7,43,61,24251},te_num=7;
inl bool isprime(ll x)
{
	if(x<2)return 0;
	for(int i=1;i<=te_num;i++)
		if(x==te_pri[i])return 1;
	for(int i=1;i<=te_num;i++)
		if(!(x%te_pri[i])||!mr(te_pri[i],x))return 0;
	return 1;
}
int main()
{
	return 0;
}