在求解POJ1811題Prime Test中應(yīng)用到的兩個重要算法是Rabin-Miller強偽素數(shù)測試和Pollard r因數(shù)分解算法。前者可以在的時間內(nèi)以很高的成功概率判斷一個整數(shù)是否是素數(shù)。后者可以在最優(yōu)的時間內(nèi)完成合數(shù)的因數(shù)分解。這兩種算法相對于試除法都顯得比較復(fù)雜。本文試圖對這兩者進行簡單的闡述，說明它們在32位計算機上限制在64位以內(nèi)的條件下的實現(xiàn)中的細節(jié)。下文提到的所有字母均表示整數(shù)。

?

一、Rabin-Miller強偽素數(shù)測試

Rabin-Miller強偽素數(shù)測試的基本思想來源于如下的Fermat小定理：

如果p是一個素數(shù)，則對任意a有。特別的，如果p不能整除a，則還有。

利用Fermat小定理可以得到一個測試合數(shù)的有力算法：對，選擇，計算，若結(jié)果不等于1則n是合數(shù)。若結(jié)果等于1則n可能是素數(shù)，并被稱為一個以a為基的弱可能素數(shù)（weak probable prime base a，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算法的計算結(jié)果都是1。最小的五個Carmichael數(shù)是561、1,105、1,729、2,465和2,801。

考慮素數(shù)的這樣一個性質(zhì)：若n是素數(shù)，則-1對模n的平方根只可能是1和。因此對模n的平方根也只可能是1和。如果本身還是一個偶數(shù)，我們可以再取一次平方根……將這些寫成一個算法：

（Rabin-Miller強偽素數(shù)測試）記，其中d是奇數(shù)而s非負。如果，或者對某個有，則認為n通過測試，并稱之為一個以a為基的強可能素數(shù)（strong probable prime base a，a-SPRP）。若n是合數(shù)，則又稱之為一個以a為基的強偽素數(shù)（strong pseudoprime）。

這個測試同時被冠以Miller的名字是因為Miller提出并證明了如下測試：如果擴展黎曼猜想（extended Riemann hypothesis）成立，那么對于所有滿足的基a，n都是a-SPRP，則n是素數(shù)。

盡管Monier和Rabin在1980年證明了這個測試的錯誤概率（即合數(shù)通過測試的概率）不超過，單個測試相對來說還是相當弱的（Pomerance、Selfridge和Wagstaff, Jr.證明了對任意都存在無窮多個a-SPRP）。但由于不存在“強Carmichael數(shù)”（任何合數(shù)n都存在一個基a試之不是a-SPRP），我們可以組合多個測試來產(chǎn)生有力的測試，以至于對足夠小的n可以用來證明其是否素數(shù)。

取前k個素數(shù)為基，并用來表示以前k個素數(shù)為基的強偽素數(shù)，Riesel在1994年給出下表：

考慮到64位二進制數(shù)能表示的范圍，只需取前9個素數(shù)為基，則對小于的所有大于1的整數(shù)測試都是正確的；對大于或等于并小于的整數(shù)測試錯誤的概率不超過。

Rabin-Miller強偽素數(shù)測試本身的形式稍有一些復(fù)雜，在實現(xiàn)時可以下面的簡單形式代替：

對，如果則認為n通過測試。

代替的理由可簡單證明如下：

仍然記，其中d是奇數(shù)而s非負。若n是素數(shù)，由可以推出或。若為前者，顯然取即可使n通過測試。若為后者，則繼續(xù)取平方根，直到對某個有，或。無論還是，n都通過測試。

function powermod(a, s, n) { ?? p := 1 ?? b := a ?? while s > 0 ?? { ????? if (s & 1) == 1 then p := p * b % n ????? b := b * b % n ????? s := s >> 1 ?? } ?? return p }

Rabin-Miller強偽素數(shù)測試的核心是冪取模（即計算）。計算冪取模有以下的算法（以Sprache偽代碼語言描述）：

