﻿ compute Mertens function
```/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE MERTENS FUNCTION (similar to E. Kuznetsov's algorithm)             C
C  01/27/16 (DKC)							      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
#include "table2.h"          // 17984 primes less than 200000
extern char *malloc();
// compute Mobius function
int newmob(unsigned int a, unsigned int b, int *out, unsigned int *table,
unsigned int tsize) {
unsigned int i,count,p,ps,beg,rem;
int t;
count=b-a;
for (i=0; i<count; i++)
out[i]=1;
for (i=0; i<tsize; i++) {
p=table[i];
ps=p*p;
if (ps>b)
rem=a-(a/ps)*ps;
if (rem!=0)
beg=ps-rem;
else
beg=0;
while (beg<count) {
out[beg]=0;
beg=beg+ps;
}
rem=a-(a/p)*p;
if (rem!=0)
beg=p-rem;
else
beg=0;
while (beg<count) {
out[beg]=-out[beg]*p;
beg=beg+p;
}
}
return(-1);
for (i=0; i<count; i++) {
if (out[i]==0)
continue;
t=out[i];
if (t<0)
t=-t;
if ((unsigned int)t<(i+a))
out[i]=-out[i];
}
for (i=0; i<count; i++) {
if (out[i]==0)
continue;
if (out[i]>0)
out[i]=1;
else
out[i]=-1;
}
return(1);
}
// compute partial sums of Mertens function
int newmert(unsigned int s, unsigned int x, int *M) {
unsigned int i,t,delta;
int sum;
t=(unsigned int)sqrt((double)x);
t=t+2;
if (s>(x/t))
return(0x7fffffff);
sum=0;
for (i=s; i<=(x/t); i++)
sum=sum+M[x/i-1];
for (i=1; i<t; i++) {		// suitable for implementation on a GPU
delta=x/i-x/(i+1);
sum=sum+delta*M[i-1];
}
return(sum);
}
//
unsigned int maxx=100000; // size of Mertens number look-up table (multiple of 1000)
unsigned int tmpmax=2000; // greater than or equal to 1000
unsigned int tsize=17984;  // size of prime look-up table

void main() {
unsigned int i,j,k,index,ta,tb,count,x,L,start;
int T[1000];
int t,sum;
char *m;
int *M,*temp;
FILE *Outfp;
Outfp = fopen("out19g.dat","w");
temp=(int*) malloc(4004);
if (temp==NULL) {
printf("not enough memory \n");
return;
}
m=(char*) malloc(maxx+1);
if (m==NULL) {
printf("not enough memory \n");
return;
}
M=(int*) malloc(4*(maxx+1));
if (M==NULL) {
printf("not enough memory \n");
return;
}
//
printf("computing Mobius function \n");
count=maxx;
count=count/1000;
index=0;
ta=1;
tb=1001;
for (i=0; i<count; i++) {
t=newmob(ta,tb,temp,table,tsize);
if (t!=1) {
printf("error \n");
goto zskip;
}
ta=tb;
tb=tb+1000;
for (j=0; j<1000; j++)
M[j+index]=temp[j];
index=index+1000;
}
//
// compute Mertens function
//
printf("computing Mertens function \n");
m[0]=(char)M[0];
for (i=1; i<=maxx; i++) {
m[i]=(char)M[i];
M[i]=M[i-1]+M[i];
}
//
// compute Mertens function
//
printf("comparing results \n");
for (i=(tmpmax-1000); i<=maxx; i++) {
x=i;
if (x<=tmpmax) {
t=newmert(2,x, M);
if (t==0x7fffffff) {
printf(" s>x/t \n");
goto zskip;
}
t=1-t;
}
else {
L=x/tmpmax;
if (L>1000) {
printf("not enough memory: 1 \n");
goto zskip;
}
for (j=1; j<=L; j++) {
start=(x/j)/tmpmax+1;
t=newmert(start,x/j,M);
if (t==0x7fffffff) {
printf(" s>x/t \n");
goto zskip;
}
T[j-1]=1-t;
}
for (j=(L/2); j>=1; j--) {
sum=0;
for (k=1; k<=(L/j-1); k++)
sum=sum+T[(k+1)*j-1];
T[j-1]=T[j-1]-sum;
}
t=T[0];
}
if (t!=M[x-1]) {
printf("error: x=%d t=%d, M[x-1]=%d \n",x,t,M[x-1]);
goto zskip;
}
if (i==(i/10000)*10000)
printf("x=%d, M(x)=%d \n",i,t);
}
zskip:
fclose(Outfp);
return;
}
```