1. 程式人生 > >Miller Rabin素數檢測與Pollard Rho演算法

Miller Rabin素數檢測與Pollard Rho演算法

一些前置知識可以看一下我的[聯賽前數學知識](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