﻿ compute lengths of sides of Farey series polygon
```/*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;
}
```