/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C COMPUTE MERTENS FUNCTION (local maxima using Deleglise & Rivat algorithm) C
C 10/31/15 (DKC) 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);
int newrivl(unsigned long long x, unsigned int u, char *mobb, int *M);
unsigned long long Msize=3000000000; // for 16 GB of RAM, 64-bit OS
//
unsigned int incnt=73;
unsigned int in[73*2]={
0x00000007, 48,
0x0000000e, 60,
0x00000015, 64,
0x0000001c, 72,
0x0000002a, 80,
0x00000038, 84,
0x00000046, 90,
0x0000004d, 96,
0x0000007e, 100,
0x0000008c, 108,
0x0000009a, 120,
0x000000e7, 128,
0x00000134, 144,
0x000001ce, 160,
0x00000268, 168,
0x00000302, 180,
0x0000039c, 192,
0x0000056a, 200,
0x00000604, 216,
0x00000738, 224,
0x000007d2, 240,
0x00000bbb, 256,
0x00000fa4, 288,
0x00001776, 320,
0x00001f48, 336,
0x0000271a, 360,
0x00002eec, 384,
0x00004662, 400,
0x00004e34, 432,
0x00005dd8, 448,
0x0000754e, 480,
0x00009c68, 504,
0x0000bbb0, 512,
0x0000ea9c, 576,
0x00015fea, 600,
0x00018ed6, 640,
0x0001d538, 672,
0x000298ba, 720,
0x00031dac, 768,
0x0004ac82, 800,
0x00053174, 864,
0x00063b58, 896,
0x0007ca2e, 960,
0x000a62e8, 1008,
0x000c76b0, 1024,
0x000f945c, 1152,
0x00175e8a, 1200,
0x001d99e2, 1280,
0x001f28b8, 1344,
0x002ebd14, 1440,
0x003b33c4, 1536,
0x0058cda6, 1600,
0x005d7a28, 1680,
0x0062ab9c, 1728,
0x00766788, 1792,
0x0094016a, 1920,
0x00c55738, 2016,
0x00eccf10, 2048,
0x012802d4, 2304,
0x01bc043e, 2400,
0x025005a8, 2688,
0x0378087c, 2880,
0x04a00b50, 3072,
0x06f010f8, 3360,
0x081813cc, 3456,
0x0aa34d38, 3584,
0x0c241db2, 3600,
0x0d4c2086, 3840,
0x10302798, 4032,
0x15469a70, 4096,
0x18483b64, 4320,
0x1a98410c, 4608,
0x27e46192, 4800};
//
void main() {
unsigned long long *oldtmp,*newtmp;
int *temp;
int *M;
char *m;
unsigned int p,tindex,prind,delta,joff,newind,count,total,h,k,u,pritab[1000];
unsigned int MAXN,g;
unsigned long long index,ta,tb,i,j,N,ut,pz,tz;
int savet,t;
double sum,f1,f2;
FILE *Outfp;
Outfp = fopen("out1ar.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, u=%d \n",(MAXN/u)*360,u);
if (((MAXN/u)*360)>Msize) {
printf("not enough memory \n");
goto zskip;
}
temp=(int*) malloc(4004);
if (temp==NULL) {
printf("not enough memory \n");
return;
}
oldtmp=(long long*) malloc(64000);
if (oldtmp==NULL) {
printf("not enough memory \n");
return;
}
newtmp=(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;
}
m=(char*) malloc(Msize+1);
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");
m[0]=(char)M[0];
for (i=1; i<=Msize; i++) {
m[i]=(char)M[i];
M[i]=M[i-1]+M[i];
}
//
printf("finding maxima \n");
for (g=0; g<incnt; g++) {
//
// factor N
//
N=(unsigned long long)in[2*g]*360;
ut=N;
prind=0;
tindex=0;
total=1;
h=(unsigned int)(sqrt((double)ut)+0.01); // for p>2, the difference between
for (i=0; i<41538; i++) { // successive primes is greater
p=table[i]; // than 1, so adding 0.01 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 (total!=in[2*g+1]) {
printf("error: total=%d %d \n",total,in[2*g+1]);
goto zskip;
}
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:
if ((newind+1)!=in[2*g+1]) {
printf("error: newind=%d %d \n",newind,in[2*g+1]);
goto zskip;
}
sum=1.0;
for (i=0; i<newind; i++) {
tz=oldtmp[i];
if (tz>Msize) {
t=newrivl(tz,u,m,M);
if (tz==N)
savet=t;
printf("x=%I64x, t=%d \n",tz,t);
}
else {
t=M[tz-1];
if (tz==N)
savet=t;
}
sum=sum+(double)t*(double)t;
}
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;
}