function mul64to128(a, b) { ?? {ah, al} := {a >> 32, a & 0xffffffff} ?? {bh, bl} := {b >> 32, b & 0xffffffff}

這個算法在32位計算機上實現(xiàn)的難點在于指令集和絕大部分編程語言的編譯器都只提供了32位相乘結(jié)果為64位的整數(shù)乘法，浮點運算由于精度的問題不能應(yīng)用于這里的乘法。唯一解決辦法是模仿一些編譯器內(nèi)建的64位整數(shù)乘法來實現(xiàn)兩個無符號64位相乘結(jié)果為128位的乘法。這個乘法可以將兩個乘數(shù)分別分割成兩個32位數(shù)來實現(xiàn)。為方便乘法之后的取模運算，運算結(jié)果應(yīng)當用連續(xù)的128個二進制位來表示。以下是其偽代碼：

?? rl := al * bl ?? c := al * bh ?? rh := c >> 32 ?? c := c << 32 ?? rl := rl + c ?? if rl < c then rh++?? // 處理進位 ?? c := ah * bl ?? rh := rh + (c >> 32) ?? c := c << 32 ?? rl := rl + c ?? if rl < c then rh++?? // 處理進位 ?? rh := rh + ah * bh ?? return {rh, rl} }

乘法之后的取模運算可以用浮點運算快速完成。具體做法是乘積的高64位和低64位分別先對除數(shù)取模，然后再利用浮點單元合并取模。這里的浮點運算要求浮點單元以最高精度運算，計算前應(yīng)先將浮點單元控制字中的精度控制位設(shè)置為64位精度。為保證精度，應(yīng)當用80位浮點數(shù)實現(xiàn)此運算。偽代碼如下：

function mod64(rh, rl, n) { ?? rh := rh % n ?? rl := rl % n ?? r := fmodl(rh * 2 ** 64, n) ?? r := r + rl ?? if r < n then r := r - n??? // 處理進位 ?? r := fmodl(r, n) ?? return r }

至此，Rabin-Miller強偽素數(shù)測試的核心已經(jīng)全部實現(xiàn)。

?

二、Pollard r因數(shù)分解算法

Pollard r因數(shù)分解算法又稱為Pollard Monte Carlo因數(shù)分解算法。它的核心思想是：選取一個隨機數(shù)a。再選一個隨機數(shù)b。檢查是否大于1。若大于1，就是n的一個因子。若不大于1，再選取隨機數(shù)c，檢查和。如此繼續(xù)，直到找到n的一個非平凡因子。

它的主要實現(xiàn)方法是從某個初值出發(fā)，通過一個適當?shù)亩囗検?/span>進行迭代，產(chǎn)生一個偽隨機迭代序列直到其對模n出現(xiàn)循環(huán)。多項式的選擇直接影響著迭代序列的長度。經(jīng)典的選擇是選取，其中。不選擇0和的原因是避免當序列中某一項時后續(xù)各項全部為1的情況。迭代序列出現(xiàn)循環(huán)的期望時間和期望長度都與成正比。設(shè)n有一個非平凡因子p。由前面可知，迭代序列關(guān)于模p出現(xiàn)循環(huán)的期望時間和期望長度與成正比。當循環(huán)出現(xiàn)時，設(shè)，記，則d是p的倍數(shù)。當時，d就是n的一個非平凡因子。

function pollard_floyd(n) { ?? x := x0 ?? y := x0 ?? d := 1 ?? while d == 1 ?? { ????? x := f(x) ????? y := f(f(y))

但是在分解成功之前，p是未知的，也就無從直接判斷循環(huán)是否已經(jīng)出現(xiàn)。仍然設(shè)迭代序列關(guān)于模p出現(xiàn)循環(huán)，并設(shè)。由取模運算的性質(zhì)可知，即。故對任意，總有。記循環(huán)部分長度為t，若，則，那么。因此，只要對迭代序列中的偶數(shù)項和對應(yīng)的計算。只要，d就是n的一個非平凡因子。以上就是Pollard原來建議使用的Floyd循環(huán)判斷算法。由此得到Pollard r因數(shù)分解算法的第一個偽代碼：

function pollard_brent(n) { ?? x := x0 ?? y := x0 ?? k := 0 ?? i := 1 ?? d := 1 ?? while d == 1 ?? { ????? k := k + 1 ????? x := f(x) ????? d := gcd(x – y, n) ????? if 1 < d AND d < n then return d ????? if d == n then return FAILURE ????? if k == i then ????? { ??????? y := x ??????? i := i << 1 ????? } ?? } }

????? d := gcd(x – y, n) ????? if 1 < d AND d < n then return d ????? if d == n then return FAILURE ?? } }

由于循環(huán)出現(xiàn)的時空期望都是，Pollard r因數(shù)分解算法的總體時間復(fù)雜度也就是。在最壞情況下，其時間復(fù)雜度可能接近，但在一般條件下，時間復(fù)雜度都可以認為是。

?

參考資料：

Chris Caldwell “Fermat, probable-primality and pseudoprimes.”

Chris Caldwell “Strong probable-primality and a practical test.”

Eric W. Weisstein “Brent’s Factorization Method.”

Eric W. Weisstein “Pollard Rho Factorization Method.”

Eric W. Weisstein “Rabin-Miller Strong Pseudoprime Test.”

Eric W. Weisstein “Strong Pseudoprime.”

Unknown Author “Pollard’s r.”

下面附一個米勒羅賓算法和布萊德算法的代碼給大家當模板：

#include <stdio.h> #include <stdlib.h> #include <algorithm> #include <functional> const int MAX_COUNT = 6; unsigned __int64 allFactor[65], nFactor;unsigned __int64 fMultiModular( unsigned __int64 a, unsigned __int64 b, const unsigned __int64 n) {unsigned __int64 back = 0, temp = a % n;while ( b > 0 ){if ( b & 0x1 ){if ( (back = back + temp) > n )back -= n;}if ( (temp <<= 1) > n )temp -= n;b >>= 1;}return back; }unsigned __int64 GCD(unsigned __int64 a, unsigned __int64 b) {if(b==0)return a;else return GCD(b,a%b); } unsigned __int64 modular_exp(const unsigned __int64 &a, unsigned __int64 b, const unsigned __int64 &n) {unsigned __int64 d(1), dTemp(a % n);while ( b > 0 ){if ( b & 0x1 )d = fMultiModular(d, dTemp, n);dTemp = fMultiModular(dTemp, dTemp, n);b >>= 1;}return d; } bool fWitNess(const unsigned __int64 &a, const unsigned __int64 &n, char t, const unsigned __int64 &u) {unsigned __int64 x, y;x = modular_exp(a, u, n);while ( t-- ){y = fMultiModular(x, x, n);if ( y == 1 && x != 1 && x != n-1 )return true; //must notx = y;}return y != 1; } bool miller_rabin(const unsigned __int64 &n, int count) {if ( n == 1 )return false;if ( n == 2 )return true;unsigned __int64 a, u;char t(0);for (u = n-1; !(u & 0x1); u>>=1)++t;for (int i = 1; i <= count; ++i){a = rand() % (n-1) + 1;if ( fWitNess(a, n, t, u) )return false;}return true; } unsigned __int64 pollard_rho(const unsigned __int64 &c, const unsigned __int64 &num) {int i(1), k(2);unsigned __int64 x = rand() % num;unsigned __int64 y = x, comDiv;do{++i;if ( (x = fMultiModular(x, x, num) - c) < 0 )x += num;if ( x == y )break;comDiv = GCD((y-x+num)%num, num);if ( comDiv > 1 && comDiv < num )return comDiv;if ( i == k ){y = x;k <<= 1;}}while ( true );return num; } void fFindFactor(const unsigned __int64 &num) {if ( miller_rabin(num, MAX_COUNT) == true ){allFactor[++nFactor] = num;return;}unsigned __int64 factor;do{factor = pollard_rho(rand()%(num-1) + 1, num);}while ( factor >= num );fFindFactor(factor);fFindFactor(num/factor); }int main() {int t,i;unsigned __int64 test_num,min ,factor;scanf("%d",&t);while(t--){scanf("%I64u",&test_num);if(miller_rabin(test_num,5))printf("Prime\n");else{min = test_num;nFactor = 0;fFindFactor(test_num);for(i = 1; i <= nFactor ; i++){if(allFactor[i] < min)min = allFactor[i];}printf("%I64u\n",min);}} }

　　

?

?

轉(zhuǎn)載于:https://www.cnblogs.com/Chinese-Coder-Clarence/articles/2155449.html

總結(jié)

以上是生活随笔為你收集整理的64位以内Rabin-Miller 强伪素数测试和Pollard rho 因数分解解析的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。