﻿ compute variants of measures
/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  GENERATE FAREY SERIES						      C
C  (compute variants of Franel's measures)                                    C
C  04/20/14 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
unsigned int haros5(unsigned int N, unsigned int M, unsigned int *R,
unsigned int H, unsigned int K, unsigned int HP,
unsigned int KP, unsigned int D);
//
unsigned int half=1;  // half-order
unsigned int out=0;   // |delta| or sqrt(A*delta^2)
unsigned int old=0;   // old measure
unsigned int flag=0;  // scaling flag (for half=1, out=0, old=0)
//
void main() {
unsigned int N,R[250000],MAXN;
unsigned int L,I,total;
double temp,temp1,sum,sum1,sum2,sum3,temp2;
double pi=3.141592654;
FILE *Outfp;
Outfp = fopen("out1b3a.dat","w");
//
// ORDER OF FAREY SERIES
//
if ((half!=0)&&(old!=0)) {
printf("parameter error \n");
goto zskip;
}
MAXN=900;	  // MAXN=900
for (N=2; N<=MAXN; N++) {
//
// GENERATE FAREY SERIES
//
L=haros5(N,0,R,0,1,1,N,1);
total=L-1;
if (half!=0)
total=total/2;
sum=0.0;		      // |delta| measure
sum1=0.0;		      // delta^2 measure
sum2=0.0;		      // check sum of cosines
sum3=0.0;		      // check sum of sines
for (I=1; I<=total; I++) {
temp2=(double)I/(double)total;
if (temp2>=1.0)
temp2=temp2-1.0;
temp1=temp2-(double)(R[I]>>16)/(double)(R[I]&0xffff);
//    printf(" %d %d %e %e \n",R[I]>>16,R[I]&0xffff,cos(2*pi*temp2),sin(2*pi*temp2));
sum2=sum2+cos(2*pi*temp2);
sum3=sum3+sin(2*pi*temp2);
temp=temp1;
if (temp<0.0)
temp=-temp;
if ((I!=total)||(old==0)) {
sum=sum+temp;
sum1=sum1+temp1*temp1;
}
//    printf(" %d %d %e %e \n",R[I]>>16,R[I]&0xffff,sum,sum1);
}
if ((old!=0)&&(1.0<=(sqrt((double)N)*sum1))) {
printf("error: N=%d, sum1=%e \n",N,sum1);
goto zskip;
}
sum1=sum1*(double)total;
sum1=sqrt(sum1);
if (out==0) {
if (flag==0)
fprintf(Outfp," %e\n",sum);
else
fprintf(Outfp," %e\n",2.0*sqrt(2.0*sum));
}
else
fprintf(Outfp," %e\n",sum1);
printf(" %d %e %e %e %e \n",total,sum,sum1,sum2,sum3);
}
zskip:
fclose(Outfp);
return;
}