﻿ compute Mertens function
```/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE MERTENS FUNCTION						      C
C  05/30/14 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <math.h>
unsigned int haros2(unsigned int N, unsigned int M, unsigned int *R,
unsigned int H, unsigned int K, unsigned int HP,
unsigned int KP);
double mertens2(unsigned int N, unsigned int *count, double *sum) {
unsigned int S[250000];
unsigned int H,K,HP,KP;
unsigned int L,I,index,count1,count2;
double temp1,temp2,sum1,sum2;
double pi;
if (N==1) {
count[0]=0;
count[1]=0;
sum[0]=0.0;
sum[1]=0.0;
return(1.0);
}
if (N==2) {
count[0]=0;
count[1]=0;
sum[0]=0.0;
sum[1]=0.0;
return(0.0);
}
if (N==3) {
count[0]=0;
count[1]=0;
sum[0]=0.0;
sum[1]=0.0;
return(-1.0);
}
if (N==4) {
count[0]=0;
count[1]=1;
sum[0]=0.0;
sum[1]=-0.5;
return(-1.0);
}
if (N>1280) {
count[0]=0;
count[1]=0;
sum[0]=0.0;
sum[1]=0.0;
return(0.0);
}
pi=3.141592654;
L=haros2(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;
}
}
count1=count1-1;   // exclude 0/1
sum1=sum1-1.0;
I=2*index+1;
HP=S[I]>>16;
KP=S[I]&0xffff;
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;
count2=count2+1;
}
count[0]=count1;
count[1]=count2;
sum[0]=sum1;
sum[1]=sum2;
return(2*(sum1+sum2));
}
```