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