/*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);
unsigned int haros3(unsigned int N, unsigned int M, unsigned int *R,
		    unsigned int H, unsigned int K, unsigned int HP,
		    unsigned int KP);
unsigned int haros4(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);
double mertens4(unsigned int N, unsigned int *count, double *sum) {
unsigned int S[250000];
unsigned int H,K;
unsigned int L,I,count1,count2,ctemp1,ctemp2;
double temp1,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==5) {
   count[0]=1;
   count[1]=2;
   sum[0]=0.309017;
   sum[1]=-1.309017;
   return(-2.0);
   }
if (N==6) {
   count[0]=2;
   count[1]=2;
   sum[0]=0.809017;
   sum[1]=-1.309017;
   return(-1.0);
   }
if (N==7) {
   count[0]=3;
   count[1]=4;
   sum[0]=1.432507;
   sum[1]=-2.432507;
   return(-2.0);
   }
if (N==8) {
   count[0]=4;
   count[1]=5;
   sum[0]=2.139614;
   sum[1]=-3.139614;
   return(-2.0);
   }
if (N>2560) {
   count[0]=0;
   count[1]=0;
   sum[0]=0.0;
   sum[1]=0.0;
   return(0.0);
   }
pi=3.141592654;
ctemp1=haros4(N,0,S,0,1,1,N);
ctemp1=ctemp1-1;		   // exclude 0/1
sum1=0.0;
for (I=1; I<ctemp1; I++) {
   H=S[I]>>16;			 // towards 1/8
   K=S[I]&0xffff;
   temp1=cos(2.0*pi*(double)H/(double)K);
   sum1=sum1+temp1;
   }
L=lagrange1(1,8,N);
count1=haros3(N,0,S,1,8,L>>16,L&0xffff);
for (I=0; I<count1; I++) {
   H=S[I]>>16;			 // towards 1/4
   K=S[I]&0xffff;
   temp1=cos(2.0*pi*(double)H/(double)K);
   sum1=sum1+temp1;
   }
count1=count1+ctemp1-2;
count[0]=count1;
sum[0]=sum1;
//
L=lagrange1(1,4,N);
ctemp2=haros4(N,0,S,1,4,L>>16,L&0xffff);
sum2=0.0;
for (I=1; I<ctemp2; I++) {
   H=S[I]>>16;		       // towards 3/8
   K=S[I]&0xffff;
   temp1=cos(2.0*pi*(double)H/(double)K);
   sum2=sum2+temp1;
   }
L=lagrange1(3,8,N);
count2=haros2(N,0,S,3,8,L>>16,L&0xffff);
count2=count2-1;
for (I=1; I<count2; I++) {
   H=S[I]>>16;		       // towards 1/2
   K=S[I]&0xffff;
   temp1=cos(2.0*pi*(double)H/(double)K);
   sum2=sum2+temp1;
   }
count2=count2+ctemp2-2;
count[1]=count2;
sum[1]=sum2;
return(2*(sum1+sum2));
}