/*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];
      addn(Q,T,m);
//
// fractional part
//
      temp1=sqrt(temp1);
      temp=temp+(temp1-(unsigned int)temp1);
      Q[0]=0;
      addn(Q,U,m);
      }
   len=(len+1.0)/(double)L;
   for (I=0; I<m; I++)
      Q[I]=0;
   Q[0]=1;
   addn(Q,T,m);
   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;
}