/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE MERTENS FUNCTION (local maxima using Deleglise and Rivat)	      C
C  10/16/15 (DKC)							      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);
int newriv(unsigned int x1, unsigned int x2, unsigned int u, char *mobb, int *M);
void div64_32(unsigned int *dividend, unsigned int *quotient,
	      unsigned int divisor);
void mul64_32(unsigned int a0, unsigned int a2, unsigned int m0,
	      unsigned int *product);
void sub64(unsigned int *a, unsigned int *b);
//
unsigned int incnt=62;
unsigned int in[62*2]={
 14, 60, 
 21, 64, 
 28, 72, 
 42, 80, 
 56, 84, 
 70, 90, 
 77, 96, 
 126, 100, 
 140, 108, 
 154, 120, 
 231, 128, 
 308, 144, 
 462, 160, 
 616, 168, 
 770, 180, 
 924, 192, 
 1386, 200, 
 1540, 216, 
 1848, 224, 
 2002, 240, 
 3003, 256, 
 4004, 288, 
 6006, 320, 
 8008, 336, 
 10010, 360, 
 12012, 384, 
 18018, 400, 
 20020, 432, 
 24024, 448, 
 30030, 480, 
 40040, 504, 
 48048, 512, 
 60060, 576, 
 90090, 600, 
 102102, 640, 
 120120, 672, 
 170170, 720, 
 204204, 768, 
 306306, 800, 
 340340, 864, 
 408408, 896, 
 510510, 960, 
 680680, 1008, 
 816816, 1024, 
 1021020, 1152, 
 1531530, 1200, 
 1939938, 1280, 
 2042040, 1344, 
 3063060, 1440, 
 3879876, 1536, 
 5819814, 1600, 
 6126120, 1680, 
 6466460, 1728, 
 7759752, 1792, 
 9699690, 1920, 
 12932920, 2016, 
 15519504, 2048, 
 19399380, 2304, 
 29099070, 2400, 
 38798760, 2688, 
 58198140, 2880, 
 77597520, 3072};
//
void main() {
int t,*temp,savet;
int *M;
char *m;
unsigned int MAXN,N[2],P[3],T[2],J[2],g,h,i,j,count,u,index,ta,tb;
double sum,f1,f2;
FILE *Outfp;
Outfp = fopen("out1af.dat","w");
MAXN=in[2*incnt-2];
f1=log(log((double)MAXN)+log(360.0));
f1=f1*f1;
f1=exp(1.0/3.0*log(f1));
f2=exp(1.0/3.0*(log((double)MAXN)+log(360.0)));
u=(unsigned int)(f1*f2);
printf("x/u=%d \n",(MAXN/u)*360);
if (((MAXN/u)*360)>400000000) {
   printf("not enough memory \n");
   goto zskip;
   }
//
temp=(int*) malloc(4004);
if (temp==NULL) {
   printf("not enough memory \n");
   return;
   }
M=(int*) malloc(1600000004);
if (M==NULL) {
   printf("not enough memory \n");
   return;
   }
m=(char*) malloc(400000001);
if (m==NULL) {
   printf("not enough memory \n");
   return;
   }
printf("computing Mobius function \n");
index=0;
ta=1;
tb=1001;
for (i=0; i<400000; 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]=temp[j];
   index=index+1000;
   }
//
// compute Mertens function
//
printf("computing Mertens function \n");
m[0]=(char)M[0];
for (i=1; i<=400000000; i++) {
   m[i]=(char)M[i];
   M[i]=M[i-1]+M[i];
   }
mul64_32(0,MAXN,360,P);
printf("computing sum: MAXN=%#010x %#010x \n",P[1],P[2]);
//
for (g=0; g<incnt; g++) {
   mul64_32(0,in[2*g],360,P);
   N[0]=P[1];
   N[1]=P[2];
   sum=0.0;
   count=0;
   h=(unsigned int)(sqrt((double)in[2*g])*sqrt(360.0));
   for (i=1; i<=h; i++) {
      div64_32(N,J,i);
      mul64_32(J[0],J[1],i,P);
      T[0]=P[1];
      T[1]=P[2];
      sub64(N,T);
      if ((T[0]==0)&&(T[1]==0)) {
	 j=J[1];
	 if ((j>400000000)||(J[0]!=0)) {
	    t=newriv(J[0],J[1],u,m,M);
	    printf("x=%#010x %#010x, t=%d \n",J[0],J[1],t);
	    }
	 else
	    t=M[j-1];
	 if (i==1)
	    savet=t;
	 sum=sum+(double)t*(double)t;
	 count=count+1;
	 if ((j!=i)||(J[0]!=0)) {
	    t=M[i-1];
	    sum=sum+(double)t*(double)t;
	    count=count+1;
	    }
	 }
      }
   if (count!=in[2*g+1]) {
      printf("error: count=%d %d \n",count,in[2*g+1]);
      goto zskip;
      }
   if (savet<0)
      savet=-savet;
   printf(" %#010x %#010x %d %d %d \n",N[0],N[1],(unsigned int)sum,count,savet);
   fprintf(Outfp," %#010x, %#010x, %d, %d, %d, \n",N[0],N[1],(unsigned int)sum,count,savet);
   }
zskip:
fclose(Outfp);
return;
}