﻿ compute measures
```/*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;
}
```