实用算法实现-第 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