素数判定 学习笔记

费马素性测试

费马小定理：

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

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

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

相关题目：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 (\bmod p)[\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(\bmod p)\) 且 \(p\ge 3\) 且 \(p\) 为素数时，\(x\) 的解为：\(x\equiv 1(\bmod p)\) 或 \(x\equiv p-1(\mod p)\)。

证明

\[\because x^2 \equiv 1(\bmod p) \]

\[\therefore x^2-1 \equiv 0(\bmod p) \]

\[\therefore (x-1)(x+1) \equiv 0 (\bmod p) \]

\[x-1 \mid p 或 x+1 \mid p \]

\[x_1\equiv 1(\bmod p),x_2\equiv p-1(\bmod p) \]

设 \(p-1=2^tu\)，那么根据费马小定理，我们有 \(a^{p-1}\equiv (a^u)^{2^t} (\bmod p)\)，用快速幂计算 \(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;
}