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