/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE MERTENS FUNCTION (sum of M(x/i)^2 where i|x, find local maxima)    C
C  09/29/15 (DKC) (maxima based on j(x))				      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
#include "table2.h"
extern char *malloc();
int newmob(unsigned int a, unsigned int b, int *out, unsigned int *table);
unsigned int MAXN=1000000000;  // maximum N
unsigned int BEGINN=2;	// beginning N
short k=18;
//
void main() {
int t,*temp;
short *M,t1;
unsigned int N,h,i,j,count,index,ta,tb;
double sum,maxt;
FILE *Outfp;
Outfp = fopen("out1bzs.dat","w");
temp=(int*) malloc(4004);
if (temp==NULL) {
   printf("not enough memory \n");
   return;
   }
M=(short*) malloc(2000000002);
if (M==NULL) {
   printf("not enough memory \n");
   return;
   }
printf("computing Mobius function \n");
index=0;
ta=1;
tb=1001;
for (i=0; i<(MAXN/1000); i++) {
   t=newmob(ta,tb,temp,table);
   if (t!=1) {
      printf("error \n");
      goto zskip;
      }
   ta=tb;
   tb=tb+1000;
   for (j=0; j<1000; j++)
      M[j+index]=(short)temp[j];
   index=index+1000;
   }
//
// compute Mertens function
//
printf("computing Mertens function \n");
for (i=1; i<MAXN; i++) {
   M[i]=M[i-1]+M[i];
   t1=(M[i]>>14)&3;
   if ((t1==1)||(t1==2)) {
      printf("possible overflow \n");
      printf(" %d %x \n",i,M[i]);
      goto zskip;
      }
   }
//
printf("finding maxima \n");
maxt=0.0;
for (N=BEGINN; N<=MAXN; N++) {
   t1=M[N-1];
   if (t1<0)
      t1=-t1;
   if (t1!=k)
      continue;
   sum=0.0;
   count=0;
   h=(unsigned int)sqrt((double)(N+1));
   for (i=1; i<=h; i++) {
      j=N/i;
      if (N==(j*i)) {
	 t=M[j-1];
	 sum=sum+(double)t*(double)t;
	 count=count+1;
	 if (j!=i) {
	    t=M[i-1];
	    sum=sum+(double)t*(double)t;
	    count=count+1;
	    }
	 }
      }
   if (sum>maxt) {
      maxt=sum;
      printf(" %d %e %d %d \n",N,maxt,count,M[N-1]);
      fprintf(Outfp," %d, %d, %d, %d, \n",N,(unsigned int)maxt,count,M[N-1]);
      }
   }
zskip:
fclose(Outfp);
return;
}