/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C COMPUTE MERTENS FUNCTION (local maxima using Deleglise & Rivat algorithm) C
C 10/18/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);
int fastriv(unsigned int x1, unsigned int u, char *mobb, int *M);
void mul64_32(unsigned int a0, unsigned int a2, unsigned int m0,
unsigned int *product);
//
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,*newtmp;
int *M;
char *m;
unsigned int MAXN,P[3],T[2],U[2],g,i,j,k,count,u,index,ta,tb;
unsigned int flag2,flag3,flag5,ut,total,tindex,p;
unsigned int pritab[100],prind,delta,joff,newind;
double sum,f1,f2;
FILE *Outfp;
Outfp = fopen("out1ah.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)>400000000) {
printf("not enough memory \n");
goto zskip;
}
//
temp=(int*) malloc(32000);
if (temp==NULL) {
printf("not enough memory \n");
return;
}
newtmp=(int*) malloc(32000);
if (newtmp==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]);
//
// factor N
//
for (g=0; g<incnt; g++) {
ut=in[2*g];
prind=0;
total=0;
flag2=3;
temp[0]=2;
temp[1]=4;
temp[2]=8;
tindex=3;
count=0;
while ((ut&1)==0) {
temp[tindex]=temp[tindex-1]*2;
tindex=tindex+1;
count=count+1;
ut=ut>>1;
}
flag2=flag2+count;
pritab[prind]=flag2;
prind=prind+1;
flag3=2;
temp[tindex]=3;
temp[tindex+1]=9;
tindex=tindex+2;
count=0;
while (ut==(ut/3)*3) {
temp[tindex]=temp[tindex-1]*3;
tindex=tindex+1;
count=count+1;
ut=ut/3;
}
flag3=flag3+count;
pritab[prind]=flag3;
prind=prind+1;
flag5=1;
temp[tindex]=5;
tindex=tindex+1;
count=0;
while (ut==(ut/5)*5) {
temp[tindex]=temp[tindex-1]*5;
tindex=tindex+1;
if (tindex>3999) {
printf("divisor table not big enough \n");
goto zskip;
}
count=count+1;
ut=ut/5;
}
flag5=flag5+count;
pritab[prind]=flag5;
prind=prind+1;
total=(flag2+1)*(flag3+1)*(flag5+1);
if (ut==1)
goto askip;
for (i=3; i<17983; i++) {
p=table[i];
count=0;
while (ut==(ut/p)*p) {
if (count==0)
temp[tindex]=p;
else
temp[tindex]=temp[tindex-1]*p;
tindex=tindex+1;
if (tindex>3999) {
printf("divisor table not big enough \n");
goto zskip;
}
ut=ut/p;
count=count+1;
}
if (count!=0) {
total=total*(count+1);
pritab[prind]=count;
prind=prind+1;
if (prind>49) {
printf("prime table not big enough \n");
goto zskip;
}
}
if (ut==1)
goto askip;
}
printf("error: prime look-up table not big enough \n");
goto zskip;
//
// compute combinations of factors
//
askip:
if (total!=in[2*g+1]) {
printf("error: total=%d %d \n",total,in[2*g+1]);
goto zskip;
}
for (i=(tindex-1); i>0; i--) {
temp[2*i+1]=temp[i];
temp[2*i]=0;
}
temp[1]=temp[0];
temp[0]=0;
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*2]=temp[2*(j+joff)];
newtmp[newind*2+1]=temp[2*(j+joff)+1];
newind=newind+1;
}
for (j=0; j<pritab[2*i+1]; j++) {
newtmp[newind*2]=temp[2*(j+joff+count)];
newtmp[newind*2+1]=temp[2*(j+joff+count)+1];
newind=newind+1;
}
for (j=0; j<count; j++) {
T[0]=temp[2*(j+joff)];
T[1]=temp[2*(j+joff)+1];
for (k=0; k<pritab[2*i+1]; k++) {
U[0]=temp[2*(k+count+joff)];
U[1]=temp[2*(k+count+joff)+1];
if (T[0]==0)
mul64_32(U[0],U[1],T[1],P);
else
mul64_32(T[0],T[1],U[1],P);
newtmp[newind*2]=P[1];
newtmp[newind*2+1]=P[2];
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++) {
temp[i*2]=newtmp[i*2];
temp[i*2+1]=newtmp[i*2+1];
}
pritab[delta]=0;
pritab[delta+1]=0;
prind=delta;
if (delta>1)
goto bskip;
if ((newind+1)!=in[2*g+1]) {
printf("error: total=%d %d \n",newind,in[2*g+1]);
goto zskip;
}
//
// compute j(x)
//
sum=1.0;
for (i=0; i<newind; i++) {
T[0]=temp[i*2];
T[1]=temp[i*2+1];
if (T[0]!=0) {
t=newriv(T[0],T[1],u,m,M);
printf("x=%#010x %#010x, t=%d \n",T[0],T[1],t);
}
else {
if (T[1]>400000000) {
t=fastriv(T[1],u,m,M);
printf("x=%#010x, t=%d \n",T[1],t);
}
else
t=M[T[1]-1];
}
sum=sum+(double)t*(double)t;
}
mul64_32(0,in[2*g],360,P);
printf(" %#010x %#010x %d %d \n",P[1],P[2],(unsigned int)sum,newind+1);
fprintf(Outfp," %#010x, %#010x, %d, %d, \n",P[1],P[2],(unsigned int)sum,newind+1);
}
zskip:
fclose(Outfp);
return;
}