博客园同步

原题链接

简要题意：

给定 $T$ 个数 $n$，判素数。

$T \leq 10^5 , n \leq 10^{18}$.

可能你判一个都有困难是不是 $\cdots \cdots$.

二次探测定理

若 $x^2 \equiv 1 \pmod p , x < p$，则 $x= 1$ 或 $x = p-1$.

简单证明：

$x^2 \equiv 1 \pmod p$

$(x-1)(x+1) \equiv 0 \pmod p$

没了。

$\text{Miller - Rabin}$ 算法

如何快速测试一个数是否为素数？这是基于 二次探测定理 的。

大家一定知道 伪素数 吧！就是那些 费马小定理逆定理的反例，以其中最小的 $341$ 为例，先以 $2$ 为底：

$2^{340} \equiv 1 \pmod {341}$，第一次通过。

$2^{170} \equiv 1 \pmod {341}$，第二次通过。

$2^{85} \equiv 32 \pmod {341}$，说明 $341$ 不是素数。当然如果这一次通过则说明 $341$ 通过了 $2$ 的底数检测，因为 $85$ 是奇数。

当然同样的道理，我们可以以其它素数为底进行判断。这样的成功效率是多少呢？

假设你判断了 $k$ 个素数，那么错误的概率是 $4^{-k}$，实际上和 $0$ 也差不多。

只用 $2,7,61$ 进行测试尽可以保证 $4 759 123 141$ 之内正确。

如果你用 $2,3,5,7,11,13,17$ 进行测试，可以保证 $341 550 071 728 320$ 之内所有数正确。

如果选用 $2, 3, 7, 61,24251$ 作为底数， $10^{16}$ 之内就存在一个反例： $46856248255981$.（当然是 $10^{16}$ 之内唯一的反例，说明了正确的概率之高）

这题是 $10^{18}$，本人采用 $2,3,5,7,11,61,24251$ 可以通过。（当然你也可以用合数进行测试，但合数效率不高）另外本人还添加了 $10$ 组随机测试（即随机生成底数进行测试），保证了正确性。

时间复杂度： $\mathcal{O}(T \log^2 n)$. 错误概率： $4^{-k}$.

实际得分： $100pts$.

细节：

快速幂中的相乘会溢出，本人的编译器不知道怎么 typedef __int128 ll 会编译错误，所以用了 long double 进行乘法的转移。（不是龟速乘啊）

#pragma GCC optimize(2)
#include<bits/stdc++.h>
using namespace std;

typedef long long ll;

inline ll read(){char ch=getchar(); int f=1; while(ch<'0' || ch>'9') {if(ch=='-') f=-f; ch=getchar();}
	ll x=0; while(ch>='0' && ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar(); return x*f;}

inline void write(int x) {
	if(x<0) {putchar('-');write(-x);return;}
	if(x<10) {putchar(char(x%10+'0'));return;}
	write(x/10);putchar(char(x%10+'0'));
}
ll prime[8]={2,3,5,7,11,61,24251};
inline ll times(ll x,ll y,ll MOD) {
	x%=MOD; y%=MOD;
	ll t=(long double) x/MOD*y;
	ll p=x*y-t*MOD;
	return ((p%MOD)+MOD)%MOD;
} //计算 x*y

inline ll pw(ll x,ll y,ll MOD) {
	ll ans=1; while(y) {
		if(y&1) ans=times(ans,x,MOD);
		x=times(x,x,MOD); y>>=1;
	} return ans;
} //快速幂

inline bool check(ll x,ll y,ll MOD) {
	ll p=pw(x,y,MOD);
	if(p!=1 && p!=MOD-1) return 0; //测试失败
	if(p==MOD-1 || (p==1 && (y&1))) return 1; //测试成功
	return check(x,y>>1,MOD); //继续测试
}

inline bool Miller_Rabin(ll n) {
	if(n<=1) return 0;
	for(int i=0;i<7;i++) {
		if(n==prime[i]) return 1;
		if(n%prime[i]==0) return 0;
		if(!check(prime[i],n-1,n)) return 0; // 7 个素数轮流测试
	} for(int i=1;i<=10;i++) {
		int t=2+rand()%(n-2);
		if(!check(t,n-1,n)) return 0; // 10 个随机
	} return 1;
}

int main() {
	ll x; while(~scanf("%lld",&x)) puts(Miller_Rabin(x)?"Y":"N"); //套路
	return 0;
}