素数判定法（Miller-Rabin素数判定）

一、依据原理
1.费马小定理：当p为素数时，且a,p互质，则有a^(p-1) ≡ 1 mod p.
2.二次探测：如果p是一个素数，0<x<p，则方程x^2 ≡ 1 mod p的解为x=1或x=p-1.

如何运用定理来完成素数的判定？

首先，1 mod p的结果有两个：0/1. 当1 mod p =0, p = 1，p既不是素数也不是合数。所以p为素数时，[a^(p-1)] mod p =1
因此：对于一个给定的与p互质的整数a，[a^(p-1)] mod p ！= 1，说明p不是素数。

第二：对于给定的整数x∈(0,p)，若x^2 mod p =1, x只能等于1或者等于p-1，如果都不等于，说明p不是素数。

具体如何运用到编程实现？

1、 对于p-1，将其进行分解：
首先，由于p是素数，那么p一定是奇数（2被特判），所以p-1为偶数。
对p-1进行拆分，即将（p-1）不断除以2，拆分为u * 2^k，一直到u为奇数时停止。

然后在[1,p-1]范围内随机生成一个整数a, 计算a^u:
(也就是x=a^u，为了保证x∈(0,p)，x=a^u mod p【这样就可以保证小于p】)
然后不断对x求平方，并进行二次探测。
这个过程要一直重复k次，此时x = (a^u)2^k， 又u*2^k=p-1，也就是此时x=a^(p-1)，然后利用费马小定理判断此时x mod p是否等于1.

这个随机生成a的过程让它重复十次，降低错误概率。
**************************************************************
至于为什么是在[1,p-1]范围内随机生成一个整数a，我的想法是如果p是一个素数，也就是说除了p自己和1，没有任何因数，所以在1~p-1范围内随机选一个数，都应该和p互质。

二、编程实现

//快速乘算法
快速乘是利用了二进制和十进制之间的相互转化，以及乘法的因式分解实现
例如：12*11 mod 5
12 * 11 = 12 * 1011(2)
1011(2) = 2^3 * 1 + 2^2 * 0 + 2^1 * 1 + 2^0 * 1 = 8+2+1
12*(8+2+1) = 96+48+12
12*11 mod 5 = (96+48+12) mod 5 = （96 mod 5 + 48 mod 5 + 12 mod 5）mod 5

#include<iostream>
using namespace std;
typedef unsigned long long ULL;
const int maxn=1e5+7;

//定义快速乘算法
//a代表乘数A，b代表乘数B，mod代表模的数，ans是要返回的答案
ULL MMul(ULL a, ULL b, ULL mod, ULL ans = 0){
	while(b){
		if(b&1) ans=(ans+a)%mod;  
		a=(a+a)%mod;
		b>>=1;}
	return ans;
}

/*代码介绍
利用b>>=1，也就是右移操作，每次都将b除以2,利用b的值控制while循环，当b除以2==0的时候，也就是十进制转为二进制操作已经结束。
当b&1为True，也就是当前二进制是1而不是0，那么就需要将当前的a加入到ans中
(a+a) --> 相当于a*(2^1) a*(2^2) a*(2^3) ...
*/

//快速幂算法
具体可以参考https://blog.csdn.net/Gsiredo/article/details/129678182

//定义快速幂算法
ULL MExp(ULL a, ULL b, ULL mod, ULL ans=0){
	while(b){
		if(b&1) MMul(ans,a,mod);
		a=MMul(a,a,mod);
		b>>=1;{
	return ans;
}

#include<cstdlib>
//实现素数判定
//参数含义： n是需要判定的数；u是对n分解后奇数的部分，t是对n分解后2的指数部分；s是重复判定的次数
bool Miller_Rabin(ULL n, ULL u=0; int t=0; int s=10){
	//首先对2进行特判
	if(n==2) return true;
	//对1或者偶数进行特判
	if(n<2||!(n&1)) return false;
	//对n进行拆解，得到奇数u，以及2的指数t
	for(t=0,u=n-1;!(u&1);t++){
		u>>=1;}
	//一共进行s次判定，每次利用rand()函数随机生成[1,p-1]范围内的随机数a，所以要加入<cstdlib>头文件的引用
	while(s--){
		ULL a=rand()%(n-1)+1;
		ULL x=MExp(a,u,n);
		for(int i=0;i<t;i++){
			ULL y=MMul(x,x,n);
			if(y==1&&x!=1&&x!=n-1) return false;
			x=y;}
		if(x!=1) return false;}
	return true;}

//构造主函数入口
int main(){
	ULL t;
	while(cin>>t){
		printf("%s\n",Miller_Rabin(t)?"Yes":"No");}
	return 0;}