﻿ compute Dirichlet products
```/*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))
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
//
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)
if (subcur==3)
goto sskip;
if (subcur==4)
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);
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);
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;
}
```