/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C COMPUTE MERTENS FUNCTION (j(x) where x is a highly composite number) C
C 01/15/16 (DKC) C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
#include <omp.h>
#include "table.h"
extern char *malloc();
// compute partial sums of Mertens function
long long newmert(unsigned long long s, unsigned long long x, int *M) {
unsigned long long t;
long long i;
long long sum;
t=(unsigned long long)sqrt((double)x);
t=t+2;
if (s>(x/t))
return(0x7fffffffffffffff);
sum=0;
#pragma omp parallel for reduction (+:sum)
for (i=(long long)s; i<=(long long)(x/t); i++)
sum=sum+M[x/i-1];
#pragma omp parallel for reduction (+:sum)
for (i=1; i<(long long)t; i++) {
sum=sum+(long long)M[i-1]*(x/(unsigned long long)i-x/(unsigned long long)(i+1));
}
return(sum);
}
// compute Mobius function
int newmobl(unsigned long long a, unsigned long long b, long long *out, unsigned int *table,
unsigned int tsize) {
unsigned int i,p,count;
unsigned long long beg,ps,rem;
long long t;
count=(unsigned int)(b-a);
for (i=0; i<count; i++)
out[i]=1;
for (i=0; i<tsize; i++) {
p=table[i];
ps=(unsigned long long)p*(unsigned long long)p;
if (ps>b)
goto askip;
rem=a-(a/ps)*ps;
if (rem!=0)
beg=ps-rem;
else
beg=0;
while (beg<count) {
out[beg]=0;
beg=beg+ps;
}
rem=a-(a/p)*p;
if (rem!=0)
beg=p-rem;
else
beg=0;
while (beg<count) {
out[beg]=-out[beg]*p;
beg=beg+p;
}
}
return(p);
askip:
for (i=0; i<count; i++) {
if (out[i]==0)
continue;
t=out[i];
if (t<0)
t=-t;
if ((unsigned long long)t<(i+a))
out[i]=-out[i];
}
for (i=0; i<count; i++) {
if (out[i]==0)
continue;
if (out[i]>0)
out[i]=1;
else
out[i]=-1;
}
return(1);
}
// compute primes
unsigned int primed(unsigned int *out, unsigned int tsize,
unsigned int *table,unsigned int limit) {
unsigned int d;
unsigned int i,j,k,l,flag,count;
count=tsize;
for (i=0; i<tsize; i++)
out[i]=table[i];
j=table[tsize-1]+1;
for (d=j; d<=limit; d++) {
if(d==(d/2)*2) continue;
if(d==(d/3)*3) continue;
if(d==(d/5)*5) continue;
if(d==(d/7)*7) continue;
if(d==(d/11)*11) continue;
l=(unsigned int)(10.0+sqrt((double)d));
k=0;
if (l>table[tsize-1])
return(0);
else {
for (i=0; i<tsize; i++) {
if (table[i]<=l)
k=i;
else
break;
}
}
flag=1;
l=k;
for (i=0; i<=l; i++) {
k=table[i];
if ((d/k)*k==d) {
flag=0;
break;
}
}
if (flag==1)
out[count]=d;
count=count+flag;
}
return(count);
}
unsigned long long Msize=15000000000; // for 64 GB of RAM, 64-bit OS
unsigned int tsize=1230; // prime look-up table size
unsigned int tmpsiz=10000; // used to compute Mobius function
unsigned int t2size=200000; // prime look-up table size
unsigned int Tsize=100000000; // used to compute partial sums of Mobius function
unsigned int oldnews=170000; // used to factor N
unsigned int prsize=1000; // used to factor N
//
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,*T,ltemp;
int *M;
unsigned int *ntable;
unsigned int *pritab;
unsigned int p,tindex,prind,delta,joff,newind,count,total,h,k,ntsize;
unsigned int g,f,L;
unsigned long long index,ta,tb,j,N,ut,pz,tz,mcount,start;
long long i;
int savet,t,ID,d,e;
FILE *Outfp;
Outfp = fopen("out1arq.dat","w");
omp_set_dynamic(0);
omp_set_num_threads(8);
#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=(long long*)malloc((Tsize+1)*8);
if (T==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<(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("computing j(x) \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;
}
L=0;
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=(N/(unsigned long long)e)/(unsigned long long)Msize+1;
ltemp=newmert(start,N/e,M);
if (ltemp==0x7fffffffffffffff) {
printf("error: s>(x/t) \n");
goto zskip;
}
T[e-1]=1-ltemp;
}
for (e=(int)(L/2); e>=1; e--) {
sum=0;
#pragma omp parallel for reduction (+:sum)
for (d=1; d<=(int)(L/(unsigned int)e-1); d++)
sum=sum+T[(d+1)*e-1];
T[e-1]=T[e-1]-sum;
}
sum=0;
savet=(int)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+T[e-1]*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 %d \n",N,sum,newind+1,savet,L);
fprintf(Outfp," %I64x, %I64x, %d, %d, \n",N,sum,newind+1,savet);
}
zskip:
fclose(Outfp);
return;
}