// hexadeca program written by Robert Gerbicz
// Searching for 16 primes: k*2^n+-1, k*2^(n+1)+-1, k*2^(n+2)+-1, k*2^(n+3)+-1, 2^n+-k, 2^(n+1)+-k, 2^(n+2)+-k, 2^(n+3)+-k
// version 1.0
// for hexadecaproth internet search see: http://mersenneforum.org/forumdisplay.php?f=63
// you can complie it under Pentium4 by the following optimization flags:
// gcc -g -O2 -fomit-frame-pointer -march=pentium4 -mcpu=pentium4 hexadeca_1_0.c -o hexadeca_1_0 .libs/libgmp.a
// or for Athlon-xp:
// gcc -g -O2 -fomit-frame-pointer -march=athlon-xp -mcpu=athlon-xp hexadeca_1_0.c -o hexadeca_1_0 .libs/libgmp.a
// or for Athlon-mp:
// gcc -g -O2 -fomit-frame-pointer -march=athlon-mp -mcpu=athlon-mp hexadeca_1_0.c -o hexadeca_1_0 .libs/libgmp.a
// and so on...

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "gmp.h"

double VERSION = 1.0;  // version number
char *progname;

void
print_usage_and_exit ()
{
  fprintf (stderr, "usage: %s startworkunit endworkunit.\n Here 0<=startworkunit<=endworkunit<3869775\n", progname);
  exit (-1);
}

main (int argc, char *argv[])

