素数判定 学习笔记

费马素性测试

费马小定理：

$$a^{p-1}\equiv 1(modp)[p\in \mathbb{Prime},gcd(a,p)=1]$$

因为费马小定理的逆命题对于绝大部分的数是正确的，所以我们可以用它来测试素性。
我们可以随机若干个 $a$，对他进行检验：

• 如果 $a^{p-1}\equiv 1(modp)$，那么继续检测。
• 否则其为合数。

相关题目：ARC017A。

代码

int p;
int qpo(int a,int b,const int mod){
	a%=mod;
	int res=1;
	while(b){
		if(b&1) res=(res*a)%mod;
		a=(a*a)%mod;
		b>>=1;
	}
	return res;
}
mt19937 rd(time(NULL));
signed main() {
	in1(p);
	For(i,1,1000000) {
		int a=abs((int)rd());
		if(a%p==0) continue;
		if(qpo(a,p-1,p)!=1) {
			cout<<"NO\n";
			return 0;
		}
	}
	cout<<"YES\n";
	return 0;
}

但是值得注意的是，有这样一类数，如 $561,1105,1729,2465,2821,\cdots$，被称作 Carmichael 数，他们虽然满足对于任意一个 $a$ 都满足 $a^{p-1}\equiv 1(modp)[gcd(a,p)=1]$，但是他们并不是素数，比如 $561=3\times 11\times 17$，所以我们在随机 $a$ 时可以通过 if(a%p==0) continue; 的方式筛掉大部分的这些数。

证明

当 $gcd(a,p)=1$ 且 $p$ 为素数时，$a=kp$，
所以一定不会出现 $(kp{)}^{p-1}$。但是 $p$ 不为素数时，
$p\mid a$ 不一定代表 $gcd(a,b)=1$，
所以可能出现一个 $a$ 可以筛掉这个数。

Miller-Rabin 素性测试

参考文章www

https://oi-wiki.org/math/number-theory/prime/#millerrabin-素性测试

https://www.luogu.com.cn/article/jr27vtto

二次探测定理：$x^2\equiv 1(modp)$ 且 $p\ge 3$ 且 $p$ 为素数时，$x$ 的解为：$x\equiv 1(modp)$ 或 $x\equiv p-1(\mathrm{mod}p)$。

证明

$$\because x^2\equiv 1(modp)$$

$$\therefore x^2-1\equiv 0(modp)$$

$$\therefore (x-1)(x+1)\equiv 0(modp)$$

$$x-1\mid p\mathrm{或}x+1\mid p$$

$$x_1\equiv 1(modp),x_2\equiv p-1(modp)$$

设 $p-1=2^{t}u$，那么根据费马小定理，我们有 $a^{p-1}\equiv (a^u{)}^{2^t}(modp)$，用快速幂计算 $a^u$。然后重复执行 $t$ 次平方。因为每一次平方都可以用二次探测定理检验，遂最多计算 $t$ 次平方。

通过选取一些固定的数字，也可以实现一定范围内的正确性。
在考场上选取 $\{2,3,5,7,11,13,17,19,23,29,31,37\}$ 可以判断 $2^{78}$ 以内的数字。
如果你想知道前 $n$ 个素数可以判断多少数字，参见A014233。

相关题目：SPOJ 288 PON - Prime or Not。

代码

#define int __int128
int t,p;
int tst[15]={0,2,3,5,7,11,13,17,19,23,29,31,37};
int qpo(int a,int b,int mod) {
	int res=1;
	for(;b;b>>=1,a=(a*a)%mod) if(b&1) res=(res*a)%mod;
	return res;
}
void divide(int p,int &x,int &u) {
	while(p%(1<<(x+1))==0) x++;
	u=p/(1<<x);
}
bool check(int p) {
	if(p<=40) {
		For(i,1,12) if(p==tst[i]) return 1;
		return 0;
	}
	int t=0,u=0;
	divide(p-1,t,u);
	For(i,1,12) {
		int a=tst[i];
		a=qpo(a,u,p);
		if(a==1) continue;
		bool flg=0;
		For(j,1,t) {
			if(a==p-1) {
				flg=1;
				break;
			}
			a=(a*a)%p;
		}
		if(!flg) return 0;
	}
	return 1;
}