﻿ compute Mertens function
```/******************************************************************************
*									      *
*  COMPUTE MERTENS FUNCTION (Deleglise and Rivat's algorithm)                 *
*  10/18/15 (dkc)							      *
*									      *
******************************************************************************/
#include <math.h>
void div64_32(unsigned int *dividend, unsigned int *quotient,
unsigned int divisor);
void div64_64(unsigned int *dividend, unsigned int *quotient,
unsigned int divisor,unsigned int divlo);
void sub64(unsigned int *a, unsigned int *b);
double newhic(unsigned int x1, unsigned int x2, unsigned int u, char *mob, int *M) {
unsigned int j,m,X[2],Q[2],S[2],I[2],U[2],index;
double temp1,sumb;
int temp;
X[0]=x1;
X[1]=x2;
sumb=0.0;
for (m=1; m<=u; m++) {
temp=mob[m-1];
div64_32(X,Q,m);
temp1=(double)Q[0]*65536.0*65536.0;
temp1=temp1+(double)Q[1];
j=(unsigned int)sqrt(temp1);
j=j+1;
I[0]=0;
I[1]=j;
S[0]=I[0];
S[1]=I[1];
sub64(Q,S);
while ((S[0]>>31)==0) {
if (I[0]!=0) {
div64_64(Q,U,I[0],I[1]);
index=U[1];
}
else {
if (Q[0]!=0) {
div64_32(Q,U,I[1]);
index=U[1];
}
else
index=Q[1]/I[1];
}
if (temp<0)
sumb=sumb-(double)M[index-1];
if (temp>0)
sumb=sumb+(double)M[index-1];
I[1]=I[1]+1;
if (I[1]==0)
I[0]=I[0]+1;
S[0]=I[0];
S[1]=I[1];
sub64(Q,S);
}
}
return(sumb);
}
```