/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  CHECK WHETHER LENGTHS, CUMULATIVE LENGTHS, AND SQRT(LENGTH) U.D. (MOD 1)   C
C  04/20/14 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
unsigned int haros1(unsigned int N, unsigned int M, unsigned int *R,
		    unsigned int H, unsigned int K, unsigned int HP,
		    unsigned int KP);
//
unsigned int sqrtflg=1;  // compute sum of square roots of lengths
unsigned int outflg=2;
unsigned int flag=1;  // 0 (linear), 1 (quadratic), 2 (cubic), 3 (quartic)
//		      // 4 (sine), 5 (cosine)
void main() {
unsigned int N,S[250000],MAXN;
unsigned int H,K,HP,KP;
unsigned int L,I;
double temp1,temp2,func,len,tmpsum,expect,lensum,lenfunc,sqrtsum;
FILE *Outfp;
Outfp = fopen("out1e.dat","w");
//
// ORDER OF FAREY SERIES
//
MAXN=900;
/***************/
/* BEGIN LOOP  */
/***************/
for (N=2; N<=MAXN; N++) {
//
// GENERATE FAREY SERIES
//
   L=haros1(N,0,S,0,1,1,N);
   len=0.0;
   tmpsum=0.0;	 // length or square root of length
   lensum=0.0;	 // cumulative length or square root of length
   sqrtsum=0.0;
   for (I=0; I<(L-1); I++) {
      H=S[I]>>16;
      K=S[I]&0xffff;
      HP=S[I+1]>>16;
      KP=S[I+1]&0xffff;
      temp1=((double)H-(double)HP)*((double)H-(double)HP);
      temp1=temp1+((double)K-(double)KP)*((double)K-(double)KP);
      len=len+sqrt(temp1);
//
// fractional part of length
//
      if (sqrtflg==0)
	 temp1=sqrt(temp1);
      else
	 temp1=sqrt(sqrt(temp1));
      sqrtsum=sqrtsum+temp1;
      temp1=temp1-(unsigned int)temp1;
      temp2=len-(unsigned int)len;
      if (flag==0) {
	 func=temp1;
	 lenfunc=temp2;
	 expect=0.5;
	 }
      if (flag==1) {
	 func=temp1*temp1;
	 lenfunc=temp2*temp2;
	 expect=0.3333333;
	 }
      if (flag==2) {
	 func=temp1*temp1*temp1;
	 lenfunc=temp2*temp2*temp2;
	 expect=0.25;
	 }
      if (flag==3) {
	 func=temp1*temp1*temp1*temp1;
	 lenfunc=temp2*temp2*temp2*temp2;
	 expect=0.20;
	 }
      if (flag==4) {
	 func=sin(temp1);
//	 printf(" %e \n",func);
	 lenfunc=sin(temp2);
	 expect=cos(1.0);
	 }
      if (flag==5) {
	 func=cos(temp1);
	 lenfunc=cos(temp2);
	 expect=sin(1.0);
	 }
      tmpsum=tmpsum+func;
//    printf(" %e \n",tmpsum);
      lensum=lensum+lenfunc;
      }
   len=(len+1.0);
   temp2=len-(unsigned int)len;
   len=len/(double)L;		      // average lengths
   if (flag==0)
      lenfunc=temp2;
   if (flag==1)
      lenfunc=temp2*temp2;
   if (flag==2)
      lenfunc=temp2*temp2*temp2;
   if (flag==3)
      lenfunc=temp2*temp2*temp2*temp2;
   if (flag==4)
      lenfunc=sin(temp2);
   if (flag==5)
      lenfunc=cos(temp2);
   lensum=lensum+lenfunc;
   lensum=lensum/(double)L;	      // average cumulative lengths
//
   if (flag==5) 		      // cos(0)=1
      tmpsum=tmpsum+1.0;
   tmpsum=tmpsum/(double)L;	      // average func(fractional part)
   sqrtsum=(sqrtsum+1.0)/(double)L;   // average square roots of lengths
   printf("%d %d %e %e %e %e %e \n",N,L,len,sqrtsum,lensum,tmpsum,expect);
   if (outflg==0)
      fprintf(Outfp," %e\n",len);
   else {
      if (outflg==1)
	 fprintf(Outfp," %e\n",sqrtsum);
      else
	 fprintf(Outfp," %e\n",tmpsum);
      }
   }
/***************/
/* END LOOP    */
/***************/
fclose(Outfp);
return;
}