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