64位以內Rabin-Miller強偽素數(shù)測試和Pollard因數(shù)分解算法的實現(xiàn)

64位以內Rabin-Miller強偽素數(shù)測試

和Pollard因數(shù)分解算法的實現(xiàn)在求解POJ1811題PrimeTest中應用到的兩個重要算法是Rabin-Miller強偽素數(shù)測試和Pollard因數(shù)分解算法。前者可以在的時間內以很高的成功概率判斷一個整數(shù)是否是素數(shù)。后者可以在最優(yōu)的時間內完成合數(shù)的因數(shù)分解。這兩種算法相對于試除法都顯得比較復雜。本文試圖對這兩者進行簡單的闡述，說明它們在32位計算機上限制在64位以內的條件下的實現(xiàn)中的細節(jié)。下文提到的所有字母均表示整數(shù)。一、Rabin-Miller強偽素數(shù)測試Rabin-Miller強偽素數(shù)測試的基本思想來源于如下的Fermat小定理：如果p是一個素數(shù)，則對任意a有。特別的，如果p不能整除a，則還有。利用Fermat小定理可以得到一個測試合數(shù)的有力算法：對，選擇，計算，若結果不等于1則n是合數(shù)。若結果等于1則n可能是素數(shù)，并被稱為一個以a為基的弱可能素數(shù)（weakprobableprimebasea，a-PRP）;若n是合數(shù)，則又被稱為一個以a為基的偽素數(shù)（pseudoprime）。這個算法的成功率是相當高的。在小于25,000,000,000的1,091,987,405個素數(shù)中，一共只用21,853個以2為基的偽素數(shù)。但不幸的是，Alford、Granville和Pomerance在1994年證明了存在無窮多個被稱為Carmichael數(shù)的整數(shù)對于任意與其互素的整數(shù)a算法的計算結果都是1。最小的五個Carmichael數(shù)是561、1,105、1,729、2,465和2,801。考慮素數(shù)的這樣一個性質：若n是素數(shù)，則1對模n的平方根只可能是1和。因此對模n的平方根也只可能是1和。如果本身還是一個偶數(shù)，我們可以再取一次平方根……將這些寫成一個算法：（Rabin-Miller強偽素數(shù)測試）記，其中d是奇數(shù)而s非負。如果，或者對某個有，則認為n通過測試，并稱之為一個以a為基的強可能素數(shù)（strongprobableprimebasea，a-SPRP）。若n是合數(shù)，則又稱之為一個以a為基的強偽素數(shù)（strongpseudoprime）。這個測試同時被冠以Miller的名字是因為Miller提出并證明了如下測試：如果擴展黎曼猜想（extendedRiemannhypothesis）成立，那么對于所有滿足的基a，n都是a-SPRP，則n是素數(shù)。盡管Monier和Rabin在1980年證明了這個測試的錯誤概率（即合數(shù)通過測試的概率）不超過，單個測試相對來說還是相當弱的（Pomerance、Selfridge和Wagstaff,Jr.證明了對任意都存在無窮多個a-SPRP）。但由于不存在“強Carmichael數(shù)”（任何合數(shù)n都存在一個基a試之不是a-SPRP），我們可以組合多個測試來產生有力的測試，以至于對足夠小的n可以用來證明其是否素數(shù)。取前k個素數(shù)為基，并用來表示以前k個素數(shù)為基的強偽素數(shù)，Riesel在1994年給出下表：考慮到64位二進制數(shù)能表示的范圍，只需取前9個素數(shù)為基，則對小于的所有大于1的整數(shù)測試都是正確的；對大于或等于并小于的整數(shù)測試錯誤的概率不超過。Rabin-Miller強偽素數(shù)測試本身的形式稍有一些復雜，在實現(xiàn)時可以下面的簡單形式代替：對，如果則認為n通過測試。代替的理由可簡單證明如下：仍然記，其中d是奇數(shù)而s非負。若n是素數(shù)，由可以推出或。若為前者，顯然取即可使n通過測試。若為后者，則繼續(xù)取平方根，直到對某個有，或。無論還是，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++ //處理進位 c:=ah*bl rh:=rh+(c>>32) c:=c<<32 rl:=rl+c ifrl<cthenrh++ //處理進位 rh:=rh+ah*bh return{rh,rl}}乘法之后的取模運算可以用浮點運算快速完成。具體做法是乘積的高64位和低64位分別先對除數(shù)取模，然后再利用浮點單元合并取模。這里的浮點運算要求浮點單元以最高精度運算，計算前應先將浮點單元控制字中的精度控制位設置為64位精度。為保證精度，應當用80位浮點數(shù)實現(xià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 //處理進位 r:=fmodl(r,n) returnr}二、Pollard因數(shù)分解算法Pollard因數(shù)分解算法又稱為PollardMonteCarlo因數(shù)分解算法。它的核心思想是：選取一個隨機數(shù)a。再選一個隨機數(shù)b。檢查是否大于1。若大于1，就是n的一個因子。若不大于1，再選取隨機數(shù)c，檢查和。如此繼續(xù)，直到找到n的一個非平凡因子。它的主要實現(xiàn)方法是從某個初值出發(fā)，通過一個適當?shù)亩囗検竭M行迭代，產生一個偽隨機迭代序列直到其對模n出現(xiàn)循環(huán)。多項式的選擇直接影響著迭代序列的長度。經(jīng)典的選擇是選取，其中。不選擇0和的原因是避免當序列中某一項時后續(xù)各項全部為1的情況。迭代序列出現(xiàn)循環(huán)的期望時間和期望長度都與成正比。設n有一個非平凡因子p。由前面可知，迭代序列關于模p出現(xiàn)循環(huán)的期望時間和期望長度與成正比。當循環(huán)出現(xiàn)時，設，記，則d是p的倍數(shù)。當時，d就是n的一個非平凡因子。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)。仍然設迭代序列關于模p出現(xiàn)循環(huán)，并設。由取模運算的性質可知，即。故對任意，總有。記循環(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)的時空期望都是，Pollard因數(shù)分解算法的總體時間復雜度也就是。在最壞情況下，其時間復雜度可能接近，但在一般條件下，時間復雜度都可以認為是。參考資料：ChrisCaldwell“Fermat,probable-primalityandpseudoprimes.”ChrisCaldwell“Strongprobable-primalityandapracticaltest.”EricW.Weisstein“Brent’sFactorizati