转素数的测试

 

一.费马小定里
if n is prime and (a,n) equals one ,then a^(n-1) = 1 (mod n)

费马小定理只是个必要条件,符合费马小定理而非素数的数叫做Carmichael.

前3个Carmichael数是561,1105,1729。

Carmichael数是非常少的。

在1~100000000范围内的整数中，只有255个Carmichael数。

为此又有二次探测定理,以确保该数为素数:

二.二次探测定理

二次探测定理 如果p是一个素数，0<x<p,则方程x^2≡1(mod p)的解为x=1,p-1

根据以上两个定理,如到Miller-Rabin算法的一般步骤:

0、先计算出m、j，使得n-1=m*2^j，其中m是正奇数，j是非负整数

1、随机取一个b，2<=b

2、计算v=b^m mod n

3、如果v==1，通过测试，返回

4、令i=1

5、如果v=n-1，通过测试，返回

6、如果i==j，非素数，结束

7、v=v^2 mod n，i=i+1

8、循环到5

注意
1 计算时会用到a^b mod n,应进行一定优化.
2 Miller-Rabin是随机算法
结果的正确率约为 75%,所以应该多次调用该函数,使正确概率提高为1-(1/4)^p

说明:

Miller-Rabin是随机算法

得到的结果的正确率为 75%,所以应该多次调用该函数,使正确概率提高为

1-(1/4)^p

 

//功能实现

 

void fFindFactor(const lltype &num) { if ( miller_rabin(num, MAX_COUNT) == true ) { allFactor[++nFactor] = num; return; } lltype factor; do { factor = pollard_rho(rand()%(num-1) + 1, num); }while ( factor >= num ); //printf("%I64d ",factor); fFindFactor(factor); fFindFactor(num/factor); } lltype pollard_rho(const lltype &c, const lltype &num) { int i(1), k(2); lltype x = rand() % num; lltype y = x, comDiv; do { ++i; if ( (x = fMultiModular(x, x, num) - c) < 0 ) x += num; if ( x == y ) break; comDiv = fGCD((y-x+num)%num, num); if ( comDiv > 1 && comDiv < num ) return comDiv; if ( i == k ) { y = x; k <<= 1; } }while ( true ); return num; } //miller-rabin法测试素数, count 测试次数 bool miller_rabin(const lltype &n, int count) { if ( n == 1 ) return false; if ( n == 2 ) return true; lltype a, u; //为了谨慎起见用了 lltype, 而没有用 int char t(0); //n-1 == (2^t)*u 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; } bool fWitNess(const lltype &a, const lltype &n, char t, const lltype &u) {//if n must not be prime return true, else return false lltype 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; } // d == a^b(mod n) lltype modular_exp(const lltype &a, lltype b, const lltype &n) { lltype d(1), dTemp(a % n);//当前二进制位表示的是进制数值 while ( b > 0 ) { if ( b & 0x1 ) d = fMultiModular(d, dTemp, n);//不直接 d*dTemp 是怕溢出 dTemp = fMultiModular(dTemp, dTemp, n); b >>= 1; } return d; } //back = (a*b)%num ( n <= 2^62) lltype fMultiModular(const lltype &a, lltype b, const lltype &n) { lltype back(0), temp(a % n); //b %= n; while ( b > 0 ) { if ( b & 0x1 ) { if ( (back = back + temp) > n ) back -= n; } if ( (temp <<= 1) > n ) temp -= n; b >>= 1; } return back; } lltype fGCD(lltype a, lltype b) { lltype c; while ( b ) { c = b; b = a%b; a = c; } return a; }