/*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;
}