/******************************************************************************
* *
* COMPUTE MERTENS FUNCTION (Deleglise and Rivat's algorithm) *
* 10/16/15 (dkc) *
* *
******************************************************************************/
#include <math.h>
void div64_32(unsigned int *dividend, unsigned int *quotient,
unsigned int divisor);
double newloc(unsigned int x1, unsigned int x2, unsigned int u, char *mob, int *M) {
unsigned int i,j,k,m,X[2],Q[2],R[2];
double temp1,sumb;
int temp;
X[0]=x1;
X[1]=x2;
sumb=0.0;
for (m=1; m<=u; m++) {
temp=mob[m-1];
j=u/m+1;
div64_32(X,Q,m);
temp1=(double)Q[0]*65536.0*65536.0;
temp1=temp1+(double)Q[1];
k=(unsigned int)sqrt(temp1);
for (i=j; i<=k; i++) {
div64_32(Q,R,i);
if (temp<0)
sumb=sumb-(double)M[R[1]-1];
if (temp>0)
sumb=sumb+(double)M[R[1]-1];
}
}
return(sumb);
}