/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  INTERPOLATE FRACTIONS AND COMPUTE FUNCTIONS				      C
C  05/17/15 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "table2.h"
#include "output2.h"  // if output5, scale=1000, MAXN=9999
		      // if output4, scale=1000, MAXN=3999
		      // if output2, scale=1000000, MAXN=999
		      // if output1, scale=1000000, MAXN=468 (computed on C64)
double numdiv(unsigned int a, unsigned int flag);
int liouvill(unsigned int a, unsigned int *t);
double mangoldt(unsigned int a, unsigned int *t);
double dirchar(unsigned int a, unsigned int p);
double interp1(double ND, double ND1, unsigned int I,
	       double init, double *out);
double scale=1000000.0;	// for output1.h or output2.h
//double scale=1000.0;	     // for output4.h or output5.h
unsigned int count=1000;	 // number of limits
unsigned int MAXN=999;   // maximum x
unsigned int delta=0;		// if set, add offset (slope for out=0)
unsigned int weight=0;      // multiply by h/i
unsigned int prime=101;     // prime for Legendre symbol
unsigned int I=1;	       // m value
unsigned int out=12;		// 0 for sum of delta(x/i)
			        // 1 for sum of delta(x/i)*log(i)
			        // 2 for sum of delta(x/i)*sigma0(i)
			        // 3 for sum of delta(x/i)*sigma1(i)
			        // 4 for sum of delta(x/i)*sigma2(i)
				// 5 for sum of delta(x/i)^2
				// 6 for sum of delta(x/i)*M(i)
				// 7 for sum of delta(x/i)*L(i)^2
				// 8 for sum of delta(x/i)^2*L(i)^2
				// 9 for sum of delta(x/i)^2*L(i)
				// 10 for sum of delta(x/i)*log(i)*d(i)
				// 11 for sum of delta(x/i)/i
				// 12 for sum of delta(x/i)*mangoldt(i)
				// 13 for sum delta(x/i), i perfect square
				// 14 for sum delta(x/i)*legendre(i)
				// 99 for interpolated limits
