64位以內(nèi)Rabin-Miller強(qiáng)偽素?cái)?shù)測試

和Pollard因數(shù)分解算法的實(shí)現(xiàn)在求解POJ1811題PrimeTest中應(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ù)分解。這兩種算法相對于試除法都顯得比較復(fù)雜。本文試圖對這兩者進(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ù)，則對任意a有。特別的，如果p不能整除a，則還有。利用Fermat小定理可以得到一個(gè)測試合數(shù)的有力算法：對，選擇，計(jì)算，若結(jié)果不等于1則n是合數(shù)。若結(jié)果等于1則n可能是素?cái)?shù)，并被稱為一個(gè)以a為基的弱可能素?cái)?shù)（weakprobableprimebasea，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ù)對于任意與其互素的整數(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對模n的平方根只可能是1和。因此對模n的平方根也只可能是1和。如果本身還是一個(gè)偶數(shù)，我們可以再取一次平方根……將這些寫成一個(gè)算法：（Rabin-Miller強(qiáng)偽素?cái)?shù)測試）記，其中d是奇數(shù)而s非負(fù)。如果，或者對某個(gè)有，則認(rèn)為n通過測試，并稱之為一個(gè)以a為基的強(qiáng)可能素?cái)?shù)（strongprobableprimebasea，a-SPRP）。若n是合數(shù)，則又稱之為一個(gè)以a為基的強(qiáng)偽素?cái)?shù)（strongpseudoprime）。這個(gè)測試同時(shí)被冠以Miller的名字是因?yàn)镸iller提出并證明了如下測試：如果擴(kuò)展黎曼猜想（extendedRiemannhypothesis）成立，那么對于所有滿足的基a，n都是a-SPRP，則n是素?cái)?shù)。盡管Monier和Rabin在1980年證明了這個(gè)測試的錯(cuò)誤概率（即合數(shù)通過測試的概率）不超過，單個(gè)測試相對來說還是相當(dāng)弱的（Pomerance、Selfridge和Wagstaff,Jr.證明了對任意都存在無窮多個(gè)a-SPRP）。但由于不存在“強(qiáng)Carmichael數(shù)”（任何合數(shù)n都存在一個(gè)基a試之不是a-SPRP），我們可以組合多個(gè)測試來產(chǎn)生有力的測試，以至于對足夠小的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ù)為基，則對小于的所有大于1的整數(shù)測試都是正確的；對大于或等于并小于的整數(shù)測試錯(cuò)誤的概率不超過。Rabin-Miller強(qiáng)偽素?cái)?shù)測試本身的形式稍有一些復(fù)雜，在實(shí)現(xiàn)時(shí)可以下面的簡單形式代替：對，如果則認(rèn)為n通過測試。代替的理由可簡單證明如下：仍然記，其中d是奇數(shù)而s非負(fù)。若n是素?cái)?shù)，由可以推出或。若為前者，顯然取即可使n通過測試。若為后者，則繼續(xù)取平方根，直到對某個(gè)有，或。無論還是，n都通過測試。functionpowermod(a,s,n){ p:=1 b:=a whiles>0 {functionpowermod(a,s,n){ p:=1 b:=a whiles>0 { if(s&1)==1thenp:=p*b%n b:=b*b%n s:=s>>1 } returnp}functionmul64to128(a,b){functionmul64to128(a,b){ {ah,al}:={a>>32,a&0xffffffff} {bh,bl}:={b>>32,b&0xffffffff} rl:=al*bl c:=al*bh rh:=c>>32 c:=c<<32 rl:=rl+c ifrl<cthenrh++ //處理進(jìn)位 c:=ah*bl rh:=rh+(c>>32) c:=c<<32 rl:=rl+c ifrl<cthenrh++ //處理進(jìn)位 rh:=rh+ah*bh return{rh,rl}}乘法之后的取模運(yùn)算可以用浮點(diǎn)運(yùn)算快速完成。具體做法是乘積的高64位和低64位分別先對除數(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)算。偽代碼如下：functionmod64(rh,rl,n){functionmod64(rh,rl,n){ rh:=rh%n rl:=rl%n r:=fmodl(rh*2**64,n) r:=r+rl ifr<nthenr:=r-n //處理進(jìn)位 r:=fmodl(r,n) returnr}二、Pollard因數(shù)分解算法Pollard因數(shù)分解算法又稱為PollardMonteCarlo因數(shù)分解算法。它的核心思想是：選取一個(gè)隨機(jī)數(shù)a。再選一個(gè)隨機(jī)數(shù)b。檢查是否大于1。若大于1，就是n的一個(gè)因子。若不大于1，再選取隨機(jī)數(shù)c，檢查和。如此繼續(xù)，直到找到n的一個(gè)非平凡因子。它的主要實(shí)現(xiàn)方法是從某個(gè)初值出發(fā)，通過一個(gè)適當(dāng)?shù)亩囗?xiàng)式進(jìn)行迭代，產(chǎn)生一個(gè)偽隨機(jī)迭代序列直到其對模n出現(xiàn)循環(huán)。多項(xiàng)式的選擇直接影響著迭代序列的長度。經(jīng)典的選擇是選取，其中。不選擇0和的原因是避免當(dāng)序列中某一項(xiàng)時(shí)后續(xù)各項(xiàng)全部為1的情況。迭代序列出現(xiàn)循環(huán)的期望時(shí)間和期望長度都與成正比。設(shè)n有一個(gè)非平凡因子p。由前面可知，迭代序列關(guān)于模p出現(xiàn)循環(huán)的期望時(shí)間和期望長度與成正比。當(dāng)循環(huán)出現(xiàn)時(shí)，設(shè)，記，則d是p的倍數(shù)。當(dāng)時(shí)，d就是n的一個(gè)非平凡因子。functionpollard_floyd(n){ x:=x0 y:=x0 d:=1 whiled==1 { x:=f(x) y:=f(f(y)) d:=gcd(x–y,n) if1<dANDd<nthenreturnd ifd==nthenreturnFAILURE }}但是在分解成功之前，p是未知的，也就無從直接判斷循環(huán)是否已經(jīng)出現(xiàn)。仍然設(shè)迭代序列關(guān)于模p出現(xiàn)循環(huán)，并設(shè)。由取模運(yùn)算的性質(zhì)可知，即。故對任意，總有。記循環(huán)部分長度為t，若functionpollard_floyd(n){ x:=x0 y:=x0 d:=1 whiled==1 { x:=f(x) y:=f(f(y)) d:=gcd(x–y,n) if1<dANDd<nthenreturnd ifd==nthenreturnFAILURE }}functionpollard_brent(n){ x:=x0 y:=x0 k:=0 i:=1 d:=1 whiled==1 { k:=k+1 functionpollard_brent(n){ x:=x0 y:=x0 k:=0 i:=1 d:=1 whiled==1 { k:=k+1 x:=f(x) d:=gcd(x–y,n) if1<dANDd<nthenreturnd ifd==nthenreturnFAILURE ifk==ithen { y:=x i:=i<<1 } }}由于循環(huán)出現(xiàn)的時(shí)空期望都是，Pollard因數(shù)分解算法的總體時(shí)間復(fù)雜度也就是。在最壞情況下，其時(shí)間復(fù)雜度可能接近，但在一般條件下，時(shí)間復(fù)雜度都可以認(rèn)為是。參考資料：ChrisCaldwell“Fermat,probable-primalityandpseudoprimes.”ChrisCaldwell“Strongprobable-primalityandapracticaltest.”EricW.Weisstein“Brent’sFactorizati