Miller Rabin素數檢測與Pollard Rho演算法
阿新 • • 發佈:2021-01-07
一些前置知識可以看一下我的[聯賽前數學知識](https://www.cnblogs.com/liuchanglc/p/13692477.html)
## 如何判斷一個數是否為質數
### 方法一:試除法
掃描$2\sim \sqrt{n}$之間的所有整數,依次檢查它們能否整除$n$,若都不能整除,則$n$是質數,否則$n$是合數。
#### 程式碼
``` cpp
bool is_prime(int n){
if(n<2) return 0;
int m=sqrt(n);
for(int i=2;i<=m;i++){
if(n%i==0) return 0;
}
return 1;
}
```
### 方法二、線性篩
用 $O(n)$ 的複雜度篩出所有的質數,然後 $O(1)$ 判斷
#### 程式碼
``` cpp
const int maxn=1e6+5;
int pri[maxn];
bool not_prime[maxn];
void xxs(int n){
not_prime[0]=not_prime[1]=1;
for(int i=2;i<=n;++i){
if(!not_prime[i]){
pri[++pri[0]]=i;
}
for (int j=1;j<=pri[0] && i*pri[j]<=n;j++){
not_prime[i*pri[j]]=1;
if(i%pri[j]== 0) break;
}
}
}
```
但是,我們來看一下這道題 [LOJ#143. 質數判定](https://loj.ac/p/143)
在 $5s$ 的時間內判斷 $10^5$ 個 $10^{18}$ 級別的數是否是質數
以上兩種方法顯然都不可行
所以我們要用到一種更高效的演算法 $Miller\ Rabin$
## Miller Rabin素數檢測
根據費馬小定理
若 $p$ 為素數,對於任意整數 $a$,都有
$a^{p-1} \equiv 1(mod\ p)$
它的逆定理在大多數情況下是成立的
所以,可以每一次隨機挑選一些數 $a$
判斷是否存在 $a^{p-1} \equiv 1(mod\ p)$
只要有一個 $a$ 不滿足條件,就將其標記為合數
否則標記為質數
### 程式碼
``` cpp
#include
#include
#include
#define rg register
typedef long long ll;
ll gsc(rg ll ds,rg ll zs,rg ll mod){
return ((unsigned long long)(ds*zs)-(unsigned long long)((long double)ds/mod*zs)*mod+mod)%mod;
}
ll ksm(rg ll ds,rg ll zs,rg ll mod){
rg ll nans=1;
while(zs){
if(zs&1LL) nans=gsc(nans,ds,mod);
ds=gsc(ds,ds,mod);
zs>>=1LL;
}
return nans;
}
int main(){
rg ll aa,bb;
rg bool jud=1;
while(scanf("%lld",&aa)!=EOF){
jud=1;
if(aa==1){
printf("N\n");
continue;
}
for (rg int i=1;i<=100;i++) {
bb=1LL*rand()*rand()%(aa-1)+1;
if(ksm(bb,aa-1,aa)!=1){
jud=0;
break;
}
}
if(jud) printf("Y\n");
else printf("N\n");
}
return 0;
}
```
但是也有例外,即存在一種極端反例卡邁克爾數(一種合數)
對於任何卡邁克爾數,費馬定理都成立
雖然這種數很少,但是還是有可能會被卡(只有$33$分)
所以我們需要藉助另一個定理來優化它
二次探測定理:若 $p$ 為質數且 $x∈(0,p)$,則方程 $x^2≡1(mod\ p)$ 的解為 $x_1=1,x_2=p-1$
對於任意一個奇素數 $p$ ,$p-1$ 一定可以寫成 $r2^t$ 的形式
因此我們可以把費馬小定理寫成 $a^{r2^t}\equiv 1(mod\ p)$ 的形式
一開始,我們先求出 $a^r mod\ p$ 的值
然後每一次給這個值平方,一共平方 $t$ 次
算一下每次得出來的結果是否滿足二次探測定理
如果不滿足,說明這個數不是質數
最後再看一下最終的值是否滿足費馬小定理即可
時間複雜度 $klog^2n$
其中 $k$ 為每次檢測的次數
事實證明,這個演算法出錯的概率非常小
為 $\frac{1}{4^k}$
### 程式碼
``` cpp
#include
#include
#include
#include
#include
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
typedef long long ll;
const int times=7;
ll gsc(rg ll ds,rg ll zs,rg ll mod){
return ((unsigned long long)(ds*zs)-(unsigned long long)((long double)ds/mod*zs)*mod+2*mod)%mod;
}
ll ksm(rg ll ds,rg ll zs,rg ll mod){
rg ll nans=1;
while(zs){
if(zs&1LL) nans=gsc(nans,ds,mod);
ds=gsc(ds,ds,mod);
zs> >=1LL;
}
return nans;
}
bool check(rg ll a,rg ll r,rg ll t,rg ll mod){
ll nans=ksm(a,r,mod),tmp=nans;
for(rg ll i=1;i<=t;i++){
nans=gsc(tmp,tmp,mod);
if(nans==1 && tmp!=1 && tmp!=mod-1) return 0;
tmp=nans;
}
if(tmp==1) return 1;
else return 0;
}
bool Miller_Robin(rg ll n){
if(n==2) return 1;
if(n<2 || (n&1LL)==0) return 0;
rg ll t=0,r=n-1;
while((r&1LL)==0){
r>>=1;
t++;
}
for(rg int i=1;i<=times;i++){
rg ll a=rand()%(n-1)+1;
if(!check(a,r,t,n)) return 0;
}
return 1;
}
ll x;
int main(){
srand(time(0));
while(scanf("%lld",&x)!=EOF){
if(Miller_Robin(x)) printf("Y\n");
else printf("N\n");
}
return 0;
}
```
## Pollard Rho演算法
[模板題](https://www.luogu.com.cn/problem/P4718)
這道題不僅要判斷是否是質數,還要求輸出最大質因子
判質數的話用 $Miller\ Rabin$ 判一下就好了
關鍵是怎麼找出最大質因子
其實還是隨機的思想
我們每一次隨機一個數 $n$,判斷 $n$ 是不是 $p$ 的因數
如果是,就分成兩個子問題 $n$ 和 $\frac{p}{n}$ 遞迴求解
如果當前的數為質數就更新答案
但是這樣隨機概率非常小
利用生日悖論,採用組合隨機取樣的方法,滿足答案的組合比單個個體要多一些.這樣可以提高正確率
具體方法是隨機兩個整數 $n$,$m$
判斷 $|n-m|$ 是不是 $p$ 的因數
看起來沒有什麼不同,實際上效率提高了不少
剩下的就在於怎麼構造隨機的 $n,m$ 了
構造的方法會影響到我們程式的效率
實踐證明,直接 $rand$ 效率極低
而構造一種形如 $f(x)=x^2+c$ 的數列效率很高
我們首先 $rand$ 一個 $x_1$
然後令 $f(x_2)=x_1^2+c$
$f(x_3)=x_2^2+c$
依次類推,每次把序列相鄰兩項作差判斷是否是 $p$ 的因數
但是因為我們是在模意義下進行運算,所以這個序列一定有迴圈節
所以當遇到迴圈節的時候我們就退出迴圈,重新構造一個數列
因為每次求 $gcd$ 的開銷太大
所以我們可以先把相鄰兩項的差連乘起來,這是不影響結果的
當累乘的次數等於我們設定的一個值時再進行求 $gcd$ 的運算
一般把這個值設為 $2^k$
### 程式碼
``` cpp
#include
#include
#include
#include
#include
#include
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
typedef long long ll;
const int times=10;
int nt;
ll x,ans;
ll gsc(rg ll ds,rg ll zs,rg ll mod){
return ((unsigned long long)(ds*zs)-(unsigned long long)((long double)ds/mod*zs)*mod+mod)%mod;
}
ll ksm(rg ll ds,rg ll zs,rg ll mod){
rg ll nans=1;
while(zs){
if(zs&1LL) nans=gsc(nans,ds,mod);
ds=gsc(ds,ds,mod);
zs> >=1LL;
}
return nans;
}
bool check(rg ll a,rg ll r,rg ll t,rg ll mod){
ll nans=ksm(a,r,mod),tmp=nans;
for(rg ll i=1;i<=t;i++){
nans=gsc(tmp,tmp,mod);
if(nans==1 && tmp!=1 && tmp!=mod-1) return 0;
tmp=nans;
}
if(tmp==1) return 1;
else return 0;
}
bool Miller_Robin(rg ll n){
if(n==2) return 1;
if(n<2 || (n&1LL)==0) return 0;
rg ll t=0,r=n-1;
while((r&1LL)==0){
r>>=1;
t++;
}
for(rg int i=1;i<=times;i++){
rg ll a=rand()%(n-1)+1;
if(!check(a,r,t,n)) return 0;
}
return 1;
}
ll gcd(rg ll aa,rg ll bb){
if(bb==0) return aa;
return gcd(bb,aa%bb);
}
ll rp(rg ll x,rg ll y){
rg int t=0,k=1;
rg ll v0=rand()%(x-1)+1,v=v0,d,s=1;
while(1){
v=(gsc(v,v,x)+y)%x;
s=gsc(s,std::abs(v-v0),x);
if(v==v0 || !s) return x;
if(++t==k){
d=gcd(s,x);
if(d!=1) return d;
s=1,v0=v,k<<=1;
}
}
}
void dfs(rg ll n){
if(n<=ans || n==1) return;
if(Miller_Robin(n)){
ans=n;
return;
}
rg ll now=n;
while(now==n) now=rp(n,rand()%n+1);
while(n%now==0) n/=now;
dfs(n),dfs(now);
}
ll solve(rg ll n){
ans=1;
dfs(n);
return ans;
}
int main(){
srand(time(0));
scanf("%d",&nt);
while(nt--){
scanf("%lld",&x);
if(Miller_Robin(x)) printf("Prime\n");
else printf("%lld\n",solve(x));
}
return