/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  GENERATE FAREY SERIES AND COMPUTE MEASURES (Franel and Landau)	      C
C  04/20/14 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
unsigned int divn_32(unsigned int *a, unsigned int *b, unsigned int d,
		  unsigned int m);
unsigned int lagrange(unsigned int N, unsigned int *R);
void setn(unsigned int *a, unsigned int b, unsigned int m);
void addn(unsigned int *a, unsigned int *b, unsigned int m);
void subn(unsigned int *a, unsigned int *b, unsigned int m);
//
unsigned int nocheck=1;   // set to 0 to use multiple-word arithmetic
unsigned int outflag=2;   // set to 0 to output all results, 1 to output
			  // |delta|, or 2 to output (A*sum(delta^2)^(1/2)
//
void main() {
unsigned int N,R[250000],MAXN;
unsigned int L,I,m;
double temp,temp1,delsq;
unsigned int A[32],B[32],C[32],T[32],U[32],Z[32];
FILE *Outfp;
Outfp = fopen("out1b.dat","w");
m=4;
setn(Z,0,m);
//
// ORDER OF FAREY SERIES
//
MAXN=600;
for (N=2; N<=MAXN; N++) {
//
// GENERATE FAREY SERIES (USING HAROS' THEOREM AND LAGRANGE'S ALGORITHM)
//
   L=lagrange(N,R);	       // generate Farey series
   temp=0.0;		       // |delta| measure
   setn(C,0,m); 	       // |delta| measure
   delsq=0.0;		       // delta^2 measure
   for (I=0; I<L; I++) {
      if (I!=0) {
	 temp1=(double)I/(double)(L-1)-(double)R[2*I]/(double)R[2*I+1];
	 if (nocheck==0) {
	    setn(A, 0, m);
	    A[0]=R[2*I];
	    divn_32(A, B, R[2*I+1], m);
	    setn(T, 0, m);
	    T[0]=I;
	    divn_32(T, U, L-1, m);
	    subn(U, B, m);
	    }
	 }
      else {
	 temp1=0.0;
	 setn(B, 0, m);
	 }
      delsq=delsq+temp1*temp1;
      if (temp1<0.0)
	 temp1=-temp1;
      temp=temp+temp1;
      if (nocheck==0) {
	 if ((B[0]&0x80000000)!=0)
	    subn(Z,B,m);
	 addn(B, C, m);
	 }
      }
   delsq=delsq*(double)(L-1);
   delsq=sqrt(delsq);	// (A*sum(delta^2)^(1/2), pg. 267 of H. M. Edwards
			// "Riemann's Zeta Function"
   printf("%d, %d, s=%e %e %e \n",N,L,temp,delsq,(double)C[0]+
     (double)C[1]/65536.0/65536.0+(double)C[2]/65536.0/65536.0/65536.0/65536.0);
   if (outflag==0) {
      fprintf(Outfp,"%d, %d, s=%e %e %e \n",N,L,temp,delsq,(double)C[0]+
	(double)C[1]/65536.0/65536.0+(double)C[2]/65536.0/65536.0/65536.0/65536.0);
      }
   else {
      if (outflag==1)
	 fprintf(Outfp," %e\n",temp);
      else
	 fprintf(Outfp," %e\n",delsq);
      }
   }
fclose(Outfp);
return;
}