﻿ compute Farey series using half-order series
```/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE MERTENS FUNCTION (using half-order Farey series)		      C
C  05/16/14 (DKC)							      C
C                                                                             C
C  This C program computes the sum of cos(2*pi*r_v) where r_v is a fraction   C
C  in the Farey series. 						      C
C                                                                             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 lagrange1(unsigned int N, unsigned int D, unsigned int O);
unsigned int out=0;   // select upper curve, lower curve, or sum
//
void main() {
unsigned int O,OP,S[150000],U[25000];
unsigned int H,K,N,D,NP,DP,NDEL,DDEL,HSAV,KSAV;
unsigned int L,I,J,f,del2,count,R,tind,maxzig,zig,ozig,tzig,oflag;
int del1;
double sum,sum2,temp,tmpsum,total,oldtot,T[25000],delta;
double pi;
FILE *Outfp;
Outfp = fopen("out1o.dat","w");
pi=3.141592654;
//
// ORDER OF FAREY SERIES
//
OP=57;			  // full order
O=OP/2; 		  // half-order
J=O;
oflag=0;
if ((O*2)!=OP) {
O=O+1;
oflag=1;
}
if (O>250) {
printf("order too big \n");
goto zskip;
}
tind=0; 		  // index of output array
maxzig=0;
ozig=0;
tzig=0;
L=haros1(O,0,S,0,1,1,O);  // generate half-order Farey series
S[L]=0;
sum=0.0;		  // half-order sum of cosines
sum2=0.0;		  // full order sum of cosines
for (I=0; I<=J; I++)
sum2=sum2+cos(2.0*pi/(double)(O+I));
oldtot=sum2;
for (I=1; I<(L-1); I++) {
H=S[I]>>16;
K=S[I]&0xffff;
HSAV=H;
KSAV=K;
tmpsum=0.0;
zig=0;
//
// begin loop
//
aloop:
f=lagrange1(H,K,OP);    // compute next fraction
N=f>>16;
D=f&0xffff;
if (f==S[I+1]) {
tmpsum=tmpsum+cos(2.0*pi*(double)N/(double)D);
goto bskip;
}
R=(OP+K)/D;	// compute next fraction
NP=R*N-H;
DP=R*D-K;
temp=cos(2.0*pi*(double)N/(double)D)+cos(2.0*pi*(double)NP/(double)DP);
if (D>DP) {
del1=(int)(2*DP)-(int)OP-1;
if (del1<=0)
del2=D-DP;
count=(unsigned int)del1/del2;
if ((unsigned int)del1!=((unsigned int)del1/del2)*del2)
count=count+1;
count=count/2;
NDEL=N-NP;
DDEL=D-DP;
for (J=0; J<count; J++) {
NP=NP-NDEL;
DP=DP-DDEL;
temp=temp+cos(2.0*pi*(double)NP/(double)DP);
}
}
else {
del1=(int)OP-(int)DP;
del2=DP-D;
count=(unsigned int)del1/del2;
NDEL=NP-N;
DDEL=DP-D;
for (J=0; J<count; J++) {
NP=NP+NDEL;
DP=DP+DDEL;
temp=temp+cos(2.0*pi*(double)NP/(double)DP);
}
}
tmpsum=tmpsum+temp;
if (((NP<<16)|DP)!=S[I+1]) {
H=NP;
K=DP;
zig=zig+1;
if (zig==2)
tzig=tzig+1;
else
ozig=ozig+1;
goto aloop;
}
//
// end loop
//
bskip:
if (zig>maxzig)
maxzig=zig;
sum2=sum2+tmpsum;
temp=cos(2.0*pi*(double)HSAV/(double)KSAV);
sum=sum+temp;
total=oldtot;
oldtot=sum2;
printf(" %d %d/%d %e %e \n",I,HSAV,KSAV,sum,total);
// fprintf(Outfp," %d %d/%d %e %e \n",I,HSAV,KSAV,sum,total);
T[tind]=total;
U[tind]=(HSAV<<16)|KSAV;
tind=tind+1;
}
if (oflag==0)
J=O;
else
J=O-1;
for (I=0; I<J; I++)
total=total+cos(2.0*pi*(double)(O+I)/(double)(O+I+1));
printf("order=%d, sums=%e, %e \n",O,sum,total);
//fprintf(Outfp,"order=%d, sums=%e, %e \n",O,sum,total);
//
// check results
//
printf("number of fractions (half order)=%d, maximum zigzags=%d \n",tind,maxzig);
printf("number of zigzag types=%d, %d \n",ozig,tzig);
J=0;
L=haros1(OP,0,S,0,1,1,OP);
sum=0.0;
for (I=0; I<L; I++) {
H=S[I]>>16;
K=S[I]&0xffff;
sum=sum+cos(2.0*pi*(double)H/(double)K);
if ((out==0)&&((double)H/(double)K<=0.5)&&(I<1000))
fprintf(Outfp," %e\n",sum);
if (S[I]!=U[J])
continue;
if ((out==1)&&((double)H/(double)K<=0.5)&&(I<1000))
fprintf(Outfp," %e\n",T[J]);
delta=sum-1.0-T[J];
if (delta<0)
delta=-delta;
if (delta>0.000001) {
printf("error: J=%d, fraction=%d/%d, sums=%e, %e \n",J,H,K,sum,T[J]);
break;
}
J=J+1;
}
printf("number of fractions (full order)=%d \n",L);
zskip:
return;
}
```