/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE MERTENS FUNCTION (sums of M(x/i)L(x), etc.)			      C
C  12/21/14 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "table2.h"
unsigned int liouvill(unsigned int a);
double numdiv(unsigned int a, unsigned int b);
double dirchar(unsigned int a, unsigned int p);
double mangoldt(unsigned int a, unsigned int *t);
unsigned int max=999; // maximum x value
unsigned int start=2;  // beginning x value, normally set to 2
unsigned int out=1;    // output flag, 0 for M(x), 1 for sum, 2 for histogram
unsigned int bins=100; // number of histogram bins (for out=2)
unsigned int n=1;      // decimation factor (normally set to 1)
unsigned int sflag=0;  // square root flag
unsigned int nunorm=0; // divide by x
unsigned int limit=0;  // partial sum flag (usually set to 0)
unsigned int delta=0;  // for Mangoldt function
unsigned int lflag=20;	// M(x/i)L(i), M(x/i)L(i)^2, M(x/i)^2L(i), M(x/i)^2L(i)^2,
			// L(i), L(i)^2, L(x) using Mertens
			// M(x/i)log(i)L(i), M(x/i)log(i)d(i)L(i), L(x),
			// M(x/i)sigma(i), |M(x/i)|sigma(i), sgn(M(x/i))sigma(i)
			// |sgn(M(x/i))|sigma(i), sigma(i), i, sigma(i) (not sum)
			// M(x/i)log(i)d(i) (#17)
			// M(x/i)sigma2(i)
			// M(x/i)char(i)
			// M(x/i)Mangoldt(i)
unsigned int fflag=1;  // select lambda or L for array (normally set to 1)
unsigned int norm=1;  // for lflag=12, 13 or 14. If set, divide by x
unsigned int nordel=0; // if set, subtract from sigsum (for norm=1)
double ratio=0.0; // if non-zero, output M(x) values (for nordel=1)
void main() {
short *m;
int *larry;
int sum,t,l,lsum;
unsigned int i,j,index,hisdif[101],width,flag,savsum,savind,numone,temp,count;
double delta,newsum,oldsum,maxdel,maxsav[10],temp1,sigsum;
FILE *Outfp;
Outfp = fopen("out7b.dat","w");
larry=(int*) malloc(500000);
if (larry==NULL)
   return;
m=(short*) malloc(500000);
if (m==NULL)
   return;
if (max>100000) {
   printf("x value too large \n");
   goto zskip;
   }
if ((out==2)&&(sflag==0)) {
   printf("parameter error \n");
   goto zskip;
   }
if ((out==6)&&(fflag!=0)) {
   printf("parameter error \n");
   goto zskip;
   }
if (lflag==6)
   goto askip;
//
// compute lambda or L
//
larry[0]=1;
lsum=1;
for (i=2; i<=max; i++) {
   l=liouvill(i);
   lsum=lsum+l;
   if (fflag!=0)
      larry[i-1]=lsum;
   else
      larry[i-1]=l;
//   printf(" %d \n",larry[i-1]);
   }
if (lflag==9)
   goto bskip;
//
// compute Mertens function
//
askip:
m[0]=1;
if (out==0)
   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)
      fprintf(Outfp," %d \n",1-sum);
   }
if (out==0)
   return;
//
// compute sum or histogram
//
if (out==2) {
   for (i=0; i<101; i++)
      hisdif[i]=0;
   for (i=0; i<10; i++)
      maxsav[i]=0.0;
   numone=0;
   oldsum=0.0;
   maxdel=0.0;
   savsum=0;
   savind=0;
   width=1;
   }
