// Program to calculate bounds on f_k that guarantee system of APs does not cover.
// See https://arxiv.org/abs/1811.03547, particularly Table 1.
//
// Portability: we need "int" to be at least 32-bit and "long long" to be at least 64-bit.


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

#define llong long long			// Need this to be at least 64-bits .
#define EPS 1e-15			// Error margin to guard against fp arithmetic errors.
#define Max ((llong)1e12)		// Bound on largest prime considered.

// The following is a list of k for which f_k will be printed.  
// Must be in increasing order, and last term should be > Max for safety.
llong disp[]={2,3,4,5,6,7,8,9,10,21,100,1000,10000,51000,100000,1000000,
              (llong)1e7,(llong)1e8,(llong)1e9,(llong)1e10,Max+1};


// Speed of program largely limited by the speed at which one can generate primes.
// The following is a modified Sieve of Eratothenes.
// Due to memory constraints, we sieve in blocks of length 30*M, stored in sv[].
// Each block is divided into shorter sub-blocks of length 30, each stored in one "int".
// Thus, if bit i of sv[j] is set then sp+30*j+i is *not* prime.
// When a prime is sieved, only multiples that are relatively prime to 30=2*3*5 are sieved,
// thus speeding the process up by a factor 30/8. 
// Only primes in pr[] are sieved (initialized by primesetup()),
// thus the largest of these should be > sqrt(Max).

#define M ((llong)1000000)		// Sieve block size is 30*M

int sv[M];				// Sieve array, "int" must be at least 32-bit
llong sp;				// Sieve array base offset, multiple of 30*M, -30*M indicates empty
llong pr[M];				// List of small primes

// msk is a bit mask containing those numbers mod 30 relatively prime to 30
const msk=(1<<1)|(1<<7)|(1<<11)|(1<<13)|(1<<17)|(1<<19)|(1<<23)|(1<<29);

// add[n mod 30] is the amount needed to be added to n so that it is relatively prime to 30 
int add[]={1,6,5,4,3,2,1,4,3,2,1,2,1,4,3,2,1,2,1,4,3,2,1,6,5,4,3,2,1,2};

// Return first prime number >p (64-bit)
//
llong nextp(llong p){
 llong i,j,l,m,r;
 if(p<5){				// deal with small numbers
  if(p<2)return 2;
  if(p<3)return 3;
  return 5;
 }
 p+=add[p%30];				// next number relatively prime to 30
 while(1){
  r=p%30;m=(p-r-sp)/30;			// m = offset in array sv[], r = p mod 30
// sv[m]|((1<<r)-1) sets bits <r in sv[m] to 1 as we only are interested in bits >=r
// sv[m]==-1 indicates all bits set
  if(m>=0 && m<M && (sv[m]|((1<<r)-1))==-1){ // no primes in this sub-block
   m++; r=1;				// go to start of next sub-block
   while(m<M && sv[m]==-1)m++;		// skip sub-blocks that have no primes
  }
  if(m>=0&&m<M){			// found something in the current array
   while((sv[m]&(1<<r))!=0)r+=add[r];	// find it
   return 30*m+sp+r;			// return it
  } 
  p=30*m+sp+r;				// regenerate p
  sp=p-(p%(30*M));			// need to calculate a new block with offset sp
  for(i=0;i<M;i++)sv[i]=~msk;		// clear sieve array, but mark those with 2,3,5 as factors
  for(m=4;pr[m]*pr[m]<sp+(30*M);m++){	// pr[4]=7 is the first prime we need to consider
   i=sp%pr[m];i=(i!=0)?pr[m]-i:0;	// Want i = -sp mod pr[m] >=0, so i+sp first multiple of pr[]>=sp
   if(i+sp<pr[m]*pr[m])i=pr[m]*pr[m]-sp; // Only filter numbers >=m*m
   for(j=0;j<30;j++){			// loop of next 30 muliples of pr[m]
    if((msk&(1<<(i%30)))!=0){		// if these are relatively prime to 30
     for(l=i/30;l<M;l+=pr[m])sv[l]|=1<<(i%30); // mark those equivalent mod 30*pr[m]
    }
    i+=pr[m];
   }
  }
 }
}

// Sets up the list pr[] of small primes needed by nextp().
// Returns 1 if error.
//
int primesetup(){
 llong k,p;
// We first set up the list pr[] to contain odd numbers.
// Running nextp() will work as it tests for divisibility
// by numbers in pr[], and these only need to include the primes.
// (Additional numbers in pr[] just slow it down.)
 pr[0]=0; pr[1]=2;			// pr[0] not used, pr[1] = first prime
 for(k=2;k<M;k++)pr[k]=2*k-1;
 if(M<2||pr[M-1]*pr[M-1]<60*M)return 1;	// pr[] list might be too short
// Now we use nextp() itself to fill the array pr[] with genuine primes.
 sp=-30*M;				// mark the sieve array as empty
 p=0;
 for(k=1;k<M;k++)pr[k]=p=nextp(p);	// fill pr[] with primes
 if(pr[M-1]*pr[M-1]<Max+30*M)return 1;	// pr[] list too short
 return 0;
}


// The following assumes a given value of f_2 and then calculates f_k, k>2
// until either the condition f_k < (\log k+\log\log k-3)^2 k holds
// or we reach a stage where it is unlikely ever to hold.
// If successful, the f_k values can be used to bound g_k
// in Table 1 of https://arxiv.org/abs/1811.03547 
//
void try(double f2){
 llong i,p,k;
 double a,b,d,f;
 time_t start,now;
 printf("\nTry f_2=%.12f\n",f2);
 f=f2;
 i=0;					// index of next f_k to print in disp[]
 k=2;p=3;				// start at prime p_2=3
 time(&start);
 printf("            k           p_k                f_k   time");
 printf(" ln k+lnln k-(f_k/k)^.5\n");	// this number >3 yields success
 while(p<=Max){
  if(k==disp[i]){
   time(&now);
   printf("%13lld %13lld %18.12g %6d    ",k,p,f,(int)(now-start));
   printf("%f\n",log(k)+log(log(k))-sqrt(f/k));
   i++;
  }
  if(k%1000==0){			// we do this test rarely as it is slow
   a=log(k)+log(log(k))-3;
   a-=a*EPS;				// ensure we have an under-estimate
   if(f<a*a*k){
    printf("Succeeded at k=%lld, p=%lld\n",k,p);
    break;
   }
   if(f>(a+1)*(a+1)*k && k>1000){
    printf("Will prrobably fail for some k>%lld, p>%lld\n",k,p);
    break;
   }
  }
  p=nextp(p);
  k++;
  a=(3.*p-1)/((double)(p-1)*(p-1));
  b=f/((double)(p-1)*(p-1));		// this is "4b_k*f_{k-1}" in the paper
  if(b>=1){
   printf("Failed at k=%lld, p=%lld\n",k,p);
   break;
  }
  d=(1+a)/(1+sqrt(1+4*a*(1+a)/b));	// the optimal value of \delta_k
  f*=(1+EPS)*(1+a/(1-d))/(1-b/(4*d*(1-d)));
 }
 if(p>Max)printf("Not able to verify f_2 small enough\n");
 time(&now);
 printf("%d seconds\n",(int)(now-start));
}

main(){
 if(primesetup()){
  printf("Need to increase M or make M even\n");
  return;
 }
 try(1.260997);				// fast, but gives slightly worse bounds
 try(1.260997375);			// slower, gives the bounds listed in Table 1
}