实用算法实现-第 28 篇 素数判别

28.1    朴素的素数判别

bool isPrime(__long n) { //简单的判断素数的确定性算法 __long i; if(n == 2|| n == 3) return true; if(n % 2 ==0) return false; for(i = 3;i < (__long)sqrt((double)n) + 1; i += 2){ if(n %i == 0) returnfalse; } return true; } int prime(int n) { int i; for(i=2;i*i<=n;i++) if(n%i==0) return0; return1; }
28.2    素数筛选

void Prime() { memset(a, 0, n*sizeof(a[0])); int num =0, i, j; for(i = 2;i < n; ++i) { if(!(a[i])) { p[num++] = i; } for(j =0; (j<num && i*p[j]<n); j++){ a[i*p[j]] = 1; if(!(i%p[j])) { break; } } } }
28.3    Miller-Rabin随机性素数测试方法

28.3.1   实例

PKU JudgeOnline, 1811, Prime Test.

28.3.2   问题描述

判断一个数N (2 < N < 2⁵⁴)是不是素数，是素数就输出 “Prim”，否则输出其最小的因子。

28.3.3   分析

这个题目考察的问题相当多，可以分为大素数的判别法、大数的分解两个问题。其中又牵涉到最大公约数的辗转相除求法（欧几里德算法）、模运算的性质（方幂模、模乘的实现）等。

这里N选取的范围恰到好处。输入的数不是很大，故此可以用__int64类型进行加、减、除、模运算。但是又要注意不能进行乘法运算，否则会溢出。

使用朴素的素数判别显然是不行的。这里需要使用概率算法：Miller-Rabin算法。其中Miller-Rabin需要求积的模和求乘方的模，这些都可以根据模的性质分解，使得计算不会溢出。

在完成素数判别之后，需要采用Pollard的rho启发式随机算法进行素数分解。该算法可能造成死循环。虽然死循环的可能性并不大，不幸碰到了就多提交几次。

最后完成的程序并不很优化，在JudgeOnline上提交只能用G++语言选项提交才能通过，耗时4000+ms。

28.3.4   程序

#include <stdio.h> #include <iostream> #include <stdlib.h> #include <time.h> #include <math.h> using namespace std; #define MAX 100000 //2147395600 = 46340^2 //2147483648 = 2^31 //2147488281 = 46341^2 #define __int64Max 3037000499 #define minNum 3037000499 //9223372030926249001 = 3037000499^2 //9223372036854775808 = 2^63 //9223372037000250000 = 3037000500^2 //注意int小于 #define maxNum 3037100499 #define __long __int64 #define repeatTime 50 __longModularMultiply(__long a, __long b, __long p) { //求积的模s=a*bmod p, 注意先求a*b可能会溢出 //a*b = a*(b/2)*2+ a*(b%2) (% c) //用同样的方法求a*(b/2)%c if(b >__int64Max) { return(ModularMultiply(a, b / 2, p) * 2 + ModularMultiply(a, b % 2, p)) % p; } if(a >__int64Max) { return(ModularMultiply(a / 2, b, p) * 2 + ModularMultiply(a % 2, b, p)) % p; // return(ModularMultiply(a / 2, b, p) * 2 + (a % 2) * b)) % p;//这个更简化，不过好像没必要 } return (a *b) % p; } __longModularExponent(__long a1, __long j, __long p) { //求方幂模s=a^jmod p, 注意先求a^j可能会溢出 __long a = a1; __long s = 1; //a为__int64就可保证正确计算(a * a)和(s * a) //故s不必为__int64 while(j> 0){ if((j %2) == 1) { s = ModularMultiply(s, a, p); } a = ModularMultiply(a, a, p); j /= 2; } return s; } bool Btest(__long a, __long n) { //n为奇数，返回真说明n是强伪素数 __long s = 0; __long t = n - 1; __long x; __long x1; //x最大不可能超过n //为了确保(x * x)计算正确采用__int64 //注意用int不一定会导致错误结果 int i; while((t %2) != 1){ s++; t /= 2; } x = ModularExponent(a,t,n); x1 = x; for(i = 1;i <= s; i++){ x1 = ModularMultiply(x, x, n); if(x1== 1 && x != 1 && x != n - 1) { returnfalse; } x = x1; } if(x1 != 1) return false; // cout<<"suspect: " << a << endl; return true; } bool Btest1(__long a, __long n) { //n为奇数，返回真说明n是强伪素数 __long s = 0; __long t = n - 1; __long x; //x最大不可能超过n //为了确保(x * x)计算正确采用__int64 //注意用int不一定会导致错误结果 __long i; while((t %2) != 1){ s++; t /= 2; } // cout << n - 1<< " = 2^" << s << " * " << t<< endl; x = ModularExponent(a, t, n); if(x ==1||x == n-1) { // cout<<"suspect1: " << a << endl; return true; } for(i = 0;i < s; i++){ x = ModularMultiply(x, x, n); if(x ==n-1) { // cout<<"suspect2: " << a << endl; returntrue; } } return false; } bool MillRab(__long n) { //奇n>4，返回真时表示素数，假表示合数 if(n == 2|| n == 3) return true; if(n % 2 ==0) return false; __long temp = rand(); __long a = (temp % (n - 3)) + 2; returnBtest(a, n); } bool RepeatMillRab(__long n, intk) { int i; for(i = 0;i < k; i++){ if(MillRab(n)== false) { returnfalse; } } return true; } __longfactorMin(__long n) { __long i; for(i=2;i*i<=n;i++){ if(n %i==0) { returni; } } return 0; } __longEuclid(__long a, __long b) { /* 欧几里德算法 GCD递归定理： 对任意非负证书a和任意正整数b gcd(a, b) = gcd(b, a mod b) */ if(b == 0){ returna; }else{ returnEuclid(b, a % b); } } __longPollardRho(__long n) { __long temp; __long x; __long y; __long i; __long k; __long d; i = 1; temp = rand(); x = temp % n; y = x; k = 2; while(1){ /* if(i % 100 ==0) { cout<< i << endl; }*/ i++; x = (ModularMultiply(x, x, n) + n - 1)% n; temp = y - x; if(temp< 0) { temp = -temp; } d = Euclid(temp, n); if(d !=1 && d != n) { returnd; } if( i== k) { y = x; k = 2 * k; } } } int main() { srand((unsignedint)time(NULL)); __long a; __long temp; __long temp1; __long min; int i, j; int num; __long queue[1000]; int top; cin >> num; for(i = 0;i < num; i++){ cin >> a; if(a ==2) { cout << "Prime" << endl; continue; } if(a %2 == 0) { cout << "2" << endl; continue; } top = 0; queue[top++] = a; min = a; for(j =0; j < top; j++){ temp = queue[j]; // cout<< "分解："<< temp << endl; if(temp< MAX){ temp1 = factorMin(temp); if(temp1!= 0){ queue[top++] = temp1; if(min> temp1) { min = temp1; } } }elseif(RepeatMillRab(temp, 50) != true){ // cout<< temp << "是一个合数" << endl; temp1 = PollardRho(temp); queue[top++] = temp1; if(min> temp1) { min = temp1; } // cout<< "得到新的因子：" << temp1; temp1 = temp / temp1; queue[top++] = temp1; // cout<< " 和" << temp1 << endl; if(min > temp1) { min = temp1; } } } if(min== a){ cout <<"Prime" <<endl; }else{ cout << min << endl; } } return 1; }
28.4    实例

PKU JudgeOnline, 3126, Prime Path.

PKU JudgeOnline, 3006, Dirichlet's Theoremon Arithmetic Progressions.

PKU JudgeOnline, 2739, Sum of ConsecutivePrime Numbers.

本文章欢迎转载，请保留原始博客链接http://blog.csdn.net/fsdev/article