/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE UPPER BOUND MINUS LOWER BOUND (AND ESTIMATED VALUES) 	      C
C  06/20/14 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
unsigned int euler(unsigned int a);
double mertens4(unsigned int, unsigned int *count, double *sum);
//
unsigned int out=15; // 0, sum of M([x/i]) when 4 divides phi(i)
		     // 1, estimated sum of M([x/i]) when 4 divides phi(i)
		     // 2, exceptions (when estimated signs not correct)
		     // 3, sum of a(x)*m(x/i)/m(x)
		     // 4, sum of b(x)*n(x/i)/n(x)
		     // 5, sum of i*a(x)*m(x/i)/m(x)
		     // 6, sum of i*b(x)*n(x/i)/n(x)
		     // 7, sum of a(x/i)
		     // 8, sum of b(x/i)
		     // 9, sum of i*a(x/i)
		     // 10, sum of i*b(x/i)
		     // 11, sum of m(x/i)
		     // 12, sum of n(x/i)
		     // 13, sum of m(x/i)-n(x/i)
		     // 14, sum of i*(m(x/i)-n(x/i))
		     // 15, 1/2 sum of i*M([x/i]) where M(1), M(3) are set to 0
		     // 16, Mertens function
		     // 17, differences in number of fractions before, after 1/4
		     //
