/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C COMPUTE DIRICHLET PRODUCTS AND MOBIUS INVERSES C
C 05/20/19 (DKC) C
C C
C Set "mobinv" to 1, "check" to 0, "outind" to 0" to do Mobius inverse. C
C If extracting a sub-curve, set "extra" to 1 and "p" and "subcur" to C
C appropriate values. Set appropriate "norm" and "factor" values ("factor" C
C depends on "MAXN"). C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
#include <stdio.h>
#include <math.h>
#include "table2.h" // prime look-up table
//#include "gue45.h" // Gaussian ensemble (4500x4500 matrix)
#include "zero1.h" // 100000 zeros (more precision)
//#include "rand.h" // 2001 uniformly distributed values
extern char *malloc();
unsigned int nucheck(unsigned int N, unsigned int p, unsigned int pprod,
unsigned int subcur, unsigned int size, unsigned int *table,
unsigned int *skip, unsigned int nope);
double interp1(double ND, double ND1, unsigned int I,
double init, double *out);
unsigned int primed(unsigned int *out, unsigned int tsize,
unsigned int *table,unsigned int limit);
int mobius(unsigned int a, unsigned int *table1, unsigned int t1size);
//
unsigned int gue=0; // if set, use Gaussian ensemble
unsigned int MAXN=64;
unsigned int I=2; // number of values to interpolate
unsigned int numlim=32;
unsigned int scale=0; // if set to 1, take logarithms of zeros // if set to 2, take logarithms of zeros
unsigned int nonorm=1; // if (set, don't normalize
unsigned int check=0; // check if N is prime, etc.
unsigned int p=3; // power of prime (up to 17 complete)
unsigned int pprod=1; // if set when p>1 and odd, do primes and combinations
unsigned int subcur=2; // sub-curves, set to 0 to do all
unsigned int nope=0; // if set when p=3, omit primes even if skip[3]==0
unsigned int outind=0; // if set, output indices too
unsigned int mobinv=1; // if set, do Mobius inversion
unsigned int comp=2; // if set to 1, take differences
// if set to 2, compute relative error
unsigned int clip=1; // if set only output upper half (comp=2 only)
unsigned int extra=0; // if set, extract sub-curve (comp=0 only)
double norm=31.713108; // I=2 (extra=0 only)
//double norm=80.453014; // I=2, n=pq (3.2)
//double norm=250.147742; // I=2, n=pqr (7.3)
//double norm=1101.339335; // I=5, p^3q^2r (17.3)
//double norm=298.432389; // I=2, n=p^2q^2 (9.3)
//double norm=1061.644516; // I=2, n=p^3q^3 (13.6)
//double norm=-1505.145075; // I=4, n=pqr (eigenvalues)
//double norm=-754.461338; // I=2, n=pq (eigenvalues)
//
double factor=1.66103262; // 64 I=2
//double factor=1.71146907; // 128
//double factor=1.73182393; // 256
//double factor=1.78701587; // 512
//double factor=1.78691539; // 1024
//double factor=1.80503140; // 2048
//double factor=1.82016313; // 4096
//double factor=1.82643984; // 8192
//double factor=1.84477360; // 16384
//double factor=1.85090537; // 32768
//double factor=1.86106885; // 65536
//double factor=1.87007977; // 131072
//double factor=1.87746605; // 262144
//double factor=1.88544670; // 524288
//double factor=1.89075898; // 1048576
//double factor=1.89658165; // 2097152
//double factor=1.90156268; // 4194304
//
//double factor=1.69176559; // 128 I=2, n=pq
//double factor=1.71556600; // 256
//double factor=1.73199042; // 512
//double factor=1.78817824; // 1024
//double factor=1.78893021; // 2048
//double factor=1.80629720; // 4096
//double factor=1.82068655; // 8192
//double factor=1.82661960; // 16384
//double factor=1.84565112; // 32768
//double factor=1.85122113; // 65536
//double factor=1.86125581; // 131072
//double factor=1.87018543; // 262144
//double factor=1.87739527; // 524288
//double factor=1.88581843; // 1048576
//double factor=1.89089448; // 2097152
//double factor=1.8970613867; // 4194304
//
//double factor=1.57416755; // 256 I=2, n=pqr
//double factor= 1.78315930; // 512
//double factor= 1.72542247; // 1024
//double factor= 1.74374660; // 2048
//double factor= 1.79300287; // 4096
//double factor= 1.78852346; // 8192
//double factor= 1.81300810; // 16348
//double factor= 1.82412852; // 32768
//double factor= 1.84018826; // 65536
//double factor= 1.84683077; // 131072
//double factor= 1.85996662; // 262144
//double factor= 1.86524311; // 524288
//double factor= 1.87329003; // 1048576
//double factor=1.8824472896; // 2097152
//double factor=1.8875485468; // 4194304
//
//double factor=3.6791295316; // 20000 I=5 n=p^3q^2r
//double factor=3.8493129276; // 100000
//double factor=4.03015276; // 500000
//double factor=4.17672296; // 2500000
//
//double factor=1.82575404; // 8192 n=p^2q^2
//double factor=1.81291183; // 16384
//double factor=1.81871522; // 32768
//double factor=1.67551833; // 65536
//double factor=1.75388948; // 131072
//double factor=1.82191754; // 262144
//double factor=1.80401095; // 524288
//double factor=1.87108078; // 1048576
//double factor=1.81071449; // 2097152
//double factor=1.8255899111; // 4194304
//
//double factor=1.64990162; // 262144, n=p^3q^3
//double factor=1.65460131; // 524288
//double factor=2.30200446; // 1048576
//double factor=1.63679995; // 2097152
//double factor=1.9350651633; // 4194304
//
//double factor=2.5989836788; // 1000, n=pqr (eigenvalues)
//double factor=1.61998712996; // 1000, n=pq (eigenvalues)
// 0, 1, 2, 3, 4 for p=3 (includes p^2 and p)
// 0, 1, 2, 3 for p=5 (includes p^4)
// 0, 1, 2, 3, 4 for p=7 (includes p^6)
// 0, 1, 2, 3, 4 for p=9 (includes p^8 [use])
// 0, 1, 2, 3, 4, 5 for p=11 (includes p^10)
// 0, 1, 2, 3, 4, 5, 6, 7, 8 for p=13 (includes p^12)
// 0, 1, 2, 3, 4, 5, 6, 7, 8 for p=15 (includes p^14)
// 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 for p=17 (includes p^16)
// 0, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18 for p=19 (also p^18)
// 0,1,2,3,4,5,6,7,8,9,10,11,12 for p=21 (includes p^20)
// 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22 for p=23
// 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27 for p=25
// 0,1,2,3, ..., 35 for p=27
//
// make sure to check skip[19], should be zero for p=3 only (for bypass)
//
//unsigned int skip[20]={1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}; // p=8
//unsigned int skip[20]={0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}; // p=9
//unsigned int skip[20]={1,1,1,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1}; // p=13
//unsigned int skip[20]={1,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1}; // p=15
//unsigned int skip[20]={1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1}; // p=15
//unsigned int skip[20]={1,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1}; // p=17
//unsigned int skip[20]={1,1,1,0,1,1,1,0,0,0,1,1,1,1,1,1,0,1,1,1}; // p=19
//unsigned int skip[20]={1,1,1,1,1,1,0,1,1,1,1,1,0,0,1,0,1,1,1,1}; // p=19
//unsigned int skip[20]={1,1,0,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1}; // p=19
//unsigned int skip[20]={1,1,1,1,0,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1}; // p=19
//unsigned int skip[20]={1,0,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1}; // p=19
//unsigned int skip[20]={1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1}; // p=19
//unsigned int skip[20]={0,1,0,1,1,1,0,1,0,1,1,0,0,0,0,0,0,0,0,1}; // p=21
//unsigned int skip[20]={1,0,1,1,1,1,1,1,1,0,1,0,0,0,0,0,0,0,0,1}; // p=21
//unsigned int skip[20]={1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1};
//unsigned int skip[27]={1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}; // p=23
//unsigned int skip[27]={1,1,1,1,1,1,0,0,1,1,1,1,0,0,0,1,1,1,0,1,1,1,1,1,1,1,1};
//unsigned int skip[27]={1,1,1,0,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,0,1,1,1,0}; // #25d
//unsigned int skip[27]={1,1,1,1,0,1,1,1,1,0,0,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1}; // #25e
unsigned int skip[35]={1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1};
//
unsigned int rnd=0; // if set, scale up the standard uniform distribution
double rndfact=5600.0; // factor
unsigned int tsize=17984; // prime look-up table size
//
void main() {
unsigned int i,N,count,t1size;
unsigned int *table1;
double *output,*dout;
double init,sum,max;
double min,temp;
FILE *Outfp;
Outfp = fopen("out2.dat","w");
if ((check!=0)&&(extra!=0)) {
printf("check and extra both set \n");
return;
}
output=(double*) malloc(800008);
if (output==NULL) {
printf("not enough memory \n");
return;
}
dout=(double*) malloc(800008);
if (dout==NULL) {
printf("not enough memory \n");
return;
}
if (output==NULL) {
printf("not enough memory \n");
return;
}
table1=(unsigned int*) malloc(800008);
if (table1==NULL) {
printf("not enough memory \n");
return;
}
t1size=primed(table1,tsize,table,100000); // make larger prime look-up table
printf("prime look-up table size=%d, last prime=%d \n",t1size,table1[t1size-1]);
if (rnd!=0) {
for (i=0; i<MAXN; i++)
riem[i]=riem[i]*rndfact;
}
if (scale==1) {
for (i=0; i<MAXN; i++)
riem[i]=log(riem[i]);
}
init=riem[0];
for (i=0; i<(numlim+1); i++)
init=interp1(riem[i],riem[i+1],I,init,&output[i*I]);
//
count=0;
min=10000000;
max=0.0;
for (N=1; N<=MAXN; N++) {
sum=0.0;
for (i=1; i<=N; i++) {
if (N==(N/i)*i)
sum=sum+output[N/i-1];
}
if (check!=0) {
if (nucheck(N,p,pprod,subcur,t1size,table1,skip,nope)==1) {
count=count+1;
if (nonorm!=0)
temp=sum;
else
temp=sum-output[0];
printf("N=%d %llf \n",N,temp);
if (temp>max)
max=temp;
if (temp<min)
min=temp;
if (outind==0)
fprintf(Outfp," %llf, \n",temp);
else
fprintf(Outfp," %llf, %d, \n",temp,N);
}
}
else {
if (N==(N/1000)*1000)
printf("N=%d %llf \n",N,sum);
if (mobinv==0)
fprintf(Outfp," %llf, \n",sum);
else
dout[N-1]=sum;
}
}
if (check!=0)
printf("count=%d, max=%llf, min=%llf \n",count,max,min);
if ((mobinv==0)||(check!=0))
goto zskip;
//
// Mobius inversion
//
printf("Mobius inversion \n");
//
// scale
//
if (gue==0)
dout[0]=14.13472514173470;
else
dout[0]=-189.1347051; // eigenvalue
for (i=1; i<MAXN; i++)
dout[i]=(dout[i]-norm)*factor+norm;
for (N=1; N<=MAXN; N++) {
sum=0.0;
for (i=1; i<=N; i++) {
if (N==(N/i)*i)
sum=sum+dout[i-1]*mobius(N/i,table1,t1size);
}
if (N==(N/1000)*1000)
printf("N=%d %llf \n",N,sum);
if (comp==0) {
if (extra==0)
fprintf(Outfp," %llf, \n",sum);
else
output[N-1]=sum;
}
else {
if (comp==1)
fprintf(Outfp," %llf, \n",sum-riem[N-1]);
else {
if (clip==0)
fprintf(Outfp," %llf, \n",sum/riem[N-1]-1.0);
else {
if (N>(MAXN/I))
fprintf(Outfp," %llf, \n",sum/riem[N-1]-1.0);
}
}
}
}
if ((extra!=0)&&(comp==0)) {
count=0;
for (N=1; N<=MAXN; N++) {
if (nucheck(N,p,pprod,subcur,t1size,table1,skip,nope)==1) {
count=count+1;
fprintf(Outfp," %llf, \n",output[N-1]);
}
}
printf("count=%d \n",count);
}
zskip:
fclose(Outfp);
return;
}