/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE MERTENS FUNCTION (sum of M(x/i)^2 where i|x, find local maxima)    C
C  10/30/15 (DKC) (maxima based on j(x))				      C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
#include "table3.h"
extern char *malloc();
int newmobl(unsigned long long a, unsigned long long b, int *out, unsigned int *table);
unsigned long long MAXN=3500000000; // (2^31-1=4294967295)
unsigned long long BEGINN=2;  // beginning N
unsigned long long Msize=3500000000;   // for 16 GB of RAM, 64-bit OS
//
void main() {
int t,*temp;
int *M;
unsigned int total,count,tindex,p,h;
unsigned int pritab[1000],prind,delta,joff,newind;
unsigned int j,k;
int savet;
unsigned long long i,tz,pz;
unsigned long long *oldtmp,*newtmp;
unsigned long long N,ta,tb,index;
unsigned long long ut;
double sum,maxt;
FILE *Outfp;
Outfp = fopen("out1ap.dat","w");
temp=(int *) malloc(4004);
if (temp==NULL) {
   printf("not enough memory \n");
   return;
   }
oldtmp=(unsigned long long*) malloc(64000);
if (oldtmp==NULL) {
   printf("not enough memory \n");
   return;
   }
newtmp=(unsigned long long*) malloc(64000);
if (newtmp==NULL) {
   printf("not enough memory \n");
   return;
   }
M=(int*) malloc((Msize+1)*4);
if (M==NULL) {
   printf("not enough memory \n");
   return;
   }
printf("computing Mobius function \n");
index=0;
ta=1;
tb=1001;
for (i=0; i<(Msize/1000); i++) {
   t=newmobl(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]=temp[j];
   index=index+1000;
   }
//
// compute Mertens function
//
printf("computing Mertens function \n");
for (i=1; i<=Msize; i++)
   M[i]=M[i-1]+M[i];
//
printf("finding maxima \n");
maxt=0.0;
for (N=BEGINN; N<=MAXN; N++) {
   if ((N>13)&&((N&1)==1))     // Note!  Maxima appear to be even for N>13.
      continue;
//
// factor N
//
   ut=N;
   prind=0;
   tindex=0;
   total=0;
   h=(unsigned int)(sqrt((double)ut)+0.001); // For p>2, the difference between
   for (i=0; i<41538; i++) {		     // successive primes is greater
      p=table[i];			     // than 1, so adding 0.001 is okay.
      if (p>h)
	 goto fskip;
      count=0;
      while (ut==(ut/p)*p) {
	 if (count==0)
	    oldtmp[tindex]=p;
	 else
	    oldtmp[tindex]=oldtmp[tindex-1]*p;
	 tindex=tindex+1;
	 if (tindex>7999) {
	    printf("divisor table not big enough: N=%d \n",N);
	    goto zskip;
	    }
	 ut=ut/p;
	 count=count+1;
	 }
      if (count!=0) {
	 total=total*(count+1);
	 pritab[prind]=count;
	 prind=prind+1;
	 if (prind>1000) {
	    printf("prime table not big enough: N=%d \n",N);
	    goto zskip;
	    }
	 }
      if (ut==1)
	 goto askip;
      }
   printf("error: prime look-up table not big enough \n");
   goto zskip;
//
//  compute combinations of factors
//
fskip:
   oldtmp[tindex]=ut;
   tindex=tindex+1;
   if (tindex>7999) {
      printf("divisor table not big enough: N=%d \n",N);
      goto zskip;
      }
   count=1;
   total=total*(count+1);
   pritab[prind]=count;
   prind=prind+1;
   if (prind>1000) {
       printf("prime table not big enough: N=%d \n",N);
       goto zskip;
       }
askip:
   if (prind==1) {
      newind=tindex;
      goto cskip;
      }
   delta=0;
   pritab[prind]=0;
   pritab[prind+1]=0;
bskip:
   joff=0;
   delta=0;
   newind=0;
   for (i=0; i<(prind+1)/2; i++) {
      count=pritab[2*i];
      for (j=0; j<count; j++) {
	 newtmp[newind]=oldtmp[j+joff];
	 newind=newind+1;
	 }
      for (j=0; j<pritab[2*i+1]; j++) {
	 newtmp[newind]=oldtmp[j+joff+count];
	 newind=newind+1;
	 }
      for (j=0; j<count; j++) {
	 tz=oldtmp[j+joff];
	 for (k=0; k<pritab[2*i+1]; k++) {
	    pz=tz*oldtmp[k+count+joff];
	    newtmp[newind]=pz;
	    newind=newind+1;
	    if (newind>3999) {
	       printf("divisor table not big enough \n");
	       goto zskip;
	       }
	    }
	 }
      joff=joff+pritab[2*i]+pritab[2*i+1];
      pritab[delta]=pritab[2*i]*pritab[2*i+1]+pritab[2*i]+pritab[2*i+1];
      delta=delta+1;
      }
   for (i=0; i<newind; i++)
      oldtmp[i]=newtmp[i];
   pritab[delta]=0;
   pritab[delta+1]=0;
   prind=delta;
   if (delta>1)
      goto bskip;
//
// compute j(x)
//
cskip:
   sum=1.0;
   for (i=0; i<newind; i++) {
      tz=oldtmp[i];
      t=M[tz-1];
      if (tz==N)
	 savet=t;
      sum=sum+(double)t*(double)t;
      }
   if (sum>maxt) {
      maxt=sum;
      printf(" %d %I64x %d %I64x %d %d \n",N,N,(unsigned int)sum,(unsigned long long)sum,newind+1,savet);
      fprintf(Outfp," %I64x, %I64x, %d, %d, \n",N,(unsigned long long)sum,newind+1,savet);
      }
   }
zskip:
fclose(Outfp);
return;
}