/******************************************************************************
* *
* COMPUTE CONVOLUTION WITH MOBIUS FUNCTION *
* 11/04/18 (dkc) *
* *
* Normally distributed. *
* *
******************************************************************************/
#include <math.h>
#include <stdio.h>
#include <errno.h>
#include "zerob.h"
#include "table2.h"
extern char *malloc();
int mobius(unsigned int a, unsigned int *table, unsigned int tsize);
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);
unsigned int primed(unsigned int *out, unsigned int tsize,
unsigned int *table,unsigned int limit);
//
unsigned int MAXN=10000000;
unsigned int MINN=1;
unsigned int newcnt=10000000; // for curves and sub-curves ("check" not equal to 0)
unsigned int wrap=0;
unsigned int scale=3; // if set to 1, divide zeros by their logarithms (not used)
// if set to 2, multiply zeros by their logarithms
// if set to 3, use Odlyzko's metric
// if set to 4, use Oklyzko's metric and compute next-to-nearest neighbors
unsigned int diff=1; // if set, take differences
unsigned int inc=0; // offset for differences (0 for next)
//
unsigned int check=1; // check if N is prime, etc.
unsigned int p=3; // power of prime (up to 27 complete)
unsigned int pprod=1; // if set when p>1 and odd, do primes and combinations
unsigned int nope=0; // if set when p=3, primes are omitted even in skip[3]=0
// if set when check=0, values at prime x are omitted
unsigned int subcur=2; // sub-curves, set to 0 to do all (if elements in skip array set to 0)
// 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)
// 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 (includes 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,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35 for p=27
//
// check if skip[19] (p=3 to 11) or skip[18] (p=13 to 21) is zero, enables bypass
//
//unsigned int skip[35]={0,0,0,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}; // out307t.dat
//unsigned int skip[35]={0,0,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}; // out309at.dat
//unsigned int skip[35]={1,1,0,0,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}; // out309bt.dat
//unsigned int skip[35]={0,0,0,0,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}; // out311t.dat
//unsigned int skip[35]={0,0,0,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}; // out313bt.dat
//unsigned int skip[35]={1,1,1,0,0,0,0,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}; // out313at.dat
//unsigned int skip[35]={1,1,1,1,0,0,0,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}; // out315at.dat
unsigned int skip[35]={1,0,0,0,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}; // out315bt.dat
//unsigned int skip[35]={1,1,0,0,0,0,0,1,0,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}; // out317a.dat
//unsigned int skip[35]={1,1,1,0,1,1,1,0,0,0,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}; // 19 out319a.dat
//unsigned int skip[35]={1,1,1,1,1,1,0,1,1,1,1,1,0,0,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}; // 19 out319b.dat
//unsigned int skip[35]={1,1,1,1,0,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}; // 19 out319c.dat
//unsigned int skip[35]={1,0,1,1,1,1,1,1,1,0,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}; // 21 out321a.dat
//unsigned int skip[35]={0,1,0,1,1,1,0,1,0,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}; // 21 out321b.dat
//unsigned int skip[35]={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,1,1,1,1,1,1,1,1}; //
//unsigned int skip[35]={0,0,0,0,1,1,0,1,0,1,0,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}; // 19 out319p.dat
//unsigned int skip[35]={1,1,1,1,0,0,1,0,1,0,1,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}; // 19 out319q.dat
//unsigned int skip[35]={1,1,1,0,1,1,1,0,0,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}; // 19 out319y.dat
unsigned int tsize=17984;
unsigned int musig=0; // mean, standard deviation, N, or all three
void main () {
errno_t err;
unsigned int i,N,count,t1size,wcount;
unsigned int *table1;
double sum,total,mean,oldmean,m2;
double max,min,pi2;
double pi=3.14159265358979323846;
double *save;
FILE *Outfp;
err = fopen_s(&Outfp,"out3.dat","w");
if (Outfp==NULL)
return;
table1=(unsigned int*) malloc(16000008);
if (table1==NULL) {
printf("not enough memory \n");
return;
}
save=(double*) malloc(400000008);
if (save==NULL) {
printf("not enough memory \n");
return;
}
pi2=pi*2.0;
t1size=primed(table1,tsize,table,10000000); // make larger prime look-up table
printf("prime look-up table size=%d, last prime=%d \n",t1size,table1[t1size-1]);
for (i=0; i<MAXN; i++) {
if (riem[i+1]-riem[i]<=0.0) {
printf("error %d %llf %llf \n",i,riem[i+1],riem[i]);
return;
}
if (riem[i+1]-riem[i]>4.0) {
printf("big gap %d %llf %llf \n",i,riem[i+1],riem[i]);
}
}
//
if ((musig==0)&&(wrap==1)&&(diff==0))
fprintf(Outfp," %llf \n",1);
mean=0.0;
m2=0.0;
count=0;
wcount=1;
max=0.0;
min=1000000.0;
for (N=MINN; N<=MAXN; N++) {
if (check!=0) {
if (nucheck(N,p,pprod,subcur,t1size,table1,skip,nope)==1) {
count=count+1;
if (count>newcnt) {
count=count-1;
goto zskip;
}
}
else
continue;
}
else {
if (nope!=0) {
if (nucheck(N,1,0,subcur,t1size,table1,skip,nope)==0) {
count=count+1;
if (count>newcnt) {
count=count-1;
goto zskip;
}
}
else
continue;
}
else
count=count+1;
}
sum=0.0;
for (i=1; i<=N; i++) {
if (N==(N/i)*i) {
if (scale==0) {
if (diff==0)
sum=sum+riem[i-1]*mobius(i,table1,t1size);
else
sum=sum+(riem[i+inc]-riem[i-1])*mobius(i,table1,t1size);
}
else {
if (diff==0) {
if (scale==1)
sum=sum+riem[i-1]/log(riem[i-1])*mobius(i,table1,t1size);
else
sum=sum+riem[i-1]*log(riem[i-1])*mobius(i,table1,t1size);
}
else {
if (scale==1)
sum=sum+(riem[i+inc]/log(riem[i+inc])-riem[i-1]/log(riem[i-1]))*mobius(i,table1,t1size);
else {
if (scale==2)
sum=sum+(riem[i+inc]*log(riem[i+inc])-riem[i-1]*log(riem[i-1]))*mobius(i,table1,t1size);
else
sum=sum+(riem[i+inc]-riem[i-1])*log(riem[i-1]/pi2)/pi2*mobius(i,table1,t1size);
}
}
}
}
}
if (sum>max)
max=sum;
if (sum<min)
min=sum;
if (count==1) {
total=sum;
mean=sum;
oldmean=sum;
m2=sum;
}
else {
total=total+sum;
mean=total/(double)count;
m2=m2+(sum-oldmean)*(sum-mean);
oldmean=mean;
if (wcount==wrap) {
printf("N=%d mean=%llf standard deviation=%llf \n",N,mean,sqrt(m2/(double)(count-1)));
if (musig==0)
fprintf(Outfp,"%llf \n",mean);
else {
if (musig==1)
fprintf(Outfp,"%llf \n",sqrt(m2/(double)(count-1)));
else {
if (musig==2)
fprintf(Outfp," %d \n",N);
else
fprintf(Outfp," %llf, %llf, %llf, %llf, %llf, \n",(double)N,mean,sqrt(m2/(double)(count-1)),min,max);
}
}
wcount=1;
}
else
wcount=wcount+1;
}
if (wrap==0) {
if (count==(count/1000)*1000)
printf("sum=%llf, N=%d \n",sum,N);
if (scale!=4)
fprintf(Outfp," %.16f, \n",sum);
else
save[count-1]=sum;
}
}
zskip:
if (scale==4) {
for (i=0; i<(count-1); i++)
fprintf(Outfp," %llf, \n",save[i]+save[i+1]);
}
printf("max=%llf min=%llf count=%d \n",max,min,count);
printf("mean=%llf standard deviation=%llf \n",mean,sqrt(m2/(double)(count-1)));
fclose(Outfp);
return;
}