/*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;
}