/*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)
	 goto askip;
      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);
	 }
      }
askip:
   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;
}