CSDN同步

原題連結

簡要題意：

給定 \(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) \pmod p \]

沒了。

\(\text{Miller - Rabin}\) 演算法

如何快速測試一個數是否為素數？這是基於 二次探測定理 的。

大家一定知道 偽素數

\(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}\)

只用 \(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;
}