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

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

double VERSION = 6.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[])

{
// define parameters for the program
	
  unsigned long int sieve_limit=100000;
  unsigned long int number_of_primes=9592;  // number of primes up to sieve_limit
  unsigned long int block_length=100000;
  unsigned long int magic_constant=100000;

// magic_constant>=block_length>=sieve_limit

  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]);

  double status;

  mpq_t status_gmp;
  mpq_t r_gmp;
  mpq_t q100_gmp;
  mpq_init(status_gmp);
  mpq_init(r_gmp);
  mpq_init(q100_gmp);
  
  mpq_set_ui(q100_gmp,100,1);  // q100=100
  
  mpz_t power_gmp;
  mpz_init(power_gmp);
  mpz_set_ui(power_gmp,10);
  mpz_pow_ui(power_gmp,power_gmp,12);  // power=10^12

     if(strchr(argv[2],'T')) {
        gmp_sscanf(argv[2],"%ZdT\n",kmin_gmp);
        mpz_mul(kmin_gmp,kmin_gmp,power_gmp);
        }
   else {
        mpz_set_str(kmin_gmp,argv[2],10);
        mpz_fdiv_q(kmin_gmp,kmin_gmp,power_gmp);
        mpz_mul(kmin_gmp,kmin_gmp,power_gmp);
        }

    if(strchr(argv[3],'T')) {
        gmp_sscanf(argv[3],"%ZdT\n",kmax_gmp);
        mpz_mul(kmax_gmp,kmax_gmp,power_gmp);
        }
   else {
        mpz_set_str(kmax_gmp,argv[3],10);
        mpz_cdiv_q(kmax_gmp,kmax_gmp,power_gmp);
        mpz_mul(kmax_gmp,kmax_gmp,power_gmp);
        }

// checking parameters
// Note that if n<25 then there is no octoproth for n

  if(n<25)
  {
  printf("n should be at least 25\n");
  print_usage_and_exit ();
  }

  if(mpz_cmp(kmin_gmp,kmax_gmp)>=0)
  {
  printf("kmax should be larger than kmin\n");
  print_usage_and_exit ();
  }
  
  if(mpz_sgn(kmin_gmp)<0)
  {
  printf("kmin should be nonnegative\n");
  print_usage_and_exit ();
  }

