/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C COMPUTE MERTENS FUNCTION (j(x) where x is a highly composite number) C
C 01/05/16 (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 newmert(unsigned int s, unsigned long long x, int *M);
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=15000000000; // for 64 GB of RAM, 64-bit OS
unsigned int tsize=1230;
unsigned int tmpsiz=10000;
unsigned int t2size=200000;
unsigned int Tsize=10000000;
unsigned int oldnews=100000;
unsigned int prsize=1000;
//
unsigned int incnt=167;
unsigned long long in[167*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,
288807105787200, 21504,
391287046550400, 23040,
577614211574400, 24576,
782574093100800, 25920,
866421317361600, 26880,
1010824870255200, 27648,
1444035528936000, 28672,
1516237305382800, 28800,
1732842634723200, 30720,
2021649740510400, 32256,
2888071057872000, 32768,
3032474610765600, 34560,
4043299481020800, 36864,
6064949221531200, 40320,
8086598962041600, 41472,
10108248702552000, 43008,
12129898443062400, 46080,
18194847664593600, 48384,
20216497405104000, 49152,
24259796886124800, 51840,
30324746107656000, 53760,
36389695329187200, 55296,
48519593772249600, 57600,
60649492215312000, 61440,
72779390658374400, 62208,
74801040398884800, 64512,
106858629141264000, 65536,
112201560598327200, 69120,
149602080797769600, 73728,
224403121196654400, 80640,
299204161595539200, 82944,
374005201994424000, 86016,
448806242393308800, 92160,
673209363589963200, 96768,
748010403988848000, 98304,
897612484786617600, 103680,
1122015605983272000, 107520,
1346418727179926400, 110592,
1795224969573235200, 115200,
2244031211966544000, 122880,
2692837454359852800, 124416,
3066842656354276800, 129024,
4381203794791824000, 131072,
4488062423933088000, 138240,
6133685312708553600, 147456,
8976124847866176000, 153600,
9200527969062830400, 161280,
12267370625417107200, 165888};
//
void main() {
unsigned long long *oldtmp,*newtmp;
long long *temp,sum,sump;
int *M,*T;
unsigned int *ntable;
unsigned int *pritab;
unsigned int p,tindex,prind,delta,joff,newind,count,total,h,k,ntsize;
unsigned int g,f,L,start;
int d,e;
unsigned long long index,ta,tb,j,N,ut,pz,tz,mcount;
long long i;
int savet,t,ID;
FILE *Outfp;
Outfp = fopen("out1arp.dat","w");
omp_set_dynamic(0); // CPU too hot for 8 threads
omp_set_num_threads(4);
#pragma omp parallel
{
ID=omp_get_thread_num();
printf(" ID=%d \n",ID);
}
ntable=(unsigned int*) malloc((t2size+1)*4);
if (ntable==NULL) {
printf("not enough memory \n");
goto zskip;
}
pritab=(unsigned int*) malloc((prsize+1)*4);
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((oldnews+1)*8);
if (oldtmp==NULL) {
printf("not enough memory \n");
return;
}
newtmp=(long long *) malloc((oldnews+1)*8);
if (newtmp==NULL) {
printf("not enough memory \n");
return;
}
M=(int*) malloc((Msize+1)*4);
if (M==NULL) {
printf("not enough memory \n");
return;
}
T = (int*)malloc((Tsize+1)*4);
if (T == NULL) {
printf("not enough memory: 5 \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<(long long)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");
for (i=1; i<=(long long)Msize; 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>oldnews) {
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>prsize) {
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>oldnews) {
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>prsize) {
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>oldnews) {
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=0;
if (N>Msize) {
L=(unsigned int)(N/(unsigned long long)Msize);
if (L>(unsigned long long)Tsize) {
printf("error: not enough memory \n");
goto zskip;
}
for (e=1; e<=(int)L; e++) {
start=(unsigned int)((N/(unsigned long long)e)/(unsigned long long)Msize)+1;
T[e-1]=1-newmert(start, N/e, M);
}
for (e=L; e>1; e--) {
#pragma omp parallel for default(none) shared(e,T) private(d)
for (d=(e-1); d>=1; d--) {
if (e==(e/d)*d)
T[d-1]=T[d-1]-T[e-1];
}
}
savet = T[0];
#pragma omp parallel for reduction (+:sum)
for (e=1; e<=(int)L; e++)
if (N==(N/(unsigned long long)e)*(unsigned long long)e)
sum=sum+(long long)T[e-1]*(long long)T[e-1];
}
sump=1;
for (f=0; f<newind; f++) {
tz=oldtmp[f];
if (tz<=Msize) {
t=M[tz-1];
sump=sump+(long long)t*(long long)t;
if (tz==N)
savet=t;
}
}
sum=sum+sump;
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;
}