﻿ compute j(x)
```/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  COMPUTE MERTENS FUNCTION (sum of M(x/i)^2 where i|x, find local maxima)    C
C  10/30/15 (DKC) (maxima based on j(x))				      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);
unsigned long long MAXN=3500000000; // (2^31-1=4294967295)
unsigned long long BEGINN=2;  // beginning N
unsigned long long Msize=3500000000;   // for 16 GB of RAM, 64-bit OS
//
void main() {
int t,*temp;
int *M;
unsigned int total,count,tindex,p,h;
unsigned int pritab[1000],prind,delta,joff,newind;
unsigned int j,k;
int savet;
unsigned long long i,tz,pz;
unsigned long long *oldtmp,*newtmp;
unsigned long long N,ta,tb,index;
unsigned long long ut;
double sum,maxt;
FILE *Outfp;
Outfp = fopen("out1ap.dat","w");
temp=(int *) malloc(4004);
if (temp==NULL) {
printf("not enough memory \n");
return;
}
oldtmp=(unsigned long long*) malloc(64000);
if (oldtmp==NULL) {
printf("not enough memory \n");
return;
}
newtmp=(unsigned 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;
}
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");
for (i=1; i<=Msize; i++)
M[i]=M[i-1]+M[i];
//
printf("finding maxima \n");
maxt=0.0;
for (N=BEGINN; N<=MAXN; N++) {
if ((N>13)&&((N&1)==1))     // Note!  Maxima appear to be even for N>13.
continue;
//
// factor N
//
ut=N;
prind=0;
tindex=0;
total=0;
h=(unsigned int)(sqrt((double)ut)+0.001); // For p>2, the difference between
for (i=0; i<41538; i++) {		     // successive primes is greater
p=table[i];			     // than 1, so adding 0.001 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)
}
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;
}
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:
sum=1.0;
for (i=0; i<newind; i++) {
tz=oldtmp[i];
t=M[tz-1];
if (tz==N)
savet=t;
sum=sum+(double)t*(double)t;
}
if (sum>maxt) {
maxt=sum;
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;
}
```