// print n, kmin, kmax and version to the screen and to results_octo.txt file
      
  FILE *found;

  mpz_t kmin_T_gmp;
  mpz_t kmax_T_gmp;
  mpz_init(kmin_T_gmp);
  mpz_init(kmax_T_gmp);

  mpz_fdiv_q(kmin_T_gmp,kmin_gmp,power_gmp);
  mpz_fdiv_q(kmax_T_gmp,kmax_gmp,power_gmp);

  printf("You can also find the k n values in results_octo.txt file ( These are 3-probable primes )\n");
  gmp_printf("   n=%d, kmin=%ZdT, kmax=%ZdT, version=%.1f, here T=10^12\n",n,kmin_T_gmp,kmax_T_gmp,VERSION);
  printf("Starting the sieve...\n");

  found=fopen("results_octo.txt","a+");
  gmp_fprintf(found,"   n=%d, kmin=%ZdT, kmax=%ZdT, version=%.1f, here T=10^12\n Starting the sieve...\n",n,kmin_T_gmp,kmax_T_gmp,VERSION);
  fclose(found);

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

  register unsigned long int J;
  register signed long int A1,A2,A3,A4,A5,A6,A7;

  register signed long int prm;

  for(i=1;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[12*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
//                     sp     1
//                     z      2
//                     first rem 3
//                     rem's diff  4,5,6,7,8,9,10
//                     start  11

// 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<=number_of_primes;i++)
  {
  q=table[12*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;
  }

  v=12*i;
  table[v+3]=h;
  table[v+4]=q-table[v+3];
  table[v+5]=2*table[v+3];
  if(table[v+5]>=q) table[v+5]-=q;
  table[v+6]=q-table[v+5];

  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+7]=h;
  table[v+8]=q-table[v+7];
  if(h&1) table[v+9]=(q+table[v+7])/2;
  else table[v+9]=table[v+7]/2;
  table[v+10]=q-table[v+9];
  }

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

  unsigned char x=3;  // we will use at least the first 3 primes

  mpz_t step_gmp;
  mpz_t s_gmp;
  mpz_t length_gmp;  // length=kmax-kmin
  mpz_t k1_gmp;  // k1=kmin/step
  mpz_init(step_gmp);
  mpz_init(s_gmp);
  mpz_init(length_gmp);
  mpz_init(k1_gmp);
  mpz_set_ui(step_gmp,30);  // 2*3*5=30

  mpz_sub(length_gmp,kmax_gmp,kmin_gmp);  

  mpz_mul_ui(s_gmp,step_gmp,table[12*x+12]);
  mpz_mul_ui(s_gmp,s_gmp,magic_constant);
    
  while(mpz_cmp(length_gmp,s_gmp)>0)
       {
       x++;
       mpz_mul_ui(step_gmp,step_gmp,table[12*x]);
       mpz_mul_ui(s_gmp,s_gmp,table[12*x+12]);
       }

// redefine the starting and ending points

  mpz_fdiv_q(kmin_gmp,kmin_gmp,step_gmp);
  mpz_mul(kmin_gmp,kmin_gmp,step_gmp);
  mpz_divexact(k1_gmp,kmin_gmp,step_gmp);
  
  mpz_cdiv_q(kmax_gmp,kmax_gmp,step_gmp);
  mpz_mul(kmax_gmp,kmax_gmp,step_gmp);

  mpz_t small_length_gmp;
  mpz_init(small_length_gmp);
  mpz_sub(small_length_gmp,kmax_gmp,kmin_gmp);
  mpz_divexact(small_length_gmp,small_length_gmp,step_gmp);
  unsigned long int small_length=mpz_get_ui(small_length_gmp);  // small_length=(kmax-kmin)/step
  
  unsigned long int bound=mpz_mod_ui(s_gmp,small_length_gmp,block_length);
  unsigned long int limit,limit_numprime;
  unsigned short small_numprime=0;  // number of primes up to bound

  while((small_numprime+1<=number_of_primes)&(table[12*small_numprime+12]<=bound)) small_numprime++;
  if(small_numprime<x+1) small_numprime=x+1;
  
// small rem table for small primes

  unsigned char small_rem[x+1][table[12*x]+1];

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

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

// computation of table[12*i+2]=modinverse(step_gmp,table[12*i]) and
//                table[12*i+1]=kmin_gmp modulo table[12*i]
// Here table[12*i] is the i-th prime number

  for(j=x+1;j<=number_of_primes;j++)
  {
  mpz_set_ui(s_gmp,table[12*j]);
  mpz_invert(s_gmp,step_gmp,s_gmp);
  table[12*j+2]=mpz_get_ui(s_gmp);
  table[12*j+1]=mpz_mod_ui(s_gmp,k1_gmp,table[12*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_set_ui(three_gmp,3);
  mpz_init(pr_gmp);
  mpz_init(num_prp_gmp);

  mpz_set_ui(num_prp_gmp,0);

  signed char f;
  signed long int z,A;
  unsigned long int F,R;

  mpz_t good_g_gmp;
  mpz_t diff_g_gmp;
  mpz_init(good_g_gmp);
  mpz_init(diff_g_gmp);
  mpz_set_ui(good_g_gmp,15);  // define the first good g as 15
  mpz_set_ui(g_gmp,15);
    
// computing the remainders for g=15

  for(i=x+1;i<=number_of_primes;i++)
  {
     v=12*i;
     prm=table[v];
     z=table[v+2];
     for(j=0;j<=7;j++)
     {
         A=(((signed long long) z) * ((signed long long) table[v+3+j]-15))%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[12*i+3] for j>v

  for(i=1;i<=number_of_primes;i++)
  {
     v=12*i+3;
     for(j=v;j<=v+7;j++)
     {
        for(h=v;h<=v+7;h++)
	     if(table[j]<table[h]) F=table[j],table[j]=table[h],table[h]=F;
      }

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

  unsigned long int v10,v3;
  signed long int M,start;

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

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

  unsigned long int s;
  signed long int diff,diff_g;

  signed char res[x];
  unsigned char sres[x];
  unsigned char sprm[x];
  
  for(i=4;i<=x;i++)  sprm[i]=table[12*i];  //  sprm[i]=prm(i)
  for(i=4;i<=x;i++)  sres[i]=30%table[12*i];  //  sres[i]=30%prm(i)
  for(i=4;i<=x;i++)  res[i]=table[12*i]-15%table[12*i];

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

  // this is the main sieve + prp

  while(mpz_cmp(step_gmp,g_gmp)>0)
  
  {
  test=1,f=4;

  for(i=4;i<=x;i++){
  res[i]+=sres[i];
  if(res[i]>=sprm[i]) res[i]-=sprm[i];
  }

  while(test&(f<=x))  test=small_rem[f][res[f]],f++;
  
  if(test)
  {
      mpz_sub(diff_g_gmp,g_gmp,good_g_gmp);
      mpz_set(good_g_gmp,g_gmp);
      v=12*x;

      if(mpz_cmp_ui(diff_g_gmp,1000000000)<0)  {
      diff_g=mpz_get_ui(diff_g_gmp);
	  for(i=x+1;i<=number_of_primes;i++)
      {
         v+=12;
         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;
      }
      }
      else {
           for(i=x+1;i<=number_of_primes;i++)
           {
	           v+=12;
	           prm=table[v];
	           mpz_mul_ui(s_gmp,diff_g_gmp,table[v+2]);
	           mpz_neg(s_gmp,s_gmp);
	           diff=mpz_mod_ui(s_gmp,s_gmp,prm);
	           table[v+3]+=diff;
	           if(table[v+3]>=prm) table[v+3]-=prm;
           } 
           }          	      

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

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

      if(R>=small_length) limit=bound,limit_numprime=small_numprime;
               else       limit=block_length,limit_numprime=number_of_primes;

      v=12*x;

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

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

          if(R==block_length)
              { 
                    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+11];

          v3=v+3,z=v10=v+10,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];

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

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

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

          if(R==block_length) 
              { 
                    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+11];

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

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

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

      mpz_add_ui(k_gmp,k1_gmp,R);
      mpz_sub_ui(k_gmp,k_gmp,block_length-i);
      mpz_mul(k_gmp,k_gmp,step_gmp);
      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-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 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+1)-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)
      {  // 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&(mpz_cmp(k_gmp,kmin_gmp)>=0)&(mpz_cmp(kmax_gmp,k_gmp)>=0))
	 {
	   printf("                                                                                                                             \r\r"),fflush(stdout);
	   gmp_printf("   %Zd  %d\n",k_gmp,n);

	   found=fopen("results_octo.txt","a+");
	   gmp_fprintf(found, "   %Zd  %d\n",k_gmp,n);
	   fclose(found);
	 }
      }
      b[i]=1;
      }
      }
      mpq_set_num(r_gmp,g_gmp);
      mpq_set_den(r_gmp,step_gmp);
      mpq_mul(r_gmp,r_gmp,q100_gmp);
      mpq_sub(r_gmp,r_gmp,status_gmp);    
    
      
      if(mpq_cmp_ui(r_gmp,1,10)>0) 
      {
      mpq_set_num(status_gmp,g_gmp);
      mpq_set_den(status_gmp,step_gmp);
      mpq_mul(status_gmp,status_gmp,q100_gmp);
      status=mpq_get_d(status_gmp);
      gmp_printf("   Status: %.1f percentage of the project is complete. Time thusfar: %ld sec. \r\r",status,time(NULL)-seconds),fflush(stdout);
      }      

  }
  mpz_add_ui(g_gmp,g_gmp,30);
  }
 	   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_octo.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
 
   mpq_clear(q100_gmp);
   mpq_clear(r_gmp);
   mpq_clear(status_gmp);
   mpz_clear(power_gmp);
   mpz_clear(kmin_gmp);
   mpz_clear(kmax_gmp);
   mpz_clear(kmin_T_gmp);
   mpz_clear(kmax_T_gmp);
   mpz_clear(k1_gmp);
   mpz_clear(step_gmp);
   mpz_clear(small_length_gmp);
   mpz_clear(s_gmp);
   mpz_clear(k_gmp);
   mpz_clear(g_gmp);
   mpz_clear(good_g_gmp);
   mpz_clear(diff_g_gmp);
   mpz_clear(three_gmp);
   mpz_clear(pr_gmp);
   mpz_clear(num_prp_gmp);

   return 0;
}