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

#define N 21		// number of small primes
#define N0 6		// index of first prime (13) to process
#define S 0.75		// Coefficient in linear expression c_5(3)-S*c_5(1)
#define LINB 9.019 	// Bound on c_5(3)-S*c_5(1)
#define STEP 1e-4	// Step size on values of c_5(1)	
#define OPT 2		// number of coordinate-wise optimization passes
#define EPS 1e-15	// Error margin to guard against floating point inaccuracies
#define DIAG 0		// 2:Print diagnostics for each (c3,c1) pair
			// 1:Just print results for each (c3,c1) pair
			// 0:Only summary data

int p[N+1];		// list of small primes
double d[N+1];		// d[i]=\delta_i parameters
double mind,maxd;	// minimum and maximum delta values

// Setup list of small primes
// Very simple algorithm to generate first few primes
// We start the array at p[1]=2 to maintain consistency with paper
//
void setupprimes(){
 int i,j,n;
 p[1]=2;i=2;
 for(n=3;i<=N;n+=2){
  for(j=1;p[j]*p[j]<=n;j++)if((n%p[j])==0)break;
  if(p[j]*p[j]>n)p[i++]=n;
 }
 if(DIAG<1)return;
 printf("Primes:");
 for(i=1;i<N0;i++)printf(" %d",p[i]);	// these primes will not be used
 printf(";");
 for(i=N0;i<=N;i++)printf(" %d",p[i]);	// these are the relevant primes
 printf("\n");
} 

// Given value of c_5(3),c_5(1) and \delta_i=d[i], bound f_N
// To prevent accumulating floating point inaccuracies giving invalid bounds,
// we overestimate c3 and underestimate c1 and mu at each step.
//
double fn(double c3,double c1){
 double mu;
 int i;
 mu=1;
 for(i=N0;i<=N;i++){
  mu-=(c3-2*c1+1)/(4*d[i]*(1-d[i])*(p[i]-1)*(p[i]-1))+EPS*mu;
  if(mu<=0)return 1e99;			// failed, so return large number
  c1*=1+1/((p[i]-1)*(1-d[i]))-EPS;
  c3*=1+3/((p[i]-1)*(1-d[i]))+EPS;
 }
 return c3/mu;				// value of f_N
}

// Optimize value of d[i]=\delta_i, given c_5(3),c_5(1)
// We will assume f_N is convex function of d[i] and use golden section search on [a,b]
//
void optimize(double c3,double c1,int i,double a,double b){
 double x,y,fx,fy,g;
 g=(sqrt(5)-1)/2;
 x=b-(b-a)*g;
 d[i]=x;
 fx=fn(c3,c1);
 y=a+(b-a)*g;
 d[i]=y;
 fy=fn(c3,c1);
 while(fabs(x-y)>EPS){
  if(fx<fy){
   y=x-(y-x)*g;
   d[i]=y;
   fy=fn(c3,c1);
  }
  else {
   x=y+(y-x)*g;
   d[i]=x;
   fx=fn(c3,c1);
  }
 }
 d[i]=(x+y)/2;
 if(d[i]<mind)mind=d[i];
 if(d[i]>maxd)maxd=d[i];
} 

// Full optimization of all \delta_i's
//
double fullopt(double c3,double c1){
 int i,j;
 for(i=N0;i<=N;i++)d[i]=0.25;	// initial choice of delta's
 if(DIAG>1){
  printf("c3=%f c1=%f\n",c3,c1);
  printf("Initial f_%d=%f\n",N,fn(c3,c1));
  printf("New d[i]: ");
 }
 for(j=0;j<OPT;j++){
  for(i=N0;i<=N;i++){
   optimize(c3,c1,i,0.1,0.3);
  if(DIAG>1)printf("%f ",d[i]);
  }
  if(DIAG>1)printf("\nNew f_%d=%f\n",N,fn(c3,c1));
 }
 return fn(c3,c1);
} 


main(){
 double c1,c3,mx,f,cm;
 int i;
 setupprimes();
 mind=1;maxd=0;
 mx=0;
 for(c1=1.;c1<=5.;c1+=STEP){
  c3=LINB+(c1+STEP)*S;
  f=fullopt(c3,c1);
  if(f>mx){mx=f;cm=c1;}
  if(DIAG)printf("c3=%f c1=%f f_%d=%f\n",c3,c1,N,f);
 }
 printf("\nMaximum f_%d=%f\n",N,mx);
 printf("Bound for line c_%d(3)-%fc_%d(1)=%f\n",N0-1,S,N0-1,LINB);
 printf("Step size for c_%d(1) = %f\n",N0-1,STEP);
 printf("Worst case at dominating point (c_%d(3),c_%d(1))=(%f,%f)\n",N0-1,N0-1,cm,LINB+(cm+STEP)*S);
 printf("All delta_i values in [%f,%f]\n",mind,maxd);
 printf("Optimization performed via %d passes of coordinate-wise optimization\n",OPT);
}