bskip:
count=0;
for (index=start; index<=max; index++) {
   newsum=0.0;
   sigsum=0.0;
   if (lflag==16) {
      if (norm==0)
	 fprintf(Outfp," %e\n",numdiv(index,4));
      else
	 fprintf(Outfp," %e\n",numdiv(index,4)/(double)index);
      continue;
      }
   for (i=n; i<=index; i+=n) {
      if ((index/i)<=limit)
	 break;
      if (nordel!=0)
	 sigsum=sigsum+numdiv(i,4);
      if (lflag==0) {
	 newsum=newsum+(double)m[index/i-1]*(double)larry[i-1];
	 }
      else {
	 if (lflag==7)
	    newsum=newsum+(double)m[index/i-1]*(double)larry[i-1]*log((double)i);
	 else {
	    if (lflag==6) {
	       temp=(unsigned int)(sqrt((double)i)+0.001);
	       if ((temp*temp)==i) {
		  newsum=newsum+(double)m[index/i-1];
		  }
	       }
	    else {
	       if (lflag==8)
		  newsum=newsum+(double)m[index/i-1]*(double)larry[i-1]*log((double)i)*numdiv(i,1);
	       else {
		  if (lflag==2)
		     newsum=newsum+(double)larry[i-1]*(double)m[index/i-1]*(double)m[index/i-1];
		  else {
		     if (lflag==5)
			newsum=newsum+(double)larry[i-1]*(double)larry[i-1];
		     else {
			if (lflag==4)
			   newsum=newsum+(double)larry[i-1];
			else {
			   if (lflag==3)
			      newsum=newsum+(double)larry[i-1]*(double)larry[i-1]*(double)m[index/i-1]*(double)m[index/i-1];
			   else {
			      if (lflag==1)
				 newsum=newsum+(double)larry[i-1]*(double)larry[i-1]*(double)m[index/i-1];
			      else {
				 if (lflag==9)
				    newsum=(double)larry[i-1];
				 else {
				    if (lflag==10) {
//				       printf(" %d %d %e \n",i,m[index/i-1],numdiv(i,4));
				       newsum=newsum+(double)m[index/i-1]*numdiv(i,4);
				       }
				    else {
				       if (lflag==11) {
					  temp1=(double)m[index/i-1];
					  if (temp1<0.0)
					     temp1=-temp1;
					  newsum=newsum+temp1*numdiv(i,4);
					  }
				       else {
					  if (lflag==12) {
					     temp1=(double)m[index/i-1];
					     if (temp1<0.0)
						newsum=newsum-numdiv(i,4);
					     if (temp1>0.0)
						newsum=newsum+numdiv(i,4);
					     }
					  else {
					     if (lflag==13) {
						if (m[index/i-1]!=0)
						   newsum=newsum+numdiv(i,4);
						}
					     else {
						if (lflag==14)
						   newsum=newsum+numdiv(i,4);
						else {
						   if (lflag==15)
						      newsum=newsum+i;
						   else {
						      if (lflag==17)
							 newsum=newsum+(double)m[index/i-1]*log((double)i)*numdiv(i,1);
						      if (lflag==18) {
//							 printf(" %d %d %e \n",i,m[index/i-1],numdiv(i,5));
							 newsum=newsum+(double)m[index/i-1]*numdiv(i,5);
							 }
						      if (lflag==19) {
//         					 printf(" %d %d %e \n",i,m[index/i-1],dirchar(i,101));
							 newsum=newsum+(double)m[index/i-1]*dirchar(i,101);
							 }
						      if (lflag==20) {
							 if (delta==0)
							    newsum=newsum+(double)m[index/i-1]*mangoldt(i,table);
							 else
							    newsum=newsum+((double)m[index/i-1]-0.10)*mangoldt(i,table);
							 }
						      }
						   }
						}
					     }
					  }
				       }
				    }
				 }
			      }
			   }
			}
		     }
		  }
	       }
	    }
	 }
      }
   if (((lflag==12)||(lflag==13)||(lflag==14))&&(norm!=0)) {
      newsum=newsum/(double)index;
      if (nordel!=0) {
	 sigsum=sigsum/(double)index;
	 newsum=sigsum-newsum;
	 }
      }
   if ((lflag!=1)&&(lflag!=3)&&(lflag!=5)&&(lflag!=10)&&(lflag!=11)&&(lflag!=12)&&(lflag!=13)&&(lflag!=14)&&(lflag!=17))
      newsum=-newsum;
   if (sflag!=0)
      newsum=sqrt(newsum);
   if (nunorm!=0)
      newsum=newsum/(double)index;
   if (nordel==0)
      printf(" %d %e \n",index,newsum);
   else
      printf(" %d %d %e %e\n",index,m[index-1],newsum,sigsum);
   if (out==1) {
      if (nordel==0)
	 fprintf(Outfp," %e\n",newsum);
      else {
	 if (ratio==0.0) {
	    if (newsum!=0.0)
	       fprintf(Outfp," %d %d %e %e %e\n",index,m[index-1],newsum,sigsum,sigsum/newsum);
	    }
	 else {
	    temp1=sigsum/newsum;
	    temp1=temp1-ratio;
	    if (temp1<0.0)
	       temp1=-temp1;
	    if ((temp1<=0.001)&&(newsum!=0.0)) {
	       fprintf(Outfp," %d \n",m[index-1]);
	       count=count+1;
	       }
	    else
	       fprintf(Outfp," %d \n",0);
	    }
	 }
      }
   else {
      delta=newsum-oldsum;
      flag=0;
      if (delta<0.0) {
	 delta=-delta;
	 flag=1;
	 }
      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);
	    goto zskip;
	    }
	 width=0;
	 j=(unsigned int)(delta*(double)bins+0.5);
	 if (j<101)
	    hisdif[j]=hisdif[j]+1;
	 }
      else
	 width=width+1;
      }
   }
printf("count=%d \n",count);
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:
free(larry);
free(m);
return;
}