double inc=0.0;      // esum+increment
unsigned int base=1;  // 0 or Mertens
unsigned int n=1;     // decimation factor
//
void main() {
unsigned int N,MAXN;  // MAXN=2560
unsigned int i,j,temp,clsave,crsave,count[2],histop[100],histopn[100];
unsigned int histon[100],histonn[100],histox,histoxn;
double sum,sum1,msum,mlsave,mrsave,lratio,rratio,tsum[2],esum,tmpsum;
double lsum,rsum,lsum2,rsum2,sum2,sum3,sum4,sum5,sum6,sum7,sum8,sum9;
double newsum,oldsum;
FILE *Outfp;
Outfp = fopen("out1ba.dat","w");
for (i=0; i<100; i++) {
   histop[i]=0;
   histopn[i]=0;
   histon[i]=0;
   histonn[i]=0;
   }
histox=0;
histoxn=0;
MAXN=900;
oldsum=0.0;
for (N=2; N<=MAXN; N++) {
   sum=0.0;
   sum1=0.0;
   sum2=0.0;
   sum3=0.0;
   sum4=0.0;
   sum5=0.0;
   sum6=0.0;
   sum7=0.0;
   sum8=0.0;
   sum9=0.0;
   for (i=1; i<=N; i++) {
      msum=mertens4(N/i, count, tsum);
      if ((out==16)&&(i==1))
	 fprintf(Outfp," %e\n",msum);
      if ((out==17)&&(i==1))
	 fprintf(Outfp," %d\n",(int)count[0]-(int)count[1]);
      if (i==1) {
	 mlsave=tsum[0];
	 mrsave=tsum[1];
	 clsave=count[0];
	 crsave=count[1];
	 if (base==0)
	    esum=msum;
	 else {
	    esum=0.0;
	    lsum=mlsave;
	    rsum=mrsave;
	    lsum2=mlsave;
	    rsum2=mrsave;
	    }
//	 printf(" %d %d %e %e \n",clsave,crsave,mlsave,mrsave);
	 }
//
// compute sums
//
      temp=euler(i);
      if ((temp/4)*4!=temp) {
	 sum=sum+msum*(double)i/(double)temp;
	 sum1=sum1+esum*(double)i/(double)temp;
	 }
      sum2=sum2+tsum[0];
      if ((N/i)==1)	     // solve weight problem
	 sum3=sum3+0.5;
      else {
	 if ((N/i)==3)
	    sum3=sum3-0.5;
	 else
	    sum3=sum3+tsum[1];
	 }
      sum4=sum4+tsum[0]*(double)i;
      sum5=sum5+tsum[1]*(double)i;
      sum6=sum6+(double)count[0];
      sum7=sum7+(double)count[1];
      sum8=sum8+(double)count[0]-(double)count[1];
      sum9=sum9+(double)i*((double)count[0]-(double)count[1]);
//    printf(" %d %d %d %e %e %e \n",i,count[0],N/i,sum4,sum5,2.0*(tsum[0]+tsum[1]));
//    printf(" %d %d %d %e %e \n",i,count[0],N/i,sum,msum);
      if (i!=1) {
	 if (clsave!=0)
	    lratio=(double)count[0]/(double)clsave;
	 else
	    lratio=1.0;
	 if (crsave!=0)
	    rratio=(double)count[1]/(double)crsave;
	 else
	    rratio=1.0;
	 esum=mlsave*lratio+mrsave*rratio;
	 esum=esum*2.0;
	 if ((i-1)==((i-1)/n)*n) {
	    lsum=lsum+mlsave*lratio;
	    rsum=rsum+mrsave*rratio;
	    }
	 lsum2=lsum2+(double)i*mlsave*lratio;
	 rsum2=rsum2+(double)i*mrsave*rratio;
//	 printf(" %e %e %e \n",lsum,rsum,lsum+rsum);
//	 fprintf(Outfp," i=%d %d %d %e %e %e %e %e %e \n",i,count[0],count[1],tsum[0],tsum[1],
//		  lratio,rratio,msum,esum);
//
// histogram (for msum>0)
//
	 if ((msum>0.0)&&(i<=(N/2))) {
	    j=(unsigned int)(msum+0.0001);
	    if ((j<100)&&(j!=0)) {
	       histop[j]=histop[j]+1;
	       if (esum+inc>0.0)
		  histopn[j]=histopn[j]+1;
	       }
	    printf("N=%d, i=%d, N/i=%d, msum=%e, esum=%e \n",N,i,N/i,msum,esum);
	    if (out==2)
	       fprintf(Outfp,"N=%d, i=%d, N/i=%d, msum=%e, esum=%e \n",N,i,N/i,msum,esum);
	    }
//
// histogram (for msum<0)
//
	 if (msum<0.0) {
	    tmpsum=-msum;
	    j=(unsigned int)(tmpsum+0.0001);
	    if ((j<100)&&(j!=0)) {
	       histon[j]=histon[j]+1;
	       if (esum-inc<0.0)
		  histonn[j]=histonn[j]+1;
	       }
	    }
//
// histogram (for msum==0)
//
	 if (msum<0.0)
	    tmpsum=-msum;
	 else
	    tmpsum=msum;
	 j=(unsigned int)(tmpsum+0.0001);
	 if ((j==0)&&(i<=(N/3))) {
	    histox=histox+1;
	    if (esum<0.0)
	       tmpsum=-esum;
	    else
	       tmpsum=esum;
	    j=(unsigned int)(tmpsum+0.1);
	    if (j==0)
	       histoxn=histoxn+1;
	    }
	 }
      }
   newsum=sum4+sum5;
   if (((int)(N-6)/12)*12==(int)(N-6)) {
      if (newsum<=oldsum) {
	 printf("error 1: N=%d, newsum=%e, oldsum=%e \n",N,newsum,oldsum);
	 goto zskip;
	 }
      }
   else {
      if (newsum>oldsum) {
	 printf("error 2: N=%d, newsum=%e, oldsum=%e \n",N,newsum,oldsum);
	 goto zskip;
	 }
      }
   oldsum=newsum;
   if (sum<0.0)
      sum=-sum;
   if (sum1<0.0)
      sum1=-sum1;
   if (out==0)
      fprintf(Outfp," %e\n",sum);
   else {
      if (out==1)
	 fprintf(Outfp," %e\n",sum1);
      else {
	 if (out==3)
	    fprintf(Outfp," %e\n",lsum);
	 if (out==4)
	    fprintf(Outfp," %e\n",rsum);
	 if (out==5)
	    fprintf(Outfp," %e\n",lsum2);
	 if (out==6)
	    fprintf(Outfp," %e\n",rsum2);
	 if (out==7)
	    fprintf(Outfp," %e\n",sum2);
	 if (out==8)
	    fprintf(Outfp," %e\n",sum3);
	 if (out==9)
	    fprintf(Outfp," %e\n",sum4);
	 if (out==10)
	    fprintf(Outfp," %e\n",sum5);
	 if (out==11)
	    fprintf(Outfp," %e\n",sum6);
	 if (out==12)
	    fprintf(Outfp," %e\n",sum7);
	 if (out==13)
	    fprintf(Outfp," %e\n",sum8);
	 if (out==14)
	    fprintf(Outfp," %e\n",sum9);
	 if (out==15)
	    fprintf(Outfp," %e\n",sum4+sum5);
	 }
      }
   }
if (out==2)
   fprintf(Outfp,"\n");
printf("  %d, %d, %d \n",0,histox,histoxn);
if (out==2)
   fprintf(Outfp,"  %d, %d, %d \n",0,histox,histoxn);
for (i=1; i<20; i++) {
   printf("  %d, %d %d %d %d \n",i,histop[i],histopn[i],histon[i],histonn[i]);
   if (out==2)
      fprintf(Outfp,"  %d, %d %d %d %d \n",i,histop[i],histopn[i],histon[i],histonn[i]);
   }
zskip:
fclose(Outfp);
return;
}