/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C COMPUTE MERTENS FUNCTION (local maxima using Deleglise & Rivat algorithm) C
C 12/06/15 (DKC) C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
#include <omp.h>
#include "table.h"
extern char *malloc();
int newmobl(unsigned long long a, unsigned long long b, long long *out, unsigned int *table,
unsigned int tsize);
int newrivl(unsigned long long x, unsigned int u, char *mobb, int *M);
unsigned int primed(unsigned int *out, unsigned int tsize, unsigned int *table, unsigned int limit);
unsigned long long Msize=12000000000; // for 64 GB of RAM, 64-bit OS
unsigned int tsize=1230;
unsigned int tmpsiz=10000;
unsigned int t2size=200000;
//
unsigned int incnt=90;
unsigned long long in[119*2]={ // highly composite numbers and their number of divisors
2, 2,
4, 3,
6, 4,
12, 6,
24, 8,
36, 9,
48, 10,
60, 12,
120, 16,
180, 18,
240, 20,
360, 24,
720, 30,
840, 32,
1260, 36,
1680, 40,
2520, 48,
5040, 60,
7560, 64,
10080, 72,
15120, 80,
20160, 84,
25200, 90,
27720, 96,
45360, 100,
50400, 108,
55440, 120,
83160, 128,
110880, 144,
166320, 160,
221760, 168,
277200, 180,
332640, 192,
498960, 200,
554400, 216,
665280, 224,
720720, 240,
1081080, 256,
1441440, 288,
2162160, 320,
2882880, 336,
3603600, 360,
4324320, 384,
6486480, 400,
7207200, 432,
8648640, 448,
10810800, 480,
14414400, 504,
17297280, 512,
21621600, 576,
32432400, 600,
36756720, 640,
43243200, 672,
61261200, 720,
73513440, 768,
110270160, 800,
122522400, 864,
147026880, 896,
183783600, 960,
245044800, 1008,
294053760, 1024,
367567200, 1152,
551350800, 1200,
698377680, 1280,
735134400, 1344,
1102701600, 1440,
1396755360, 1536,
2095133040, 1600,
2205403200, 1680,
2327925600, 1728,
2793510720, 1792,
3491888400, 1920,
4655851200, 2016,
5587021440, 2048,
6983776800, 2304,
10475665200, 2400,
13967553600, 2688,
20951330400, 2880,
27935107200, 3072,
41902660800, 3360,
48886437600, 3456,
64250746560, 3584,
73329656400, 3600,
80313433200, 3840,
97772875200, 4032,
128501493120, 4096,
146659312800, 4320,
160626866400, 4608,
240940299600, 4800,
293318625600, 5040,
321253732800, 5376,
481880599200, 5760,
642507465600, 6144,
963761198400, 6720,
1124388064800, 6912,
1606268664000, 7168,
1686582097200, 7200,
1927522396800, 7680,
2248776129600, 8064,
3212537328000, 8192,
3373164194400, 8640,
4497552259200, 9216,
6746328388800, 10080,
8995104518400, 10368,
9316358251200, 10752,
13492656777600, 11520,
18632716502400, 12288,
26985313555200, 12960,
27949074753600, 13440,
32607253879200, 13824,
46581791256000, 14336,
48910880818800, 14400,
55898149507200, 15360,
65214507758400, 16128,
93163582512000, 16384,
97821761637600, 17280,
130429015516800, 18432,
195643523275200, 20160,
260858031033600, 20736};
//
void main() {
unsigned long long *oldtmp,*newtmp,MAXN;
long long *temp,sum;
int *M;
char *m;
unsigned int *ntable;
unsigned int *pritab;
unsigned int p,tindex,prind,delta,joff,newind,count,total,h,k,u,ntsize;
unsigned int g,f;
unsigned long long i,index,ta,tb,j,N,ut,pz,tz,mcount;
int savet,t,ID;
double f1,f2;
FILE *Outfp;
Outfp = fopen("out1ar.dat","w");
omp_set_dynamic(0); // CPU too hot for 8 threads
omp_set_num_threads(8);
#pragma omp parallel
{
ID=omp_get_thread_num();
printf(" ID=%d \n",ID);
}
MAXN=in[2*incnt-2];
f1=log(log((double)MAXN));
f1=f1*f1;
f1=exp(1.0/3.0*log(f1));
f2=exp(1.0/3.0*(log((double)MAXN)));
u=(unsigned int)(f1*f2);
printf("x/u=%I64, u=%d \n",MAXN/(unsigned long long)u,u);
if ((MAXN/(unsigned long long)u)>Msize) {
printf("not enough memory \n");
goto zskip;
}
ntable=(unsigned int*) malloc((t2size+1)*4);
if (ntable==NULL) {
printf("not enough memory \n");
goto zskip;
}
pritab=(unsigned int*) malloc(4000);
if (pritab==NULL) {
printf("not enough memory \n");
return;
}
temp=(long long*) malloc((tmpsiz+1)*8);
if (temp==NULL) {
printf("not enough memory \n");
return;
}
oldtmp=(long long*) malloc(160000);
if (oldtmp==NULL) {
printf("not enough memory \n");
return;
}
newtmp=(long long *) malloc(160000);
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;
}
ntsize=primed(ntable,tsize,table,2000000);
printf("prime look-up table size=%d, largest prime=%d \n",ntsize,ntable[ntsize-1]);
printf("computing Mobius function \n");
index=0;
ta=1;
tb=(unsigned long long)(tmpsiz+1);
mcount=Msize/(unsigned long long)tmpsiz;
for (i=0; i<mcount; i++) {
t=newmobl(ta,tb,temp,ntable,ntsize);
if (t!=1) {
printf("error \n");
goto zskip;
}
ta=tb;
tb=tb+(unsigned long long)tmpsiz;
for (j=0; j<tmpsiz; j++)
M[j+index]=(int)temp[j];
index=index+(unsigned long long)tmpsiz;
}
//
// 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=in[2*g];
ut=N;
prind=0;
tindex=0;
total=1;
h=(unsigned int)(sqrt((double)ut)+0.01); // for p>2, the difference between
for (f=0; f<ntsize; f++) { // successive primes is greater
p=table[f]; // 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>19999) {
printf("divisor table not big enough (1): 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>19999) {
printf("divisor table not big enough (2): 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!=(unsigned int)in[2*g+1]) {
printf("error: total=%d %d \n",total,(unsigned int)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 (f=0; f<(prind+1)/2; f++) {
count=pritab[2*f];
for (j=0; j<count; j++) {
newtmp[newind]=oldtmp[j+joff];
newind=newind+1;
}
for (j=0; j<pritab[2*f+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*f+1]; k++) {
pz=tz*oldtmp[k+count+joff];
newtmp[newind]=pz;
newind=newind+1;
if (newind>19999) {
printf("divisor table not big enough (3): N=%d \n",N);
goto zskip;
}
}
}
joff=joff+pritab[2*f]+pritab[2*f+1];
pritab[delta]=pritab[2*f]*pritab[2*f+1]+pritab[2*f]+pritab[2*f+1];
delta=delta+1;
}
for (f=0; f<newind; f++)
oldtmp[f]=newtmp[f];
pritab[delta]=0;
pritab[delta+1]=0;
prind=delta;
if (delta>1)
goto bskip;
//
// compute j(x)
//
cskip:
if ((newind+1)!=(unsigned int)in[2*g+1]) {
printf("error: newind=%d %d \n",newind,(unsigned int)in[2*g+1]);
goto zskip;
}
sum=1;
for (f=0; f<newind; f++) {
tz=oldtmp[f];
if (tz>Msize) {
f1=log(log((double)tz));
f1=f1*f1;
f1=exp(1.0/3.0*log(f1));
f2=exp(1.0/3.0*(log((double)tz)));
u=(unsigned int)(f1*f2);
t=newrivl(tz,u,m,M);
printf("x=%I64x, t=%d, u=%d x/u=%I64x \n",tz,t,u,N/(unsigned long long)u);
}
else {
t=M[tz-1];
}
if (tz==N)
savet=t;
sum=sum+(long long)t*(long long)t;
}
printf(" %I64x %I64x %d %d \n",N,sum,newind+1,savet);
fprintf(Outfp," %I64x, %I64x, %d, %d, \n",N,sum,newind+1,savet);
}
zskip:
fclose(Outfp);
return;
}