/*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>
double numdiv(unsigned int a, unsigned int flag);
double interp(int N, int D, int N1, int D1, unsigned int I, double init,
	      double *out);
unsigned int MAXN=76;		// maximum x
unsigned int delta=0;		// if set, add offset (slope for out=0)
unsigned int I=4;		// n value
unsigned int out=0;		// 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)
//
int frac[40]={0,1,1,2,1,4,1,3,1,6,2,5,2,15,31,105,29,140,19,42,41,420,
   76,385,201,1540,751,1430,1109,4004,803,2718,857,13411,3577,11807,
   721,17163,738,2897};
double slope[8]={.1653,.1544,.1399,.1278,.1176,.1093,.1022,.09607};
void main() {
double output[10000],init,sum11,temp;
unsigned int N,h,i,J;
short M[60001];
int t,sum;
FILE *Outfp;
Outfp = fopen("out5.dat","w");
if ((19*I)>10000) {
   printf(" array not big enough \n");
   goto zskip;
   }
if (MAXN>10000) {
   printf(" MAXN too big \n");
   goto zskip;
   }
if (MAXN>(19*I)) {
   printf(" MAXN too big \n");
   goto zskip;
   }
//
// compute Mertens function
//
if (out==6) {
   if (MAXN>60000) {
      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<10000; i++)
   output[i]=0.0;
J=(I+1)/2;
init=0.0;
for (i=0; i<19; i++)
   init=interp(frac[2*i],frac[2*i+1],frac[2*i+2],frac[2*i+3],I,init,
		&output[i*I]);
//for (i=0; i<19*I; i++) {
//printf(" %e \n",output[i]);
// fprintf(Outfp," %e \n",output[i]);
// }
if (I<9)
   temp=slope[I-1];
else
   temp=0.0;
for (h=2; h<=MAXN; h++) {
   sum11=0.0;
   for (i=1; i<=h; i++) {
      if (out==0)
	 sum11=sum11+output[h/i-1];
      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)
	 sum11=sum11+output[h/i-1]*numdiv(i,4);
      if (out==4)
	 sum11=sum11+output[h/i-1]*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
	    sum11=sum11+(output[h/i-1]-temp)*(double)M[i-1];
	 }
      }
      printf(" %d %e \n",h,sum11);
      fprintf(Outfp," %e\n",sum11);
   }
zskip:
fclose(Outfp);
return;
}