64位以内Rabin-Miller 强伪素数测试和Pollard rho 因数分解解析
在求解POJ1811题Prime Test中应用到的两个重要算法是Rabin-Miller强伪素数测试和Pollard r因数分解算法。前者可以在
一、Rabin-Miller强伪素数测试
Rabin-Miller强伪素数测试的基本思想来源于如下的Fermat小定理:
如果p是一个素数,则对任意a有
利用Fermat小定理可以得到一个测试合数的有力算法:对
这个算法的成功率是相当高的。在小于25,000,000,000的1,091,987,405个素数中,一共只用21,853个以2为基的伪素数。但不幸的是,Alford、Granville和Pomerance在1994年证明了存在无穷多个被称为Carmichael数的整数对于任意与其互素的整数a算法的计算结果都是1。最小的五个Carmichael数是561、1,105、1,729、2,465和2,801。
考虑素数的这样一个性质:若n是素数,则1对模n的平方根只可能是1和
(Rabin-Miller强伪素数测试)记
这个测试同时被冠以Miller的名字是因为Miller提出并证明了如下测试:如果扩展黎曼猜想(extended Riemann hypothesis)成立,那么对于所有满足
尽管Monier和Rabin在1980年证明了这个测试的错误概率(即合数通过测试的概率)不超过
取前k个素数为基,并用
考虑到64位二进制数能表示的范围,只需取前9个素数为基,则对小于
Rabin-Miller强伪素数测试本身的形式稍有一些复杂,在实现时可以下面的简单形式代替:
对
代替的理由可简单证明如下:
仍然记
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 } |
function
mul64to128(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 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} } |
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 } |
二、Pollard r因数分解算法
Pollard r因数分解算法又称为Pollard Monte Carlo因数分解算法。它的核心思想是:选取一个随机数a。再选一个随机数b。检查
它的主要实现方法是从某个初值
function
pollard_floyd(n) { x := x0 y := x0 d := 1 while d == 1 { x := f(x) y
:= f(f(y)) |
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 } } |
由于循环出现的时空期望都是
参考资料:
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 not x = 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); } } }
posted on 2011-08-27 13:16 _Clarence 阅读(1106) 评论(0) 编辑 收藏 举报