/*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 lagrange(unsigned int N, unsigned int *R);
//
unsigned int sqrtflg=1;  // compute sum of square roots of lengths
unsigned int flag=5;  // 0 (linear), 1 (quadratic), 2 (cubic), 3 (quartic)
//		      // 4 (sine), 5 (cosine)
void main() {
unsigned int N,S[250000],MAXN;
unsigned int L,I;
double temp1,temp2,func,len,tmpsum,expect,lensum,lenfunc,sqrtsum;
FILE *Outfp;
Outfp = fopen("out1c.dat","w");
//
// ORDER OF FAREY SERIES
//
MAXN=600;
/***************/
/* BEGIN LOOP  */
/***************/
for (N=2; N<=MAXN; N++) {
//
// GENERATE FAREY SERIES (USING HAROS' THEOREM AND LAGRANGE'S ALGORITHM)
//
   L=lagrange(N,S);
   len=0.0;
   tmpsum=0.0;
   lensum=0.0;
   sqrtsum=0.0;
   for (I=0; I<(L-1); I++) {
      temp1=((double)S[2*I]-(double)S[2*I+2])*((double)S[2*I]-(double)S[2*I+2]);
      temp1=temp1+((double)S[2*I+1]-(double)S[2*I+3])*((double)S[2*I+1]-(double)S[2*I+3]);
      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);
	 lenfunc=sin(temp2);
	 expect=cos(1.0);
	 }
      if (flag==5) {
	 func=cos(temp1);
	 lenfunc=cos(temp2);
	 expect=sin(1.0);
	 }
      tmpsum=tmpsum+func;
      lensum=lensum+lenfunc;
      }
   len=(len+1.0);
   temp2=len-(unsigned int)len;
   len=len/(double)L;
   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;
//
   if (flag==5) 	     // cos(0)=1
      tmpsum=tmpsum+1.0;
   tmpsum=tmpsum/(double)L;
   sqrtsum=(sqrtsum+1.0)/(double)L;
   printf("%d %d %e %e %e %e %e \n",N,L,len,sqrtsum,lensum,tmpsum,expect);
   if (sqrtflg==0)
      fprintf(Outfp," %e\n",len);
   else
      fprintf(Outfp," %e\n",sqrtsum);
   }
/***************/
/* END LOOP    */
/***************/
fclose(Outfp);
return;
}