/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE MERTENS FUNCTION (cosines only, first two quadrants) 	      C
C  05/12/14 (DKC)							      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);
//
void main() {
unsigned int N,S[250000];
unsigned int H,K,HP,KP,MAXN;
unsigned int L,I,index,count1,cnt1sav,count2,cnt2sav;
double temp1,temp2,sum,sum1,sum2,sum1sav,sum2sav;
double pi,ratio1,ratio2;
FILE *Outfp;
Outfp = fopen("out1p.dat","w");
pi=3.141592654;
//
// ORDER OF FAREY SERIES
//
MAXN=900;
N=MAXN;
while (N>4) {
   L=haros1(N,0,S,0,1,1,N);
   index=0;
   for (I=0; I<L; I++) {	 // find index of 1/4
      H=S[I]>>16;
      K=S[I]&0xffff;
      if ((double)H/(double)K>=0.25)
	 break;
      index=index+1;
      }
   count1=index;
   count2=0;
   sum1=0.0;
   sum2=0.0;
   for (I=index; I>0; I--) {
      H=S[I-1]>>16;		      // towards 0/1
      K=S[I-1]&0xffff;
      HP=S[index+1+count2]>>16; 	 // towards 1/2
      KP=S[index+1+count2]&0xffff;
      temp1=cos(2.0*pi*(double)H/(double)K);
      sum1=sum1+temp1;
      if (((double)HP/(double)KP)<0.5) {
	 temp2=cos(2.0*pi*(double)HP/(double)KP);
	 sum2=sum2+temp2;
	 count2=count2+1;
	 }
//    printf(" %d/%d %e %d/%d %e \n",H,K,sum1,HP,KP,sum2);
      }
   count1=count1-1;   // exclude 0/1
   sum1=sum1-1.0;
   I=2*index+1;
   HP=S[I]>>16;
   KP=S[I]&0xffff;
// printf(" %d %d \n",HP,KP);
   while (((double)HP/(double)KP)<0.5) {
      temp2=cos(2.0*pi*(double)HP/(double)KP);
      sum2=sum2+temp2;
      I=I+1;
      HP=S[I]>>16;
      KP=S[I]&0xffff;
//    printf(" %d %d \n",HP,KP);
      count2=count2+1;
      }
   sum=sum1+sum2;
   ratio1=(sum1/sum1sav)*(double)cnt1sav/(double)count1;
   ratio2=(sum2/sum2sav)*(double)cnt2sav/(double)count2;
   if (N!=MAXN) {
      printf("N=%d, s1=%e %e, c1=%d %d \n",N,sum1sav,sum1,cnt1sav,count1);
      printf("      s2=%e %e, c2=%d %d \n",sum2sav,sum2,cnt2sav,count2);
      printf("       r=%e %e, s=%e \n",ratio1,ratio2,2.0*sum);
      printf("\n");
      }
   if (N!=(N/2)*2)
      N=(N+1)/2;
   else
      N=N/2;
   sum1sav=sum1;
   sum2sav=sum2;
   cnt1sav=count1;
   cnt2sav=count2;
   }
return;
}