发现自己竟然菜到不太会龟速乘，所以把 \(Miller-Rabin\) 算法所需要用到的算法全学了一遍……

龟速乘

龟速乘是一种 \(O(\log n)\) 的乘法计算方法。

考虑有时普通乘法取模会爆 \(long\ long\)，因此我们考虑用类似快速幂的方式进行乘法运算。

int mul(int x,int y,int c){
	x%=c,y%=c;
	int re=0;
	while(y){
		if(y&1) re+=x,re-=(re>=c)?c:0;
		x+=x,x-=(x>=c)?c:0,y>>=1;
	}return re;
}

快速乘

快速乘是一种 \(O(1)\) 的乘法计算方法。

发现 \(long\ double\) 可以存储 \(10^{300}\) 的数，所以考虑使用 \(long\ double\) 进行乘法运算。我还没有遇到过出锅的情况。

int mul(int a,int b,int p){
	return (a*b-(int)((long double)a/p*b)*p+p)%p;
}

\(Miller-Rabin\) 算法

哈，这玩意错的跟对的一样。

和快速幂求逆元一样，我们考虑费马小定理。

我们可以得到：

设 \(f(x)=a^{x-1}\bmod x\)，则当 \(f(p)\ne 1\) 时，\(p\) 一定是合数；当 \(p\) 是质数时，\(f(p)=1\)。

但是正确率不是人类所能接受的，所以考虑优化。

考虑下面这个性质：

若 \(b^2\bmod p=1\)，\(p\) 为质数，则 \(b-1|p\) 或 \(b+1|p\)。

这样我们就可以进行二次探测，正确率大大提高。

我们设用于测试的集合为 \(A\)，则 \(A=\{2,3,5,7,11,13,17,19,23\}\) 时，对于 \(x\le 10^{18}\) 的情况，都可以正确判断 \(x\) 的素性。

int mr[9]={2,3,5,7,11,13,17,19,23};
int mul(int a,int b,int p){
	return (a*b-(int)((long double)a/p*b)*p+p)%p;
}int qpow(int x,int y,int c){
	int re=1;
	while(y){
		if(y&1) re=mul(re,x,c);
		x=mul(x,x,c),y>>=1;
	}return re;
}int check(int x,int md){
	int c=md-1,mid=qpow(x,c,md);
	if(mid!=1) return 0;
	while(!(c&1)&&mid==1)
		c>>=1,mid=qpow(x,c,md);
	return (mid==1||mid==md-1);
}int miller_rabin(int x){
	if(x<2) return 0;
	if(x<=23){
		for(int i=0;i<9;i++)
			if(mr[i]==x) return 1;
		return 0;
	}for(int i=0;i<9;i++)
		if(!check(mr[i],x)) return 0;
	return 1;
}