/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  CHECK RESULTS GIVEN BY FAREY SERIES ALGORITHM                              C
C  06/11/07 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
double mertens2(unsigned int N, unsigned int *count, double *sum);
double mertens3(unsigned int N, unsigned int *count, double *sum);
double merten4b(unsigned int N, unsigned int *count, double *sum,
		unsigned int half);
//
unsigned int outflag=1; // delta or delta^2
unsigned int horder=0;	// half-order flag
//
void main() {
unsigned int N,count1[2],count2[2];
double sum1[2],sum2[4],suma,sumb,temp;
FILE *Outfp;
Outfp = fopen("out4f.dat","w");
//
// ORDER OF FAREY SERIES
//
suma=0.0;
/***************/
/* BEGIN LOOP  */
/***************/
for (N=3; N<=2560; N++) {
// suma=mertens3(N, count1, sum1);
   suma=0.0;
   sumb=merten4b(N, count2, sum2, horder);
   printf("N=%d, sums=%e %e %e %e \n",N,suma,sumb,sum2[2],sum2[3]);
   if (outflag==0)
      fprintf(Outfp," %e\n",sum2[2]);
   else
      fprintf(Outfp," %e\n",sum2[3]);
/* if ((sum1[0]!=sum2[0])||(sum1[1]!=sum2[1])||(suma!=sumb)) {
      temp=sum1[0]-sum2[0];
      if (temp<0.0)
	 temp=-temp;
      if (temp>0.000001) {
	 printf(" %d %e %e %e %e %e %e\n",N,suma,sumb,sum1[0],sum2[0],sum1[1],sum2[1]);
	 break;
	 }
      temp=sum1[1]-sum2[1];
      if (temp<0.0)
	 temp=-temp;
      if (temp>0.000001) {
	 printf(" %d %e %e %e %e %e %e\n",N,suma,sumb,sum1[0],sum2[0],sum1[1],sum2[1]);
	 break;
	 }
      temp=suma-sumb;
      if (temp<0.0)
	 temp=-temp;
      if (temp>0.000001) {
	 printf(" %d %e %e %e %e %e %e\n",N,suma,sumb,sum1[0],sum2[0],sum1[1],sum2[1]);
	 break;
	  }
      } 
   if ((count1[0]!=count2[0])||(count1[1]!=count2[1])) {
      printf(" %d %d %d %d \n",count1[0],count2[0],count1[1],count2[1]);
      break;
      } */
   }
fclose(Outfp);
/***************/
/* END LOOP    */
/***************/
return;
}