// dodeca program written by Robert Gerbicz
// Searching for 12 primes: k*2^n+-1, 2^n+-k, k*2^(n+1)+-1, 2^(n+1)+-k, k*2^(n+2)+-1, 2^(n+2)+-k
// version 2.0
// for octoproth and dodecaproth 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 dodeca_2_0.c -o dodeca_2_0 .libs/libgmp.a
// or for Athlon-xp:
// gcc -g -O2 -fomit-frame-pointer -march=athlon-xp -mcpu=athlon-xp dodeca_2_0.c -o dodeca_2_0 .libs/libgmp.a
// or for Athlon-mp:
// gcc -g -O2 -fomit-frame-pointer -march=athlon-mp -mcpu=athlon-mp dodeca_2_0.c -o dodeca_2_0 .libs/libgmp.a
// and so on...

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

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

void
print_usage_and_exit ()
{
  fprintf (stderr, "usage: %s n kmin kmax\n", progname);
  exit (-1);
}

main (int argc, char *argv[])

{
  progname = argv[0];

  time_t seconds;

  if (argc != 4)  print_usage_and_exit ();

  mpz_t kmin_gmp;
  mpz_t kmax_gmp;
  mpz_init(kmin_gmp);
  mpz_init(kmax_gmp);

  signed long int n = atoi (argv[1]);
  signed long long int s,kmin,kmax;
  double status=0.0;

  mpz_set_str(kmin_gmp,argv[2],10);
  mpz_set_str(kmax_gmp,argv[3],10);

// checking parameters

  if(n<25)
  {
  printf("n should be at least 25\n");
  print_usage_and_exit ();
  }
  if(n>=100)
  {
  printf("n is too large\n");
  print_usage_and_exit ();
  }
  if(mpz_sgn(kmin_gmp)<0)
  {
  printf("kmin should be nonnegative\n");
  print_usage_and_exit ();
  }
  s=mpz_sizeinbase(kmax_gmp,2);
  if(s>60)
  {
  printf("kmax should be less than 2^60\n");
  print_usage_and_exit ();
  }
  if(mpz_cmp(kmin_gmp,kmax_gmp)>0)
  {
  printf("kmax should be at least kmin\n");
  print_usage_and_exit ();
  }

  FILE *found;

  printf("You can also find the k n values in results_dodeca.txt file ( These are 3-probable primes )\n");
  gmp_printf("   n=%d, kmin=%Zd, kmax=%Zd, version=%.1f\n",n,kmin_gmp,kmax_gmp,VERSION);
  printf("Starting the sieve...\n");

  found=fopen("results_dodeca.txt","a+");
  gmp_fprintf(found,"   n=%d, kmin=%Zd, kmax=%Zd, version=%.1f\n Starting the sieve...\n",n,kmin_gmp,kmax_gmp,VERSION);
  fclose(found);

// kmax=kmax_gmp and kmin=kmin_gmp

  mpz_t q_gmp;
  mpz_t r_gmp;
  mpz_init(q_gmp);
  mpz_init(r_gmp);

  signed long long int q,r,d=1<<30;

  mpz_fdiv_qr_ui(q_gmp,r_gmp,kmin_gmp,d);
  q=mpz_get_ui(q_gmp);
  r=mpz_get_ui(r_gmp);
  kmin=d*q+r;

  mpz_fdiv_qr_ui(q_gmp,r_gmp,kmax_gmp,d);
  q=mpz_get_ui(q_gmp);
  r=mpz_get_ui(r_gmp);
  kmax=d*q+r;

  unsigned long int j;
  register unsigned short J;
  unsigned long int i,num_primes=0;
  unsigned short table[16*3433];
  unsigned char b[32000];

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

// sieving up to 32000
// there are 3432 primes up to 32000

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

// positions in table: prm    0
//                     sp     1
//                     z      2
//                     first rem 3
//                     rem's diff  4,5,6,7,8,9,10,11,12,13,14
//                     start  15

// create rem table

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

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

  for(i=2;i<=3432;i++)
  {
  q=table[16*i],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[16*i+3]=h;
  table[16*i+4]=q-h;
  table[16*i+5]=2*table[16*i+3];
  if(table[16*i+5]>=q) table[16*i+5]-=q;
  table[16*i+6]=q-table[16*i+5];
  table[16*i+7]=2*table[16*i+5];
  if(table[16*i+7]>=q) table[16*i+7]-=q;
  table[16*i+8]=q-table[16*i+7];

  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[16*i+9]=h;
  table[16*i+10]=q-h;
  if(table[16*i+9]&1) table[16*i+11]=(table[16*i+9]+q)/2;
  else table[16*i+11]=table[16*i+9]/2;
  table[16*i+12]=q-table[16*i+11];
  if(table[16*i+11]&1) table[16*i+13]=(table[16*i+11]+q)/2;
  else table[16*i+13]=table[16*i+11]/2;
  table[16*i+14]=q-table[16*i+13];
  }

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

  unsigned long int F;
  signed long int v;

  unsigned long long int step=210,length=kmax-kmin;
  signed char x=4;  // we will use at least the first 4 primes
  unsigned long long int magic_constant=32000;  // to reduce the number of modular multiplications

  // we will use the first x primes

  while(length/(((unsigned long long) table[16*(x+1)])*step)>magic_constant)  x++,step*=table[16*x];

// redefine the starting and ending points

  kmin=kmin/step*step;
  kmax=kmax/step*step+step;

  unsigned long long int k1=kmin/step,k2=kmax/step;
  unsigned long long int small_length=k2-k1;  // One block length=32000 ( we are sieving also up to 32000 )
  unsigned long int bound=(k2-k1)%32000+1;
  unsigned long int limit,limit_numprime;
  unsigned short s_num=0;  // number of primes up to bound

  while((s_num+1<=3432)&(table[16*s_num+16]<=bound)) s_num++;

  mpz_t k1_gmp;
  mpz_init(k1_gmp);

  mpz_set_ui(k1_gmp,k1>>30);
  mpz_mul_2exp(k1_gmp,k1_gmp,30);
  mpz_add_ui(k1_gmp,k1_gmp,k1&1073741823);

// small rem table=srem ( for small primes )

  signed char srem[x+1][table[16*x]+1];

  for(i=2;i<=x;i++)
  {
  for(j=0;j<table[16*i];j++) srem[i][j]=1;
  for(j=1;j<=12;j++) srem[i][table[16*i+j+2]]=0;
  }

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

  mpz_t step_gmp;  // = step
  mpz_t u;
  mpz_init(step_gmp);
  mpz_init(u);

  mpz_set_ui(step_gmp,step>>30);
  mpz_mul_2exp(step_gmp,step_gmp,30);
  mpz_add_ui(step_gmp,step_gmp,step&1073741823);

// computation of z[i]=modinverse(step,prm[i]) and
//               sp[i]=k1%prm[i]

  for(i=x+1;i<=3432;i++)
  {
  mpz_set_ui(u,table[16*i]);
  mpz_invert(u,step_gmp,u);
  table[16*i+2]=mpz_get_ui(u);
  table[16*i+1]=k1%table[16*i];
  }

  signed long long int g,start0;
  unsigned char test;
  signed long int start,s0;
  register signed long int prm;

  mpz_t k;      // this is the multiplier
  mpz_t three;  // this is 3
  mpz_t g_gmp;  // g_gmp is equal to g
  mpz_t pr;
  mpz_t num_prp_gmp;
  mpz_init(k);
  mpz_init(g_gmp);
  mpz_init(three);
  mpz_set_ui(three,3);
  mpz_init(pr);
  mpz_init(num_prp_gmp);

  mpz_set_ui(num_prp_gmp,0);

  signed char f;
  signed long int z;
  unsigned long long int R;

  mpz_t R_gmp;
  mpz_init(R_gmp);

  signed long int A;

  register unsigned short A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11;

  signed long long int good_g=105,diff_g; // define the first good g as 0
  signed long int diff;

// computing the remainders for g=105

  for(i=x+1;i<=3432;i++)
  {
     v=16*i;
     prm=table[v];
     z=table[v+2];
     for(j=0;j<=11;j++)
     {
         A=(((signed long long) z) * ((signed long long) table[v+3+j]-105))%prm;
	 if(A<0)  table[v+3+j]=prm+A;
	     else table[v+3+j]=A;
     }
  }

// sorting the remainders in the table, and then save rem's difference table[j]=table[j]-table[16*i+3] for j>v

  for(i=1;i<=3432;i++)
  {
     v=16*i+3;
     for(j=v;j<=v+11;j++)
     {
        for(h=v;h<=v+11;h++)
	     if(table[j]<table[h]) F=table[j],table[j]=table[h],table[h]=F;
      }

      for(j=v+1;j<=v+11;j++)  table[j]-=table[v];
  }

  signed long int sp;
  signed char H;
  unsigned long int v14,v3;
  signed long int M;

  printf("Using the first %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);

  seconds = time (NULL);  // start the time

  // this is the main sieve + prp

  for(g=105;g<step;g+=210)
  {
  test=1,f=5;

  while(test&(f<=x)) test=srem[f][g%table[16*f]],f++;

  if(test)
  {
      diff_g=good_g-g;
      good_g=g;

      mpz_set_ui(g_gmp,g>>30);
      mpz_mul_2exp(g_gmp,g_gmp,30);
      mpz_add_ui(g_gmp,g_gmp,g&1073741823);

      v=16*x;

      for(i=x+1;i<=3432;i++)
      {

         v+=16;
         prm=table[v];

	   diff=(((signed long long) table[v+2])*diff_g)%prm;
	   if(diff<0) diff+=prm;

         table[v+3]+=diff;
         if(table[v+3]>=prm) table[v+3]-=prm;
      }

      mpz_set_ui(R_gmp,0);

      for(R=32000;R-32000<=small_length;R+=32000)
      {

      if(R==32000) {
            for(i=bound;i<=32000;i++) b[i]=1;
      }

      if(R>=small_length) limit=bound,limit_numprime=s_num;
               else       limit=32000,limit_numprime=3432;

      mpz_add_ui(R_gmp,R_gmp,32000);

      v=16*x;

      for(i=x+1;i<=limit_numprime;i++)
      {
          v+=16;
          prm=table[v];

          if(R==32000) 
              { 
                    if(table[v+3]>=table[v+1]) start=table[v+3]-table[v+1];
                                    else       start=prm+table[v+3]-table[v+1];
              }
          else  start=table[v+15];

          v3=v+3,z=v14=v+14,s0=prm-start;
	         while((z>v3)&(table[z]>=s0))  b[table[z]-s0]=0,z--;

          A1=table[v+4];
          A2=table[v+5];
          A3=table[v+6];
          A4=table[v+7];
          A5=table[v+8];
          A6=table[v+9];
          A7=table[v+10];
          A8=table[v+11];
          A9=table[v+12];
          A10=table[v+13];
          A11=table[v+14];

             for(J=start;J<limit-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]=0,M=limit-J,table[v+15]=prm-M;

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

      for(i=limit_numprime+1;i<=3432;i++)
      {
          v+=16;
          prm=table[v];

          if(R==32000) 
              { 
                    if(table[v+3]>=table[v+1]) start=table[v+3]-table[v+1];
                                    else       start=prm+table[v+3]-table[v+1];
              }
          else  start=table[v+15];

          v3=v+3,z=v14=v+14,s0=prm-start;
	         while((z>v3)&(table[z]>=s0))  b[table[z]-s0]=0,z--;

        if(start<limit) b[start]=0;
	  while((v3<v14)&(table[v3+1]+start<limit))  v3++,b[table[v3]+start]=0;
      }

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

      mpz_add(k,k1_gmp,R_gmp);
      mpz_sub_ui(k,k,32000-i);
      mpz_mul(k,k,step_gmp);
      mpz_add(k,k,g_gmp);


         mpz_set_ui(pr,1);
         mpz_mul_2exp(pr,pr,n);
         mpz_add(pr,pr,k);
         mpz_powm(q_gmp,three,pr,pr);
         if(mpz_cmp_ui(q_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);

      if(test)
      {
         mpz_mul_2exp(pr,pr,1);
         mpz_sub(pr,pr,k);
         mpz_powm(q_gmp,three,pr,pr);
	   if(mpz_cmp_ui(q_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }

      if(test)
      {
         mpz_mul_2exp(pr,pr,1);
         mpz_sub(pr,pr,k);
         mpz_powm(q_gmp,three,pr,pr);
	   if(mpz_cmp_ui(q_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }

      if(test)
      {
         mpz_set_ui(pr,1);
         mpz_mul_2exp(pr,pr,n);
         mpz_sub(pr,pr,k);
         mpz_powm(q_gmp,three,pr,pr);
         if(mpz_cmp_ui(q_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }

      if(test)
      {
         mpz_mul_2exp(pr,pr,1);
         mpz_add(pr,pr,k);
         mpz_powm(q_gmp,three,pr,pr);
         if(mpz_cmp_ui(q_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }
  
      if(test)
      {
         mpz_mul_2exp(pr,pr,1);
         mpz_add(pr,pr,k);
         mpz_powm(q_gmp,three,pr,pr);
         if(mpz_cmp_ui(q_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }

      if(test)
      {
         mpz_mul_2exp(pr,k,n);
         mpz_add_ui(pr,pr,1);
         mpz_powm(q_gmp,three,pr,pr);
         if(mpz_cmp_ui(q_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }
 
      if(test)
      {
         mpz_mul_2exp(pr,pr,1);
         mpz_sub_ui(pr,pr,1);
         mpz_powm(q_gmp,three,pr,pr);
         if(mpz_cmp_ui(q_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }

      if(test)
      {
         mpz_mul_2exp(pr,pr,1);
         mpz_sub_ui(pr,pr,1);
         mpz_powm(q_gmp,three,pr,pr);
         if(mpz_cmp_ui(q_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }

      if(test)
      {
         mpz_mul_2exp(pr,k,n);
         mpz_sub_ui(pr,pr,1);
         mpz_powm(q_gmp,three,pr,pr);
         if(mpz_cmp_ui(q_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }
 
      if(test)
      {
         mpz_mul_2exp(pr,pr,1);
         mpz_add_ui(pr,pr,1);
         mpz_powm(q_gmp,three,pr,pr);
         if(mpz_cmp_ui(q_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }

      if(test)
      {
         mpz_mul_2exp(pr,pr,1);
         mpz_add_ui(pr,pr,1);
         mpz_powm(q_gmp,three,pr,pr);
         if(mpz_cmp_ui(q_gmp,3)!=0) test=0;
         mpz_add_ui(num_prp_gmp,num_prp_gmp,1);
      }

	 if(test&(mpz_cmp(k,kmin_gmp)>=0)&(mpz_cmp(kmax_gmp,k)>=0))
	 {
	   printf("                                                                                                                             \r\r"),fflush(stdout);
	   gmp_printf("   %Zd  %d\n",k,n);

	   found=fopen("results_dodeca.txt","a+");
	   gmp_fprintf(found, "   %Zd  %d\n",k,n);
	   fclose(found);
	 }
      }
      b[i]=1;
      }
      }
      if((double) g/step-status/100.0>0.001) status=(double) 100*g/step,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_dodeca.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(kmin_gmp);
   mpz_clear(kmax_gmp);
   mpz_clear(q_gmp);
   mpz_clear(r_gmp);
   mpz_clear(k);
   mpz_clear(u);
   mpz_clear(g_gmp);
   mpz_clear(step_gmp);
   mpz_clear(k1_gmp);
   mpz_clear(R_gmp);
   mpz_clear(three);
   mpz_clear(pr);
   mpz_clear(num_prp_gmp);

   return 0;
}