{
// define parameters for the program, don't change these parameters!
// we'll searching for hexadecaproth only for n=76
// the expected number of hexadecaproths for this n is 16
// There are 3869775 workunits.
// We'll find also about 6728 dodecaproths for this n, but not all of them
// ( there are about 773157 dodecaproths for this n ) 
// because this is optimized for hexadecaproth search.
	
  unsigned long int sieve_limit=1000;
  unsigned long int number_of_primes=168;  // number of primes up to sieve_limit
  unsigned long int block_length=23223;
  unsigned char n=76;
  unsigned char x=15;  // we will using 15 primes: 2,3,5,7,11,13,17,19,23,37,41,43,47,67,71

  progname = argv[0];

  time_t seconds;

  if (argc!=3)  print_usage_and_exit ();

  double status=0.0,status_new;
  
  unsigned long int start_workunit=atoi(argv[1]);
  unsigned long int end_workunit=atoi(argv[2]);

// checking parameters
 
  if(start_workunit>end_workunit)
  print_usage_and_exit ();
  
  if(end_workunit>=3869775)
  print_usage_and_exit();

// print n, start_workunit, end_workunit and version to the screen and to results_hexadeca.txt file
      
  FILE *found;

  printf("You can also find the results in results_hexadeca.txt file ( These are 3-probable primes )\n");
  printf("If the program find some dodecaproths then these will be saved to results_ex.txt file, but we won't display these\n");
  gmp_printf("   n=%d, start_workunit=%d, end_workunit=%d, version=%.1f\n",n,start_workunit,end_workunit,VERSION);
  printf("Starting the sieve...\n");

  found=fopen("results_hexadeca.txt","a+");
  gmp_fprintf(found,"   n=%d, start_workunit=%d, end_workunit=%d, version=%.1f\n Starting the sieve...\n",n,start_workunit,end_workunit,VERSION);
  fclose(found);

  unsigned char mul_rem[7][25]={{0,1,10,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
	                         {0,5,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
	                         {0,3,5,6,7,10,11,12,14,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
	                         {0,1,2,9,10,17,18,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
	                         {0,3,5,7,8,9,10,13,14,15,16,18,20,0,0,0,0,0,0,0,0,0,0,0,0},
	                         {0,1,2,3,4,6,8,9,12,14,18,19,23,25,28,29,31,33,34,35,36,0,0,0,0},
	                         {0,1,3,6,7,9,11,12,13,14,15,17,19,22,24,26,27,28,29,30,32,34,35,38,40}};

  unsigned char mul[7]={11,13,17,19,23,37,41};
  
  unsigned char length[7]={3,3,9,7,13,21,25};

  unsigned long int i,q,v,num_primes=0;
  unsigned short table[18*number_of_primes+18];
  unsigned char b[block_length];
  unsigned long long int j;

  register unsigned short J;
  register unsigned short A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12,A13,A14,A15;
  register unsigned short prm;

  for(i=0;i<=block_length;i++) b[i]=1;

// sieving up to sieve_limit
// there are number_of_primes up to sieve_limit

  for(i=2;i<=sieve_limit;i++)
  {
  if(b[i]==1) {
  num_primes++;
  table[18*num_primes]=i;
  for(j=i*i;j<=sieve_limit;j+=i) b[j]=0;
  }
  }

  if(num_primes!=number_of_primes) 
  {
	  printf("An error occured in the program!\n");
	  exit(-1);
  }
  
// positions in table: prm    0
//                     z      1
//                     first rem 2
//                     rem's diff  3,4,5,6,7,8,9,10,11,12,13,14,15,16,17

// create rem table

  signed long long int h;
  unsigned long int m=n;
  unsigned char bits=0;

  while(m>0)  bits++,b[bits]=m%2,m=m/2;

  for(i=2;i<=number_of_primes;i++)
  {
  v=18*i;
  q=table[v],h=2;

  for(j=bits-1;j>=1;j--)
  {
  h=(h*h)%q;
  if(b[j]) h=2*h;
  if(h>=q) h-=q;
  }

  table[v+2]=h;
  table[v+3]=q-table[v+2];
  table[v+4]=2*table[v+2];
  if(table[v+4]>=q) table[v+4]-=q;
  table[v+5]=q-table[v+4];
  table[v+6]=2*table[v+4];
  if(table[v+6]>=q) table[v+6]-=q;
  table[v+7]=q-table[v+6];
  table[v+8]=2*table[v+6];
  if(table[v+8]>=q) table[v+8]-=q;
  table[v+9]=q-table[v+8];
  
  h=(q+1)/2;
  for(j=bits-1;j>=1;j--)
  {
  h=(h*h)%q;
  if(b[j])
  {
  if(h&1) h=(h+q)/2;
  else h=h/2;
  }
  }
  table[v+10]=h;
  table[v+11]=q-table[v+10];
  if(table[v+10]&1) table[v+12]=(table[v+10]+q)/2;
  else table[v+12]=table[v+10]/2;
  table[v+13]=q-table[v+12];
  if(table[v+12]&1) table[v+14]=(table[v+12]+q)/2;
  else table[v+14]=table[v+12]/2;
  table[v+15]=q-table[v+14];
  if(table[v+14]&1) table[v+16]=(table[v+14]+q)/2;
  else table[v+16]=table[v+14]/2;
  table[v+17]=q-table[v+16];
  }
  
// replace 29 by 67 and 31 by 71:
  j=18*19-1;
  for(i=18*10;i<18*11;i++) j++,v=table[j],table[j]=table[i],table[i]=v;
  
  j=18*20-1;
  for(i=18*11;i<18*12;i++) j++,v=table[j],table[j]=table[i],table[i]=v;

// replace 67 by 37 and 71 by 41:
  j=18*12-1;
  for(i=18*10;i<18*11;i++) j++,v=table[j],table[j]=table[i],table[i]=v;
  
  j=18*13-1;
  for(i=18*11;i<18*12;i++) j++,v=table[j],table[j]=table[i],table[i]=v;

// so now the order: 2,3,5,7,11,13,17,19,23,37,41,67,71,43,47  
      
  for(j=0;j<=block_length;j++) b[j]=1;

  mpz_t step_gmp;
  mpz_t s_gmp;
  mpz_t small_step_gmp;
  mpz_init(step_gmp);
  mpz_init(s_gmp);
  mpz_init(small_step_gmp);
  
  mpz_set_ui(step_gmp,1);
  for(i=1;i<=15;i++) mpz_mul_ui(step_gmp,step_gmp,table[18*i]);  // step_gmp=2*3*5*7*11*13*17*19*23*37*41*43*47*67*71
  
  mpz_divexact_ui(small_step_gmp,step_gmp,43*47*67*71);  // small_step_gmp=2*3*5*7*11*13*17*19*23*37*41
    
// small rem table for small primes

  unsigned char srem[x+1][table[18*x]+1];

  for(i=2;i<=x;i++)
  {
  for(j=0;j<table[18*i];j++) srem[i][j]=1;
  for(j=1;j<=16;j++) srem[i][table[18*i+j+1]]=0;  // we'll sieving 16 forms
  }

  srem[1][0]=0;  // for the first prime=2
  srem[1][1]=1;

// computation of table[18*i+1]=(-small_step_gmp*modinverse(step_gmp,table[18*i])) mod table[18*i]
// Here table[18*i] is the i-th prime number

  for(j=x+1;j<=number_of_primes;j++)
  {
  mpz_set_ui(s_gmp,table[18*j]);
  mpz_invert(s_gmp,step_gmp,s_gmp);
  mpz_mul(s_gmp,s_gmp,small_step_gmp);
  mpz_neg(s_gmp,s_gmp);
  table[18*j+1]=mpz_mod_ui(s_gmp,s_gmp,table[18*j]);
  }

  unsigned char test;
  signed long int s0;
  
  mpz_t k_gmp;      // this is the multiplier
  mpz_t three_gmp;  // this is 3
  mpz_t g_gmp;  // g_gmp is equal to g
  mpz_t pr_gmp;
  mpz_t num_prp_gmp;
  mpz_init(k_gmp);
  mpz_init(g_gmp);
  mpz_init(three_gmp);
  mpz_init(pr_gmp);
  mpz_init(num_prp_gmp);

  mpz_set_ui(three_gmp,3);
  mpz_set_ui(num_prp_gmp,0);

  unsigned char f;
  signed long int z;
  unsigned long int F,R;
  signed long int A;
    
  unsigned long int v17,v2;
  signed long int M,start;

  printf("Using %d primes to reduce the size of the sieve array\n",x);

  found=fopen("results_dodeca.txt","a+");
  gmp_fprintf(found,"Using the first %d primes to reduce the size of the sieve array\n",x);
  fclose(found);

  unsigned char res[x];
  unsigned char sres[x];
  unsigned char sprm[x];
  unsigned long int work;

  unsigned long int s,y,work_s;
  unsigned char D;
  unsigned long int L,good_L,diff;
  unsigned long long int diff_L;
        
  for(i=12;i<=x;i++) sprm[i]=table[18*i];  // sprm[i]=prm(i)
  for(i=12;i<=x;i++) sres[i]=mpz_mod_ui(s_gmp,small_step_gmp,sprm[i]);  // sres[i]=small_step mod prm(i)
  
  mpz_t prod_gmp;
  mpz_t r_gmp;
  mpz_init(prod_gmp);
  mpz_init(r_gmp);
  
  seconds = time (NULL);  // start the time

// this is the main sieve + prp

  for(work=start_workunit;work<=end_workunit;work++) {

//  initialize work
	  
  mpz_set_ui(g_gmp,105);
  mpz_set_ui(prod_gmp,210);
  s=3869775;
  work_s=work;
  good_L=0;

  for(i=0;i<=6;i++)  {
  s=s/length[i];
  D=mul_rem[i][work_s/s];
  work_s=work_s%s;
  prm=mul[i];
  mpz_set_ui(s_gmp,prm);
  mpz_invert(s_gmp,prod_gmp,s_gmp);
  mpz_set_ui(r_gmp,D);
  mpz_sub(r_gmp,r_gmp,g_gmp);
  mpz_mul(s_gmp,s_gmp,r_gmp);
  y=mpz_mod_ui(s_gmp,s_gmp,prm);
  mpz_addmul(g_gmp,s_gmp,prod_gmp);
  mpz_mul_ui(prod_gmp,prod_gmp,prm);
  }

// determine starting remainders
    
  for(i=x+1;i<=number_of_primes;i++)
  {
     v=18*i;
     prm=table[v];
     for(j=0;j<=15;j++)
     {
	     mpz_set_ui(s_gmp,table[v+2+j]);
	     mpz_sub(s_gmp,s_gmp,g_gmp);
	     mpz_set_ui(r_gmp,prm);
	     mpz_invert(r_gmp,step_gmp,r_gmp);
	     mpz_mul(s_gmp,s_gmp,r_gmp);
	     table[v+2+j]=mpz_mod_ui(s_gmp,s_gmp,prm);
     }
  }

  mpz_sub(g_gmp,g_gmp,small_step_gmp);

  for(i=12;i<=15;i++)  res[i]=mpz_mod_ui(s_gmp,g_gmp,sprm[i]);
      
// sorting the remainders in the table, and then save rem's difference table[j]=table[j]-table[18*i+2] for j>v

  for(i=x+1;i<=number_of_primes;i++)
  {
     v=18*i+2;
     for(j=v;j<=v+15;j++)
     {
        for(h=v;h<=v+15;h++)
     	     if(table[j]<table[h]) F=table[j],table[j]=table[h],table[h]=F;
     }
     for(j=v+1;j<=v+15;j++)  table[j]-=table[v];
  }
	  
  for(L=0;L<9613897;L++)  // 9613897=43*47*67*71
  {
  test=1,f=12;

  mpz_add(g_gmp,g_gmp,small_step_gmp);
  
  for(i=12;i<=x;i++) {
  res[i]+=sres[i];
  if(res[i]>=sprm[i]) res[i]-=sprm[i];
  }

  while(test&(f<=x))  test=srem[f][res[f]],f++;

  if(test)
  {
        diff_L=L-good_L;
	  good_L=L;
        v=18*x;
	  for(i=x+1;i<=number_of_primes;i++)
      {
         v+=18;
         prm=table[v];
         diff=(((unsigned long long) table[v+1])*diff_L)%prm;
         table[v+2]+=diff;
         if(table[v+2]>=prm) table[v+2]-=prm;
      }

      v=18*x;
      for(i=x+1;i<=number_of_primes;i++)
      {
          v+=18;
          prm=table[v];
          start=table[v+2];

          v2=v+2,z=v17=v+17,s0=prm-start;
          while((z>v2)&(table[z]>=s0))  b[table[z]-s0]=0,z--;

          A1=table[v+3];
          A2=table[v+4];
          A3=table[v+5];
          A4=table[v+6];
          A5=table[v+7];
          A6=table[v+8];
          A7=table[v+9];
          A8=table[v+10];
          A9=table[v+11];
          A10=table[v+12];
          A11=table[v+13];
          A12=table[v+14];
          A13=table[v+15];
          A14=table[v+16];
          A15=table[v+17];

             for(J=start;J<=block_length-prm;J+=prm)   b[J]=0,b[J+A1]=0,b[J+A2]=0,b[J+A3]=0,b[J+A4]=0,b[J+A5]=0,b[J+A6]=0,b[J+A7]=0,b[J+A8]=0,b[J+A9]=0,b[J+A10]=0,b[J+A11]=0,b[J+A12]=0,b[J+A13]=0,b[J+A14]=0,b[J+A15]=0;

             b[J]=0,M=block_length-J;

	  while((v2<v17)&(table[v2+1]<M)) v2++,b[table[v2]+J]=0;
      }

      for(i=0;i<block_length;i++)
      {
      if(b[i])
      {
      test=1;

      mpz_mul_ui(k_gmp,step_gmp,i);
      mpz_add(k_gmp,k_gmp,g_gmp);

         // 3-prp test for 2^n+k
         mpz_set_ui(pr_gmp,1);
         mpz_mul_2exp(pr_gmp,pr_gmp,n);
         mpz_add(pr_gmp,pr_gmp,k_gmp);
         mpz_powm(s_gmp,three_gmp,pr_gmp,pr_gmp);
         if(mpz_cmp_ui(s_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);

      if(test)
      {  // 3-prp test for 2^(n+1)+k
         mpz_mul_2exp(pr_gmp,pr_gmp,1);
         mpz_sub(pr_gmp,pr_gmp,k_gmp);
         mpz_powm(s_gmp,three_gmp,pr_gmp,pr_gmp);
	   if(mpz_cmp_ui(s_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }

      if(test)
      {  // 3-prp test for 2^(n+2)+k
         mpz_mul_2exp(pr_gmp,pr_gmp,1);
         mpz_sub(pr_gmp,pr_gmp,k_gmp);
         mpz_powm(s_gmp,three_gmp,pr_gmp,pr_gmp);
	   if(mpz_cmp_ui(s_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }
      
      if(test)
      {  // 3-prp test for 2^n-k
         mpz_set_ui(pr_gmp,1);
         mpz_mul_2exp(pr_gmp,pr_gmp,n);
         mpz_sub(pr_gmp,pr_gmp,k_gmp);
         mpz_powm(s_gmp,three_gmp,pr_gmp,pr_gmp);
         if(mpz_cmp_ui(s_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }

      if(test)
      {  // 3-prp test for 2^(n+1)-k
         mpz_mul_2exp(pr_gmp,pr_gmp,1);
         mpz_add(pr_gmp,pr_gmp,k_gmp);
         mpz_powm(s_gmp,three_gmp,pr_gmp,pr_gmp);
         if(mpz_cmp_ui(s_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }
      
     if(test)
      {  // 3-prp test for 2^(n+2)-k
         mpz_mul_2exp(pr_gmp,pr_gmp,1);
         mpz_add(pr_gmp,pr_gmp,k_gmp);
         mpz_powm(s_gmp,three_gmp,pr_gmp,pr_gmp);
         if(mpz_cmp_ui(s_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }
      
      if(test)
      {  // 3-prp test for k*2^n+1
         mpz_mul_2exp(pr_gmp,k_gmp,n);
         mpz_add_ui(pr_gmp,pr_gmp,1);
         mpz_powm(s_gmp,three_gmp,pr_gmp,pr_gmp);
         if(mpz_cmp_ui(s_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }
 
      if(test)
      {  // 3-prp test for k*2^(n+1)+1
         mpz_mul_2exp(pr_gmp,pr_gmp,1);
         mpz_sub_ui(pr_gmp,pr_gmp,1);
         mpz_powm(s_gmp,three_gmp,pr_gmp,pr_gmp);
         if(mpz_cmp_ui(s_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }

      if(test)
      {  // 3-prp test for k*2^(n+2)+1
         mpz_mul_2exp(pr_gmp,pr_gmp,1);
         mpz_sub_ui(pr_gmp,pr_gmp,1);
         mpz_powm(s_gmp,three_gmp,pr_gmp,pr_gmp);
         if(mpz_cmp_ui(s_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }
      
      if(test)
      {  // 3-prp test for k*2^n-1
         mpz_mul_2exp(pr_gmp,k_gmp,n);
         mpz_sub_ui(pr_gmp,pr_gmp,1);
         mpz_powm(s_gmp,three_gmp,pr_gmp,pr_gmp);
         if(mpz_cmp_ui(s_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }
      
      if(test)
      {  // 3-prp test for k*2^(n+1)-1
         mpz_mul_2exp(pr_gmp,pr_gmp,1);
         mpz_add_ui(pr_gmp,pr_gmp,1);
         mpz_powm(s_gmp,three_gmp,pr_gmp,pr_gmp);
         if(mpz_cmp_ui(s_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }

      if(test)
      {  // 3-prp test for k*2^(n+2)-1
         mpz_mul_2exp(pr_gmp,pr_gmp,1);
         mpz_add_ui(pr_gmp,pr_gmp,1);
         mpz_powm(s_gmp,three_gmp,pr_gmp,pr_gmp);
         if(mpz_cmp_ui(s_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }
      
	 if(test)
	 { // we have found a dodecaproth! Save the result in results_ex.txt file.
	   found=fopen("results_ex.txt","a+");
	   gmp_fprintf(found,"   %Zd  %d\n",k_gmp,n);
	   fclose(found);
	 }
	 
	 if(test)
	 {  // 3-prp test for 2^(n+3)+k
	   mpz_set_ui(pr_gmp,1);
	   mpz_mul_2exp(pr_gmp,pr_gmp,n+3);
	   mpz_add(pr_gmp,pr_gmp,k_gmp);
         mpz_powm(s_gmp,three_gmp,pr_gmp,pr_gmp);
         if(mpz_cmp_ui(s_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }
      
      if(test)
      { // 3-prp test for 2^(n+3)-k	 
         mpz_set_ui(pr_gmp,1);
	   mpz_mul_2exp(pr_gmp,pr_gmp,n+3);
	   mpz_sub(pr_gmp,pr_gmp,k_gmp);
         mpz_powm(s_gmp,three_gmp,pr_gmp,pr_gmp);
         if(mpz_cmp_ui(s_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
       }
       
       if(test)
       { // 3-prp test for k*2^(n+3)+1
         mpz_mul_2exp(pr_gmp,k_gmp,n+3);
         mpz_add_ui(pr_gmp,pr_gmp,1);
         mpz_powm(s_gmp,three_gmp,pr_gmp,pr_gmp);
         if(mpz_cmp_ui(s_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
       }
 
	   if(test)
       { // 3-prp test for k*2^(n+3)-1
	   mpz_sub_ui(pr_gmp,pr_gmp,2);
	   mpz_powm(s_gmp,three_gmp,pr_gmp,pr_gmp);
         if(mpz_cmp_ui(s_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
       }
	 
       if(test)
       { // we have found a hexadecaproth!
	     found=fopen("results_hexadeca.txt","a+");
	     gmp_fprintf(found,"   %Zd  %d\n",k_gmp,n);
	     fclose(found);
	     found=fopen("results_magic.txt","a+");
	     gmp_fprintf(found,"   %Zd  %d\n",k_gmp,n);
	     fclose(found);     
	     gmp_printf("   %Zd  %d\n",k_gmp,n);
       }	 
      }
      b[i]=1;
      }

      status_new=((double)((work-start_workunit)*9613897+L)/((double)((end_workunit-start_workunit+1))*9613897.0)*100.0);

      if(status_new-status>0.1) 
      {
      status=status_new;
      gmp_printf("   Status: %.1f percentage of the project is complete. Time thusfar: %ld sec. \r\r",status,time(NULL)-seconds),fflush(stdout);
      }      
  }
  }
  }
 	   printf("                                                                                                                               \r\r"),fflush(stdout);

// The sieving is complete

   gmp_printf("The sieving is complete.\nNumber of Prp tests=%Zd\nTime=%ld sec.\n",num_prp_gmp,time(NULL)-seconds);

   found=fopen("results_hexadeca.txt","a+");
   gmp_fprintf(found,"The sieving is complete.\nNumber of Prp tests=%Zd\nTime=%ld sec.\n",num_prp_gmp,time(NULL)-seconds);
   fclose(found);

// clear gmp variables
 
   mpz_clear(step_gmp);
   mpz_clear(small_step_gmp);
   mpz_clear(s_gmp);
   mpz_clear(r_gmp);
   mpz_clear(k_gmp);
   mpz_clear(g_gmp);
   mpz_clear(prod_gmp);
   mpz_clear(three_gmp);
   mpz_clear(pr_gmp);
   mpz_clear(num_prp_gmp);

   return 0;
}