millrab.c

高概率找到正确解 基本思想：为了增加一个一致的P正确算法成功的概率

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#ifndef MAX_VOLUME

#define MAX_VOLUME 10000

#endif

/*** definite algorithm for prime number test ***/
int DefPrime(int n)
{
 int i;
 for(i = 2; i <= (int)sqrt(n); i++){
   if((n % i) == 0) 
     return 0;
 }
 return 1;
}

/*** get s = pow(a, t) % n, prevent the overflow ***/
int ModularExponent(int a, int t, int n)
{
 int s = 1;
 
 while(t > 0){
   if((t % 2) == 1)
      s = (s * a) % n;
   a = (a * a) % n;
   t = t / 2; 
 }
 return s;
}

/*** random algorithm for prime number test ***/
int Btest(int a, int n)
{
 int s, t, x, i;

 s = 0, t = n - 1;
 do{
   s++; t = t/2;
 }while(t % 2 != 1);
 
 x = ModularExponent(a, t, n);
 if(x == 1 || x == n - 1)
   return 1;
 for(i = 1; i <= s - 1; i++){
   x = (x * x) % n;
   if(x == n - 1)
     return 1;
 }
 return 0;
}

int MillRab(int n)
{
 int a;

 a = rand() % (n - 3) + 2;
 return Btest(a, n); 
}

int RepeatMillRab(int n, int k)
{
 int i;

 srand((unsigned)time( NULL ));
 for(i = 1; i <= k; i++)
   if(MillRab(n) == 0)
     return 0;
 return 1;
}

/*** store prime numbers lower than <m> in arrays <primes1> <primes2> got 
     by algorthms DefPrime and RepeatMillRab respectly ***/
int GetPrimes(int m, int primes1[], int primes2[], int *count1, int *count2)
{
 int n, i;

 primes1[0] = primes2[0] = 2, primes1[1] = primes2[1] = 3;
 *count1 = *count2 = 2;

 n = 5, i = 2;
 do{
   if(DefPrime(n)){
     primes1[i] = n;
     (*count1)++;
     i++;
   }
   n = n + 2;
 }while(n < m);

 n =5, i = 2;
 do{
   if(RepeatMillRab(n, (int)log(n))){
     primes2[i] = n;
     (*count2)++;
     i++;
   }
   n = n + 2;

 }while(n < m);
 return 0;
}

int main(int argc, char **argv)
{
 int primes1[MAX_VOLUME], primes2[MAX_VOLUME];
 int count1 = 0, count2 = 0, m, i, j;
 double prob = 0.0;

 if(argc != 2){
   printf("Usage: %s <n>          get primes lower than <n>\n", argv[0]);
   exit(1);
 }
 m = atoi(*(argv+1));

 GetPrimes(m, primes1, primes2, &count1, &count2);
/*
 for(i = 0; i < count2; i++){
   for(j = 0; j < count1; j++){
     if(primes1[j] == primes2[i])
       break;
   }
   if(j >= count1)
     error++;
 }

 for(i= 0; i < count1; i++){
   printf("%d ", primes1[i]);
   if((i + 1) % 10 == 0)
     printf("\n");
 } 
 printf("\n");
 for(i= 0; i < count2; i++){
   printf("%d ", primes2[i]);
   if((i + 1) % 10 == 0)
     printf("\n");
 } 
 printf("\n");
*/
 prob = (double)(count2-count1) / count2;
 printf("m: %d\n", m);
 printf("count1: %d, count2: %d\n", count1, count2);
 printf("The error probability: %f\n", prob);
 if(prob < 0)
     printf("Overflowed during computing!\n");
 return 0;
}