/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE MERTENS FUNCTION (sum of sgn(M[x/i])) or sgn(M[x/i]))i)	      C
C  12/20/14 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <math.h>
#include <stdio.h>
unsigned int out=1;    // output flag, 0 for Mertens, 1 for sum, 2 for histogram
unsigned int max=100000; // maximum x value
unsigned int n=1;      // decimation factor (normally set to 1)
unsigned int bins=200; // number of histogram bins (for out=2)
unsigned int flag1=0;  // select plus or minus 1 or plus or minus i
unsigned int sflag=1;  // square root flag
unsigned int limit=1;  // partial sum flag (for flag1=1, normally set to 0)
//
void main() {
short m[500000];
int sum,t;
unsigned int i,j,index,hisdif[1001],width,savsum,savind,numone;
double delta,newsum,oldsum,maxdel,maxsav[10],temp;
FILE *Outfp;
Outfp = fopen("out7sy.dat","w");
if (max>500000) {
   printf("x too large \n");
   goto zskip;
   }
if ((out==2)&&(sflag==0)) {
   printf("parameter error \n");
   goto zskip;
   }
//
// compute Mertens function
//
m[0]=1;
if (out==0) {
// printf(" %d \n",1);
   fprintf(Outfp," %d \n",1);
   }
for (index=2; index<=max; index++) {
   sum=0;
   for (i=2; i<=(index/3); i++)
      sum=sum+m[index/i-1];
   sum=sum+(index+1)/2;
   t=1-sum;
   if ((t>32767)||(t<-32768)) {
      printf("overflow \n");
      goto zskip;
      }
   m[index-1]=(short)t;
   if (out==0) {
//    printf(" %d \n",1-sum);
      fprintf(Outfp," %d \n",1-sum);
      }
   }
if (out==0)
   return;
//
// compute sums or histograms
//
if (out==2) {
   for (i=0; i<1001; i++)
      hisdif[i]=0;
   for (i=0; i<10; i++)
      maxsav[i]=0.0;
   numone=0;
   if (flag1!=0)
      oldsum=1.0;
   else
      oldsum=0.0;
   maxdel=0.0;
   savsum=0;
   savind=0;
   width=1;
   }
for (index=2; index<=max; index++) {
   newsum=0.0;
   for (i=n; i<=index; i+=n) {
      if ((index/i)<=limit)
	 break;
      temp=(double)m[index/i-1];
      if (flag1!=0) {
	 if (temp>0.0)
	    newsum=newsum+(double)i;
	 if (temp<0.0)
	    newsum=newsum-(double)i;
	 }
      else {
	 if (temp>0.0)
	    newsum=newsum+1.0;
	 if (temp<0.0)
	    newsum=newsum-1.0;
//	 printf(" %d %d %d \n",i,index/i,m[index/i-1]);
	 }
      }
   if (limit!=0)
      newsum=-newsum;
   if (sflag!=0)
      newsum=sqrt(newsum);
   if (out==1) {
      printf(" %d %e\n",index,newsum);
      fprintf(Outfp," %e\n",newsum);
      }
   else {
      delta=newsum-oldsum;
      if (delta<0.0) {
	 delta=-delta;
	 if ((flag1!=0)&&(limit==0)) {
	    printf("negative step: index=%d \n",index);
	    goto zskip;
	    }
	 }
      if (delta>1.0)
	 numone=numone+1;
      if (delta>maxdel) {
	 for (j=9; j>0; j--)
	    maxsav[j]=maxsav[j-1];
	 maxsav[0]=delta;
	 maxdel=delta;
	 savsum=sum;
	 savind=index;
	 }
      oldsum=newsum;
      if (delta>0.00001) {
	 width=width+1;
	 if (width!=(width/n)*n) {
	    printf("error: width=%d \n",width);
	    if (index>n)
	       goto zskip;
	    }
	 width=0;
	 j=(unsigned int)(delta*(double)bins+0.5);
	 if (j<=bins)
	    hisdif[j]=hisdif[j]+1;
	 }
      else
	 width=width+1;
      }
   }
//
// output histogram
//
if (out==2) {
   printf("maximum difference=%e, sum=%d, index=%d \n",maxdel,savsum,savind);
   printf("number greater than one=%d \n",numone);
   for (i=0; i<10; i++)
      printf(" %e \n",maxsav[i]);
   j=0;
   for (i=0; i<=bins; i++) {
      fprintf(Outfp," %d \n",hisdif[i]);
      j=j+hisdif[i];
      }
   printf("sum=%d \n",j);
   }
zskip:
return;
}