实用算法实践-第 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 < 254)是不是素数，是素数就输出 “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