/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C CHECK DIRICHLET PRODUCTS C C 01/30/18 (DKC) C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/ #include <stdio.h> #include <math.h> //#include "zeros1.h" // 20001 zeros of the Riemann zeta function //#include "rand.h" // 2001 uniformly distributed values #include "zerosdp.h" // Dirichlet products divided by their logarithms /****************************************************************************** * * * INTERPOLATE BETWEEN FRACTIONS * * 04/02/10 (dkc) * * * ******************************************************************************/ #include <math.h> double interp1(double ND, double ND1, unsigned int I, double init, double *out) { unsigned int i; double F,inc; inc=(ND1-ND)/(double)I; F=init; for (i=0; i<I; i++) { out[i]=F; F=F+inc; } return ND1; } /****************************************************************************** * * * CHECK IF PRIME, ETC. * * 01/30/18 (dkc) * * * ******************************************************************************/ #include <math.h> #include "table2.h" unsigned int pcheck(unsigned int N, unsigned int power, unsigned int extra, unsigned int subcur, unsigned int pfact, unsigned int tsize) { unsigned int i,j,k,l,p,q,r,s,pth,qth; if ((power==17)&&(extra!=0)) goto acskip; if ((power==15)&&(extra!=0)) goto gskip; if ((power==13)&&(extra!=0)) goto fskip; if ((power==11)&&(extra!=0)) goto eskip; if ((power==9)&&(extra!=0)) goto dskip; if ((power==7)&&(extra!=0)) goto cskip; if ((power==5)&&(extra!=0)) goto bskip; if ((power==3)&&(extra!=0)) goto askip; for (i=0; i<tsize; i++) { p=table[i]; pth=p; for (j=1; j<power; j++) pth=pth*p; if (N==pth) return(1); if (N<pth) break; } return(0); // // 3 // askip: if (pfact!=1) goto xskip; if (subcur==2) goto agskip; for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; if ((p*q)==N) return(1); if (N<(p*q)) break; } } if (subcur==1) return(0); agskip: for (i=0; i<tsize; i++) { p=table[i]; p=p*p*p; if (p==N) return(1); if (N<p) break; } return(0); xskip: // fixed p p=pfact; for (i=0; i<tsize; i++) { q=table[i]; if ((p*q)==N) return(1); if (N<(p*q)) break; } return(0); // // 5 // bskip: if (pfact!=1) goto wskip; if (subcur==2) goto ahskip; for (i=0; i<tsize; i++) { p=table[i]; p=p*p; for (j=i+1; j<tsize; j++) { q=table[j]; if ((p*q)==N) return(1); if (N<(p*q)) break; } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q; if ((p*q)==N) return(1); if (N<(p*q)) break; } } if (subcur==1) return(0); ahskip: for (i=0; i<tsize; i++) { p=table[i]; p=p*p*p*p*p; if (p==N) return(1); if (N<p) break; } return(0); wskip: // fixed p p=pfact; for (i=0; i<tsize; i++) { q=table[i]; q=q*q; if ((p*q)==N) return(1); if (N<(p*q)) break; } return(0); // // 7 // cskip: if (subcur==2) goto afskip; if (subcur==3) goto aiskip; for (i=0; i<tsize; i++) { p=table[i]; p=p*p*p; for (j=i+1; j<tsize; j++) { q=table[j]; if ((p*q)==N) return(1); if (N<(p*q)) break; } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q*q; if ((p*q)==N) return(1); if (N<(p*q)) break; } } if (subcur==1) return(0); afskip: for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } if (subcur==2) return(0); aiskip: for (i=0; i<tsize; i++) { p=table[i]; pth=p*p*p*p*p; if (N<pth) break; p=p*p*pth; if (p==N) return(1); if (N<p) break; } return(0); // // 9 // dskip: if (subcur==2) goto nskip; if (subcur==3) goto ajskip; for (i=0; i<tsize; i++) { p=table[i]; p=p*p*p*p; for (j=i+1; j<tsize; j++) { q=table[j]; if ((p*q)==N) return(1); if (N<(p*q)) break; } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q*q*q; if ((p*q)==N) return(1); if (N<(p*q)) break; } } if (subcur==1) return(0); nskip: for (i=0; i<tsize; i++) { p=table[i]; p=p*p; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q; if ((p*q)==N) return(1); if (N<(p*q)) break; } } if (subcur==2) return(0); ajskip: for (i=0; i<tsize; i++) { p=table[i]; pth=p*p*p*p*p; if (N<pth) break; p=p*p*p*p*pth; if (p==N) return(1); if (N<p) break; } return(0); // // 11 // eskip: if (subcur==4) goto akskip; if (subcur==3) goto pskip; if (subcur==2) goto oskip; for (i=0; i<tsize; i++) { p=table[i]; pth=p*p*p; if (N<pth) break; p=p*p*pth; if (N<p) break; for (j=i+1; j<tsize; j++) { q=table[j]; if ((p*q)==N) return(1); if (N<(p*q)) break; } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; qth=q*q*q; if (N<qth) break; q=q*q*qth; if (N<q) break; if ((p*q)==N) return(1); if (N<(p*q)) break; } } if (subcur==1) return(0); oskip: for (i=0; i<tsize; i++) { p=table[i]; p=p*p*p; if (N<p) break; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q; if ((p*q)==N) return(1); if (N<(p*q)) break; } } for (i=0; i<tsize; i++) { p=table[i]; p=p*p; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q*q; if (N<q) break; if ((p*q)==N) return(1); if (N<(p*q)) break; } } if (subcur==2) return(0); pskip: for (i=0; i<tsize; i++) { p=table[i]; p=p*p; for (j=i+1; j<tsize; j++) { q=table[j]; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; r=r*r; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } if (subcur==3) return(0); akskip: for (i=0; i<tsize; i++) { p=table[i]; pth=p*p*p*p; if (N<pth) break; pth=p*p*p*p*pth; if (N<pth) break; p=p*p*p*pth; if (p==N) return(1); if (N<p) break; } return(0); // // 13 // fskip: if (subcur==6) goto alskip; if (subcur==5) goto uskip; if (subcur==4) goto tskip; if (subcur==3) goto rskip; if (subcur==2) goto qskip; for (i=0; i<tsize; i++) { p=table[i]; pth=p*p*p; if (N<pth) break; p=p*p*p*pth; if (N<p) break; for (j=i+1; j<tsize; j++) { q=table[j]; if ((p*q)==N) return(1); if (N<(p*q)) break; } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; qth=q*q*q; if (N<qth) break; q=q*q*q*qth; if (N<q) break; if ((p*q)==N) return(1); if (N<(p*q)) break; } } if (subcur==1) return(0); qskip: for (i=0; i<tsize; i++) { p=table[i]; p=p*p*p*p; if (N<p) break; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q; if ((p*q)==N) return(1); if (N<(p*q)) break; } } for (i=0; i<tsize; i++) { p=table[i]; p=p*p; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q*q*q; if (N<q) break; if ((p*q)==N) return(1); if (N<(p*q)) break; } } if (subcur==2) return(0); rskip: for (i=0; i<tsize; i++) { p=table[i]; p=p*p*p; for (j=i+1; j<tsize; j++) { q=table[j]; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q*q; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; r=r*r*r; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } if (subcur==3) return(0); tskip: for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; if (N<(p*q*r)) break; for (l=k+1; l<tsize; l++) { s=table[l]; if ((p*q*r*s)==N) return(1); if (N<(p*q*r*s)) break; } } } } if (subcur==4) return(0); uskip: for (i=0; i<tsize; i++) { p=table[i]; p=p*p*p; if (N<p) break; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q*q; if (N<q) break; if ((p*q)==N) return(1); if (N<(p*q)) break; } } if (subcur==5) return(0); alskip: for (i=0; i<tsize; i++) { p=table[i]; pth=p*p*p*p*p; if (N<pth) break; pth=p*p*p*p*pth; if (N<pth) break; p=p*p*p*p*pth; if (p==N) return(1); if (N<p) break; } return(0); // // 15 // gskip: if (subcur==2) goto kskip; if (subcur==3) goto mskip; if (subcur==4) goto lskip; if (subcur==5) goto abskip; if (subcur==6) goto amskip; for (i=0; i<tsize; i++) { p=table[i]; pth=p*p*p*p; if (N<pth) break; p=p*p*p*pth; if (N<p) break; for (j=i+1; j<tsize; j++) { q=table[j]; if ((p*q)==N) return(1); if (N<(p*q)) break; } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; qth=q*q*q*q; if (N<qth) break; q=q*q*q*qth; if (N<q) break; if ((p*q)==N) return(1); if (N<(p*q)) break; } } if (subcur==1) return(0); kskip: for (i=0; i<tsize; i++) { p=table[i]; pth=p*p*p; if (N<pth) break; p=p*p*pth; if (N<p) break; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q; if ((p*q)==N) return(1); if (N<(p*q)) break; } } for (i=0; i<tsize; i++) { p=table[i]; p=p*p; for (j=i+1; j<tsize; j++) { q=table[j]; qth=q*q*q; if (N<qth) break; q=q*q*qth; if (N<q) break; if ((p*q)==N) return(1); if (N<(p*q)) break; } } if (subcur==2) return(0); mskip: for (i=0; i<tsize; i++) { p=table[i]; p=p*p; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; r=r*r; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } for (i=0; i<tsize; i++) { p=table[i]; p=p*p; for (j=i+1; j<tsize; j++) { q=table[j]; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; r=r*r; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } if (subcur==3) return(0); lskip: for (i=0; i<tsize; i++) { p=table[i]; p=p*p*p*p; for (j=i+1; j<tsize; j++) { q=table[j]; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q*q*q; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; r=r*r*r*r; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } if (subcur==4) return(0); abskip: for (i=0; i<tsize; i++) { p=table[i]; p=p*p*p*p; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q*q; if ((p*q)==N) return(1); if (N<(p*q)) break; } } for (i=0; i<tsize; i++) { p=table[i]; p=p*p*p; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q*q*q; if ((p*q)==N) return(1); if (N<(p*q)) break; } } if (subcur==5) return(0); amskip: for (i=0; i<tsize; i++) { p=table[i]; pth=p*p*p*p*p; if (N<pth) break; pth=p*p*p*p*p*pth; if (N<pth) break; p=p*p*p*p*p*pth; if (p==N) return(1); if (N<p) break; } return(0); // // 17 // acskip: if (subcur==2) goto adskip; if (subcur==3) goto sskip; if (subcur==4) goto aaskip; if (subcur==5) goto aeskip; if (subcur==6) goto anskip; for (i=0; i<tsize; i++) { p=table[i]; pth=p*p*p*p*p; if (N<pth) break; p=p*p*p*pth; if (N<p) break; for (j=i+1; j<tsize; j++) { q=table[j]; if ((p*q)==N) return(1); if (N<(p*q)) break; } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; qth=q*q*q*q*q; if (N<qth) break; q=q*q*q*qth; if (N<q) break; if ((p*q)==N) return(1); if (N<(p*q)) break; } } if (subcur==1) return(0); adskip: for (i=0; i<tsize; i++) { p=table[i]; pth=p*p*p*p; if (N<pth) break; p=p*p*pth; if (N<p) break; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q; if ((p*q)==N) return(1); if (N<(p*q)) break; } } for (i=0; i<tsize; i++) { p=table[i]; p=p*p; for (j=i+1; j<tsize; j++) { q=table[j]; qth=q*q*q*q; if (N<qth) break; q=q*q*qth; if (N<q) break; if ((p*q)==N) return(1); if (N<(p*q)) break; } } if (subcur==2) return(0); sskip: for (i=0; i<tsize; i++) { p=table[i]; p=p*p*p; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } for (i=0; i<tsize; i++) { p=table[i]; p=p*p*p; for (j=i+1; j<tsize; j++) { q=table[j]; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; r=r*r; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } for (i=0; i<tsize; i++) { p=table[i]; p=p*p; for (j=i+1; j<tsize; j++) { q=table[j]; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; r=r*r*r; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } for (i=0; i<tsize; i++) { p=table[i]; p=p*p; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q*q; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q*q; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; r=r*r; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; r=r*r*r; if ((p*q*r)==N) return(1); if (N<(p*q*r)) break; } } } if (subcur==3) return(0); aaskip: for (i=0; i<tsize; i++) { p=table[i]; p=p*p; for (j=i+1; j<tsize; j++) { q=table[j]; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; if (N<(p*q*r)) break; for (l=k+1; l<tsize; l++) { s=table[l]; if ((p*q*r*s)==N) return(1); if (N<(p*q*r*s)) break; } } } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; if (N<(p*q*r)) break; for (l=k+1; l<tsize; l++) { s=table[l]; if ((p*q*r*s)==N) return(1); if (N<(p*q*r*s)) break; } } } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; r=r*r; if (N<(p*q*r)) break; for (l=k+1; l<tsize; l++) { s=table[l]; if ((p*q*r*s)==N) return(1); if (N<(p*q*r*s)) break; } } } } for (i=0; i<tsize; i++) { p=table[i]; for (j=i+1; j<tsize; j++) { q=table[j]; if (N<(p*q)) break; for (k=j+1; k<tsize; k++) { r=table[k]; if (N<(p*q*r)) break; for (l=k+1; l<tsize; l++) { s=table[l]; s=s*s; if ((p*q*r*s)==N) return(1); if (N<(p*q*r*s)) break; } } } } if (subcur==4) return(0); aeskip: for (i=0; i<tsize; i++) { p=table[i]; pth=p*p*p; if (N<pth) break; p=p*p*pth; if (N<p) break; for (j=i+1; j<tsize; j++) { q=table[j]; q=q*q*q; if ((p*q)==N) return(1); if (N<(p*q)) break; } } for (i=0; i<tsize; i++) { p=table[i]; p=p*p*p; for (j=i+1; j<tsize; j++) { q=table[j]; qth=q*q*q; if (N<qth) break; q=q*q*qth; if (N<q) break; if ((p*q)==N) return(1); if (N<(p*q)) break; } } if (subcur==5) return(0); anskip: for (i=0; i<tsize; i++) { p=table[i]; pth=p*p*p*p*p; if (N<pth) break; pth=p*p*p*p*pth; if (N<pth) break; pth=p*p*p*p*pth; if (N<pth) break; p=p*p*p*p*pth; if (p==N) return(1); if (N<p) break; } return(0); } /*****************************************************************************/ /* */ /* MAIN PROGRAM */ /* 02/03/18 (dkc) */ /* */ /*****************************************************************************/ #include <stdio.h> #include <math.h> extern char *malloc(); unsigned int pcheck(unsigned int N, unsigned int p, unsigned int pprod, unsigned int subcur, unsigned int pfact, unsigned int size); double interp1(double ND, double ND1, unsigned int I, double init, double *out); // unsigned int MAXN=10000; unsigned int I=1; // number of values to interpolate unsigned int numlim=10000; unsigned int check=0; // check if N is prime, etc. unsigned int nolog=1; // if set, don't take logarithm unsigned int p=1; // power of prime (up to 17 complete) unsigned int pprod=1; // if set when p>1 and odd, do primes and combinations unsigned int norm=1; // if set, normalize unsigned int pfact=1; // if greater than 1, do pfact times primes or pfact // times square of primes (p=2 or 3, pprod=1) unsigned int subcur=0; // sub-curves, set to 0 to do all // 0, 1, 2 for p=3 // 0, 1, 2 for p=5 // 0, 1, 2, 3 for p=7 // 0, 1, 2, 3 for p=9 // 0, 1, 2, 3, 4 for p=11 // 0, 1, 2, 3, 4, 5, 6 for p=13 // 0, 1, 2, 3, 4, 5, 6 for p=15 // 0, 1, 2, 3, 4, 5, 6 for p=17; unsigned int sort=1; // sort zeros (for zerosdp) unsigned int limit=10001; // for sorting unsigned int scale=0; // if set, multiply zeros by their logarithms unsigned int outscale=0; // if set, divide output by logarithms unsigned int tsize=17984; // void main() { unsigned int i,N,count,first; double *output; double init,sum,savsum; double min,temp; FILE *Outfp; Outfp = fopen("out2.dat","w"); output=(double*) malloc(100008); if (output==NULL) { printf("not enough memory \n"); return; } if (sort!=0) { for (count=0; count<limit; count++) { N=count; min=riem[count]; for (i=count+1; i<limit; i++) { if (riem[i]<min) { min=riem[i]; N=i; } } temp=riem[count]; riem[count]=riem[N]; riem[N]=temp; } printf("input sorted \n"); } if (scale!=0) { for (i=0; i<10001; i++) riem[i]=riem[i]*log(riem[i]); } if (nolog==0) { init=log(riem[0]); for (i=0; i<(numlim+1); i++) init=interp1(log(riem[i]),log(riem[i+1]),I,init,&output[i*I]); } else { init=riem[0]; for (i=0; i<(numlim+1); i++) init=interp1(riem[i],riem[i+1],I,init,&output[i*I]); } first=1; count=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 (pcheck(N,p,pprod,subcur,pfact,tsize)==1) { if (first!=0) { first=0; savsum=sum-output[0]; } count=count+1; printf("i=%d %f %f \n",N,output[N-1],sum-output[0]); if (norm==0) fprintf(Outfp," %llf, \n",sum-output[0]); else fprintf(Outfp," %llf, \n",sum-output[0]-savsum); } } else { printf("i=%d %f %f \n",N,output[N-1],sum-output[0]); if (outscale==0) fprintf(Outfp," %llf, \n",sum); else fprintf(Outfp," %llf, \n",sum/log(sum)); } } printf("count=%d \n",count); fclose(Outfp); return; }