double slope[8]={.1704,.1548,.1402,.128,.1179,.1095,.1023,.09617};
//
void main() {
double init,sum11,temp,tsum;
double *output;
unsigned int N,h,i,J,tmp;
int t,sum;
short M[100001];
FILE *Outfp;
Outfp = fopen("out6.dat","w");
if ((count*I)>40000000) {
   printf(" array not big enough \n");
   goto zskip;
   }
if (MAXN>40000000) {
   printf(" MAXN too big \n");
   goto zskip;
   }
if (MAXN>(count*I)) {
   printf(" MAXN too big \n");
   goto zskip;
   }
if ((MAXN>199999)&&(out==12)) {
   printf(" MAXN too big \n");
   goto zskip;
   }
output=(double*) malloc(40000000);
if (output==NULL)
   return;
//
// compute Mertens function
//
if (out==6) {
   if (MAXN>100000) {
      printf("not enough memory allocated");
      goto zskip;
      }
   M[0]=1;
   for (N=2; N<=(MAXN+1); N++) {
      sum=0;
      for (i=2; i<=(N/3); i++)
	 sum=sum+M[N/i-1];
      sum=sum+(N+1)/2;
      t=1-sum;
      M[N-1]=(short)t;
//	  printf(" %d \n",M[N-1]);
      }
   }
for (i=0; i<100000; i++)
   output[i]=0.0;
if (I<9)
   temp=slope[I-1];
else
   temp=0.0;
J=(I+1)/2;
init=0.0;
for (i=0; i<count; i++)
   init=interp1((double)frac[i]/scale,(double)frac[i+1]/scale,I,init,
		&output[i*I]);
if (out==99) {
   for (i=0; i<(count*I); i++) {
      printf(" %e \n",output[i]);
      fprintf(Outfp," %e \n",output[i]);
      }
   goto zskip;
   }
for (h=2; h<=MAXN; h++) {
   sum11=0.0;
   tsum=0;
   for (i=1; i<=h; i++) {
      if (out==0) {
	 if (weight==0)
	    sum11=sum11+output[h/i-1];
	 else
	    sum11=sum11+output[h/i-1]*(double)(h/i);
	 }
      if (out==1) {
	 if (delta==0)
	    sum11=sum11+output[h/i-1]*log((double)i);
	 else
	    sum11=sum11+(output[h/i-1]+temp)*log((double)i);
	 }
      if (out==2) {
	 if (delta==0)
	    sum11=sum11+output[h/i-1]*2.0*numdiv(i,1);
	 else
	    sum11=sum11+(output[h/i-1]+temp)*2.0*numdiv(i,1);
	 }
      if (out==3) {
	 if (delta==0)
	    sum11=sum11+output[h/i-1]*numdiv(i,4);
	 else
	    sum11=sum11+(output[h/i-1]+temp)*numdiv(i,4);
	 }
      if (out==4) {
	 if (delta==0)
	    sum11=sum11+output[h/i-1]*numdiv(i,5);
	 else
	    sum11=sum11+(output[h/i-1]+temp)*numdiv(i,5);
	 }
      if (out==5) {
	 if (delta==0)
	    sum11=sum11+output[h/i-1]*output[h/i-1];
	 else
	    sum11=sum11+(output[h/i-1]+temp)*(output[h/i-1]+temp);
	 }
      if (out==6) {
	 if (delta==0)
	    sum11=sum11+output[h/i-1]*(double)M[i-1];
	 else {
	    if (weight==0)
	       sum11=sum11+(output[h/i-1]+temp)*(double)M[i-1];
	    else
	       sum11=sum11+(output[h/i-1]+temp)*(double)M[i-1]*(double)(h/i);
	    }
	 }
      if (out==7) {
	 t=liouvill(i, table);
	 tsum=tsum+t;
	 if (delta==0)
	    sum11=sum11+output[h/i-1]*(double)tsum*(double)tsum;
	 else
	    sum11=sum11+(output[h/i-1]+temp)*(double)tsum*(double)tsum;
	 }
      if (out==8) {
	 t=liouvill(i, table);
	 tsum=tsum+t;
	 if (delta==0)
	    sum11=sum11+output[h/i-1]*output[h/i-1]*(double)tsum*(double)tsum;
	 else
	    sum11=sum11+(output[h/i-1]+temp)*(output[h/i-1]+temp)*(double)tsum*(double)tsum;
	 }
      if (out==9) {
	 t=liouvill(i, table);
	 tsum=tsum+t;
	 if (delta==0)
	    sum11=sum11+output[h/i-1]*output[h/i-1]*(double)tsum;
	 else
	    sum11=sum11+(output[h/i-1]+temp)*(output[h/i-1]+temp)*(double)tsum;
	 }
      if (out==10) {
	 if (delta==0)
	    sum11=sum11+output[h/i-1]*log((double)i)*numdiv(i,1);
	 else
	    sum11=sum11+(output[h/i-1]+temp)*log((double)i)*numdiv(i,1);
	 }
      if (out==11) {
	 if (delta==0)
	    sum11=sum11+output[h/i-1]/(double)i;
	 else
	    sum11=sum11+(output[h/i-1]+temp)/(double)i;
	 }
      if (out==12) {
	 if (delta==0) {
	    if (weight==0)
	       sum11=sum11+output[h/i-1]*mangoldt(i,table);
	    else
	       sum11=sum11+output[h/i-1]*mangoldt(i,table)*(double)(h/i);
	    }
	 else
	    sum11=sum11+(output[h/i-1]+temp)*mangoldt(i,table);
	 }
      if (out==13) {
	 tmp=(unsigned int)(sqrt((double)i)+0.001);
	 if ((tmp*tmp)==i) {
	    if (delta==0)
	       sum11=sum11+output[h/i-1];
	    else {
	       if (weight==0)
		  sum11=sum11+output[h/i-1]+temp;
	       else
		  sum11=sum11+(output[h/i-1]+temp)*(double)(h/i);
	       }
	    }
	 }
      if (out==14) {
//	 printf(" %d %d %e \n",i,output[h/i-1],dirchar(i,prime));
	 if (delta==0)
	    sum11=sum11+output[h/i-1]*dirchar(i,prime);
	 else
	    sum11=sum11+(output[h/i-1]+temp)*dirchar(i,prime);
	 }
      }
   if (out!=5) {
      printf(" %d %e \n",h,-sum11);
      fprintf(Outfp," %e\n",-sum11);
      }
   else {
      printf(" %d %e \n",h,sum11);
      fprintf(Outfp," %e\n",sum11);
      }
   }
zskip:
free(output);
fclose(Outfp);
return;
}