﻿ compute lengths of sides of Farey series polygon
```/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  CHECK MULTIPLE-WORD ARITHMETIC, COMPUTE LENGTHS, FRACTIONAL PARTS	      C
C  06/11/07 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
void divn_32(unsigned int *a, unsigned int *b, unsigned int d, unsigned int n);
unsigned int sqrtn(unsigned int a, unsigned int *b, unsigned int n);
void addn(unsigned int *a, unsigned int *b, unsigned int n);
unsigned int lagrange(unsigned int N, unsigned int *R);
//
unsigned int outflg=0; // 0 for average lengths, 1 for average fractional parts
//
void main() {
unsigned int N,S[90000],MAXN;
unsigned int L,I,J,o;
unsigned int P[64],Q[32],T[32],U[32];
double temp,temp1,len;
unsigned int m=8;	// number of words
FILE *Outfp;
Outfp = fopen("out1a.dat","w");
//
// ORDER OF FAREY SERIES
//
MAXN=100;
/***************/
/* BEGIN LOOP  */
/***************/
for (N=2; N<=MAXN; N++) {
//
// GENERATE FAREY SERIES
//
L=lagrange(N,S);
len=0.0;
temp=0.0;
for (I=0; I<m; I++) {
T[I]=0;
U[I]=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);
o=sqrtn((unsigned int)temp1,P, m*2);
if (o==0) {
printf("error \n");
goto L400;
}
for (J=m; J<2*m; J++)
Q[J-m]=P[J];
//
// fractional part
//
temp1=sqrt(temp1);
temp=temp+(temp1-(unsigned int)temp1);
Q[0]=0;
}
len=(len+1.0)/(double)L;
for (I=0; I<m; I++)
Q[I]=0;
Q[0]=1;
divn_32(T,Q,L,m);
temp=temp/(double)L;
divn_32(U,U,L,m);
temp1=(double)U[0]+(double)U[1]/65536.0/65536.0;
printf("n=%d %d   l=%e %e   %e %e \n",N,L,len,(double)Q[0]+(double)Q[1]/65536.0/65536.0,temp,temp1);
if (outflg==0)
fprintf(Outfp," %e\n",(double)Q[0]+(double)Q[1]/65536.0/65536.0);
else
fprintf(Outfp," %e\n",temp1);
}
/***************/
/* END LOOP    */
/***************/
L400:
fclose(Outfp);
return;
}
```