POJ 1811 Prime Test (Rabin-Miller强伪素数测试 和Pollard-rho 因数分解)
題目鏈接
Description
Given a big integer number, you are required to find out whether it's a prime number.
Input
The first line contains the number of test cases T (1 <= T <= 20 ), then the following T lines each contains an integer number N (2 <= N < 254).
Output
For each test case, if N is a prime number, output a line containing the word "Prime", otherwise, output a line containing the smallest prime factor of N.
Sample Input
2
5
10
Sample Output
Prime
2
分析:
給定一個(gè)小于2^54的整數(shù),判斷該數(shù)是不是素?cái)?shù),如果是素?cái)?shù)的話,輸出“Prime”,否則輸出該數(shù)的最小的素因子。熟悉的有關(guān)素?cái)?shù)的算法在這道題看來都還是太弱了,所以得額外的找方法來解決。
應(yīng)用到的兩個(gè)重要算法是Rabin-Miller強(qiáng)偽素?cái)?shù)測試和Pollard ?因數(shù)分解算法。前者可以在 的時(shí)間內(nèi)以很高的成功概率判斷一個(gè)整數(shù)是否是素?cái)?shù)。后者可以在最優(yōu) 的時(shí)間內(nèi)完成合數(shù)的因數(shù)分解。這兩種算法相對(duì)于試除法都顯得比較復(fù)雜。本文試圖對(duì)這兩者進(jìn)行簡單的闡述,說明它們在32位計(jì)算機(jī)上限制在64位以內(nèi)的條件下的實(shí)現(xiàn)中的細(xì)節(jié)。下文提到的所有字母均表示整數(shù)。
一、Rabin-Miller強(qiáng)偽素?cái)?shù)測試
Rabin-Miller強(qiáng)偽素?cái)?shù)測試的基本思想來源于如下的Fermat小定理:
如果p是一個(gè)素?cái)?shù),則對(duì)任意a有 (a^p)%p=a。特別的,如果p不能整除a,則還有(a^(p-1))%p=1 。
利用Fermat小定理可以得到一個(gè)測試合數(shù)的有力算法:對(duì)n>1,選擇a>1,計(jì)算(a^(n-1))%n,若結(jié)果不等于1則n是合數(shù)。若結(jié)果等于1則n可能是素?cái)?shù),并被稱為一個(gè)以a為基的弱可能素?cái)?shù)(weak probable prime base a,a-PRP);若n是合數(shù),則又被稱為一個(gè)以a為基的偽素?cái)?shù)(pseudoprime)。
這個(gè)算法的成功率是相當(dāng)高的。在小于25,000,000,000的1,091,987,405個(gè)素?cái)?shù)中,一共只用21,853個(gè)以2為基的偽素?cái)?shù)。但不幸的是,Alford、Granville和Pomerance在1994年證明了存在無窮多個(gè)被稱為Carmichael數(shù)的整數(shù)對(duì)于任意與其互素的整數(shù)a算法的計(jì)算結(jié)果都是1。最小的五個(gè)Carmichael數(shù)是561、1,105、1,729、2,465和2,801。
考慮素?cái)?shù)的這樣一個(gè)性質(zhì):若n是素?cái)?shù),則1對(duì)模n的平方根只可能是1和-1 。因此a^(n-1) 對(duì)模n的平方根 a^((n-1)/2)也只可能是1和-1 。如果(n-1)/2本身還是一個(gè)偶數(shù),我們可以再取一次平方根……將這些寫成一個(gè)算法:
(Rabin-Miller強(qiáng)偽素?cái)?shù)測試)記n-1=(2^s)d,其中d是奇數(shù)而s非負(fù)。如果(a^d)%n=1 ,或者對(duì)某個(gè) 0<=r<s有(a^((2^r)d))%n=-1 ,則認(rèn)為n通過測試,并稱之為一個(gè)以a為基的強(qiáng)可能素?cái)?shù)(strong probable prime base a,a-SPRP)。若n是合數(shù),則又稱之為一個(gè)以a為基的強(qiáng)偽素?cái)?shù)(strong pseudoprime)。
這個(gè)測試同時(shí)被冠以Miller的名字是因?yàn)镸iller提出并證明了如下測試:如果擴(kuò)展黎曼猜想(extended Riemann hypothesis)成立,那么對(duì)于所有滿足1<a<2(log n)^2 的基a,n都是a-SPRP,則n是素?cái)?shù)。
盡管Monier和Rabin在1980年證明了這個(gè)測試的錯(cuò)誤概率(即合數(shù)通過測試的概率)不超過 1/4,單個(gè)測試相對(duì)來說還是相當(dāng)弱的(Pomerance、Selfridge和Wagstaff, Jr.證明了對(duì)任意a>1都存在無窮多個(gè)a-SPRP)。但由于不存在“強(qiáng)Carmichael數(shù)”(任何合數(shù)n都存在一個(gè)基a試之不是a-SPRP),我們可以組合多個(gè)測試來產(chǎn)生有力的測試,以至于對(duì)足夠小的n可以用來證明其是否素?cái)?shù)。
取前k個(gè)素?cái)?shù)為基,并用 來表示以前k個(gè)素?cái)?shù)為基的強(qiáng)偽素?cái)?shù),Riesel在1994年給出下表:
考慮到64位二進(jìn)制數(shù)能表示的范圍,只需取前9個(gè)素?cái)?shù)為基,則對(duì)小于φ8 的所有大于1的整數(shù)測試都是正確的;對(duì)大于或等于 φ8并小于2^64 的整數(shù)測試錯(cuò)誤的概率不超過1/(2^18) 。
Rabin-Miller強(qiáng)偽素?cái)?shù)測試本身的形式稍有一些復(fù)雜,在實(shí)現(xiàn)時(shí)可以下面的簡單形式代替:
對(duì) n>1,如果(a^(n-1))%n=1則認(rèn)為n通過測試。
代替的理由可簡單證明如下:
仍然記n-1=(2^s)d ,其中d是奇數(shù)而s非負(fù)。若n是素?cái)?shù),由(a^(n-1))%n=1可以推出(a^(2^(s-1)d))%n=1=(a^((n-1)/2))%n或(a^(2^(s-1)d))%n=-1=(a^((n-1)/2))%n 。若為前者,顯然取r=s-1 即可使n通過測試。若為后者,則繼續(xù)取平方根,直到對(duì)某個(gè) 1<=r<s有 (a^((2^r)d))%n=-1,或(a^(2d))%n=1 。無論 (a^d)%n=1還是 a^d)%n=-1,n都通過測試。
Rabin-Miller強(qiáng)偽素?cái)?shù)測試的核心是冪取模(即計(jì)算 )(a^s)%n。計(jì)算冪取模有以下的 O(log n)算法(以Sprache偽代碼語言描述):
這個(gè)算法在32位計(jì)算機(jī)上實(shí)現(xiàn)的難點(diǎn)在于指令集和絕大部分編程語言的編譯器都只提供了32位相乘結(jié)果為64位的整數(shù)乘法,浮點(diǎn)運(yùn)算由于精度的問題不能應(yīng)用于這里的乘法。唯一解決辦法是模仿一些編譯器內(nèi)建的64位整數(shù)乘法來實(shí)現(xiàn)兩個(gè)無符號(hào)64位相乘結(jié)果為128位的乘法。這個(gè)乘法可以將兩個(gè)乘數(shù)分別分割成兩個(gè)32位數(shù)來實(shí)現(xiàn)。為方便乘法之后的取模運(yùn)算,運(yùn)算結(jié)果應(yīng)當(dāng)用連續(xù)的128個(gè)二進(jìn)制位來表示。以下是其偽代碼:
乘法之后的取模運(yùn)算可以用浮點(diǎn)運(yùn)算快速完成。具體做法是乘積的高64位和低64位分別先對(duì)除數(shù)取模,然后再利用浮點(diǎn)單元合并取模。這里的浮點(diǎn)運(yùn)算要求浮點(diǎn)單元以最高精度運(yùn)算,計(jì)算前應(yīng)先將浮點(diǎn)單元控制字中的精度控制位設(shè)置為64位精度。為保證精度,應(yīng)當(dāng)用80位浮點(diǎn)數(shù)實(shí)現(xiàn)此運(yùn)算。偽代碼如下:
至此,Rabin-Miller強(qiáng)偽素?cái)?shù)測試的核心已經(jīng)全部實(shí)現(xiàn)。
二、Pollard-rho 因數(shù)分解算法
Pollard-rho因數(shù)分解算法又稱為Pollard Monte Carlo因數(shù)分解算法。它的核心思想是:選取一個(gè)隨機(jī)數(shù)a。再選一個(gè)隨機(jī)數(shù)b。檢查 gcd(a-b,n)是否大于1。若大于1, gcd(a-b,n)就是n的一個(gè)因子。若不大于1,再選取隨機(jī)數(shù)c,檢查 gcd(c-b,n)和 gcd(c-a,n)。如此繼續(xù),直到找到n的一個(gè)非平凡因子。
代碼:
#include <cstdio> #include <cstring> #include <cstdlib> #include <time.h> #include <iostream> #include <algorithm> using namespace std; #define ll __int64//**************************************************************** // Miller_Rabin 算法進(jìn)行素?cái)?shù)測試 //速度快,而且可以判斷 <2^63的數(shù) //**************************************************************** const int S=20;//隨機(jī)算法判定次數(shù),S越大,判錯(cuò)概率越小//計(jì)算 (a*b)%c. a,b都是long long的數(shù),直接相乘可能溢出的 /* ll Mult_mod (ll a,ll b,ll c) // a,b,c <2^63 {a%=c;b%=c;ll ret=0;while (b){if (b&1) {ret+=a;ret%=c;}a<<=1;if (a>=c)a%=c;b>>=1;}return ret; }*/ll Mult_mod (ll a,ll b,ll c) //減法實(shí)現(xiàn)比取模速度快 {//返回(a*b) mod c,a,b,c<2^63a%=c;b%=c;ll ret=0;while (b){if (b&1){ret+=a;if (ret>=c) ret-=c;}a<<=1;if (a>=c) a-=c;b>>=1;}return ret; }//計(jì)算 x^n %c ll Pow_mod (ll x,ll n,ll mod) //x^n%c {if (n==1) return x%mod;x%=mod;ll tmp=x;ll ret=1;while (n){if (n&1) ret=Mult_mod(ret,tmp,mod);//(a*b)%ctmp=Mult_mod(tmp,tmp,mod);n>>=1;}return ret; }//以a為基,n-1=x*2^t a^(n-1)=1(mod n) 驗(yàn)證n是不是合數(shù) //一定是合數(shù)返回true,不一定返回false bool Check (ll a,ll n,ll x,ll t) {ll ret=Pow_mod(a,x,n);//(a^x)%nll last=ret;for (int i=1; i<=t; i++){ret=Mult_mod(ret,ret,n);if(ret==1&&last!=1&&last!=n-1) return true; //合數(shù)last=ret;}if (ret!=1) return true;return false; }// Miller_Rabin()算法素?cái)?shù)判定 //是素?cái)?shù)返回true.(可能是偽素?cái)?shù),但概率極小) //合數(shù)返回false;bool Miller_Rabin (ll n) {if (n<2) return false;if (n==2) return true;if ((n&1)==0) return false;//偶數(shù)ll x=n-1;ll t=0;while ((x&1)==0)//不斷的對(duì)于x進(jìn)行右移操作{x>>=1;t++;}for (int i=0; i<S; i++){ll a=rand()%(n-1)+1; //rand()需要stdlib.h頭文件if (Check(a,n,x,t))return false;//合數(shù)}return true; }//************************************************ //pollard_rho 算法進(jìn)行質(zhì)因數(shù)分解 //************************************************ll factor[100];//質(zhì)因數(shù)分解結(jié)果(剛返回時(shí)是無序的) int tol;//質(zhì)因數(shù)的個(gè)數(shù)。數(shù)組下標(biāo)從0開始ll Gcd (ll a,ll b) {if (a==0) return 1; if (a<0) return Gcd(-a,b);while (b){ll t=a%b;a=b;b=t;}return a; }ll Pollard_rho (ll x,ll c) {ll i=1,k=2;ll x0=rand()%x;ll y=x0;while (true){i++;x0=(Mult_mod(x0,x0,x)+c)%x;ll d=Gcd(y-x0,x);if (d!=1 && d!=x) return d;if (y==x0) return x;if (i==k){y=x0;k+=k;}} } //對(duì)n進(jìn)行素因子分解 void Findfac (ll n) {if (Miller_Rabin(n)) //素?cái)?shù){factor[tol++]=n;return;}ll p=n;while (p>=n) p=Pollard_rho(p,rand()%(n-1)+1);Findfac(p);Findfac(n/p); }int main () // Poj 1811 交G++ 比c++ 快很多 {// srand(time(NULL));//需要time.h頭文件 //POJ上G++要去掉這句話int T;scanf("%d",&T);while (T--){ll n;scanf("%I64d",&n);if (Miller_Rabin(n)){printf("Prime\n");continue;}tol=0;Findfac(n);ll ans=factor[0];for (int i=1; i<tol; i++)if (factor[i]<ans)ans=factor[i];printf("%I64d\n",ans);}return 0; }轉(zhuǎn)載于:https://www.cnblogs.com/cmmdc/p/8779801.html
總結(jié)
以上是生活随笔為你收集整理的POJ 1811 Prime Test (Rabin-Miller强伪素数测试 和Pollard-rho 因数分解)的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 建行手机银行怎么开通?如何查询手机银行开
- 下一篇: bzoj 4898: [Apio2017