﻿ approximate Mertens function
```/*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;
}
```