﻿ generate Farey series and compute measures
```/*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);
}
}
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;
}
```