/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  CHECK RESULTS GIVEN BY FAREY SERIES ALGORITHM                              C
C  (compute Franel measures with full and half orders [alpha and gamma])      C
C  (compute Franel measure with full order [beta])			      C
C  08/11/14 (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);
double merten4d(unsigned int N, unsigned int *count, double *sum,
		unsigned int half, unsigned int flag);
//
unsigned int outflag=1; // sum(|delta|), sqrt(A*sum(delsq^2)), sqrt(A+A*sum(delsq^2)),
			// when horder=0 uses total # fractions, old measure for
			// outflag=0 or 1, "beta" measure for outflag=2
unsigned int horder=0;	// half-order flag, "gamma" measure
unsigned int scale=0;	// scale when outflag=0 and horder=1
unsigned int scalep=0;	// scale when outflag=1 and horder=0
unsigned int scaleb=0;	// scale when outflag=2 and horder=0
unsigned int flag=1;	// for outflag=1 and horder=0
//
void main() {
unsigned int N,count1[2],count2[3];
double sum1[4],sum2[5],suma,sumb,temp,pi;
FILE *Outfp;
Outfp = fopen("out4h.dat","w");
//
// ORDER OF FAREY SERIES
//
pi=3.141592654;
sum1[0]=0.0;
sum1[1]=0.0;
sum1[2]=0.0;
sum1[3]=0.0;
count1[0]=0;
count1[1]=0;
suma=0.0;
/***************/
/* BEGIN LOOP  */
/***************/
for (N=2; N<=2500; N++) {    // MAXN=5128
// suma=merten4b(N, count1, sum1, horder);
   suma=0.0;
   sumb=merten4d(N, count2, sum2, horder, flag);
/* if (sum1[2]!=sum2[2]) {
      temp=sum1[2]-sum2[2];
      if (temp<0.0)
	 temp=-temp;
      if (temp>0.000001) {
	 printf("error 1: N=%d, sums=%e %e \n",N,sum1[2],sum2[2]);
	 if (N!=2)
	    break;
	 }
      }
   if (sum1[3]!=sum2[3]) {
      temp=sum1[2]-sum2[2];
      if (temp<0.0)
	 temp=-temp;
      if (temp>0.000001) {
	 printf("error 2: N=%d, sums=%e %e \n",N,sum1[3],sum2[3]);
	 if (N!=2)
	    break;
	 }
      } */
   printf("N=%d, sums=%e %e %e %e \n",N,suma,sumb,sum2[2],sum2[3]);
   if (outflag==0) {
      if ((horder==1)&&(scale==1)) {
	 temp=2.0*sqrt(2.0*sum2[2]);
	 fprintf(Outfp," %e\n",temp);
	 }
      else
	 fprintf(Outfp," %e\n",sum2[2]);
      }
   else {
      if (outflag==1) {
	 if (scalep==0)
	    fprintf(Outfp," %e\n",sum2[3]);
	 else
	    fprintf(Outfp," %e\n",2*pi*sum2[3]);
	 if ((horder==0)&&(flag==0)) {
	    temp=sqrt((double)N);
	    temp=0.4542*temp-0.4412;
	    if (sum2[3]<(temp-0.50)) {
	       printf("error 1: N=%d, temp=%e, sum=%e \n",N,temp-0.50,sum2[3]);
	       break;
	       }
	    if (sum2[3]>(temp+0.50)) {
	       printf("error 2: N=%d, temp=%e, sum=%e \n",N,temp+0.50,sum2[3]);
//	       break;
	       }
	    }
	 }
      else {
	 if (outflag==2) {
	    if (scaleb==0)
	       fprintf(Outfp," %e\n",sum2[4]);
	    else
	       fprintf(Outfp," %e\n",sum2[4]/sqrt(count2[2])-exp(1.0/(double)(3*N)));
	    }
	 else
	    fprintf(Outfp," %d\n",count2[2]);
	 }
      }
/* 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;
}