/******************************************************************************
*									      *
*  FACTOR								      *
*  01/30/18 (dkc)							      *
*									      *
******************************************************************************/
#include <math.h>
#include "table2.h"
unsigned int factor(unsigned int N, unsigned int *pfact, unsigned int tsize) {
unsigned int i,p,index,limit;
limit=(unsigned int)sqrt((double)N);
limit=limit+10;
index=0;
for (i=0; i<tsize; i++) {
   p=table[i];
   while (N==((N/p)*p)) {
      pfact[index]=p;
      index=index+1;
      N=N/p;
      }
   if (N==1)
      break;
   }
return(index);
}
/*****************************************************************************/
//									     //
// HISTOGRAM SPACINGS OF ZEROS AND DIRICHLET PRODUCTS OF ZEROS		     //
// 02/04/18 (dkc)							     //
//									     //
/*****************************************************************************/
#include <math.h>
#include <stdio.h>
#include "zerosi.h"     // 100000 zeros
//#include "rand.h"     // standard uniform distribution
//#include "zerosdpb.h"  // Dirichlet products, divided by their logarithms
unsigned int factor(unsigned int N, unsigned int *tfact, unsigned int tsize);
//
long double delta=2.0;	 // upper bound of spacings
long double hisfac=100;   // to make integers (for zerosdp, 200 bins)
unsigned int numzeros=100000; // number of zeros or Dirichlet products
unsigned int scale=0;  // normally set to 0, used for zerosi
unsigned int sort=0;   // if set, sort (for zerosdp)
unsigned int compl=0;	// normally set to 0, used for zerosdpb
unsigned int limit=100000;  // for zerosdpb
unsigned int tsize=17984;  // prime look-up table size
unsigned int outhis=1;   // if set, output histogram only
//
void main () {
unsigned int i,j,k,N,count;
long double temp,min;
unsigned int tfact[1000];
unsigned int histo[1000];
FILE *Outfp;
Outfp = fopen("out5.dat","w");
if (Outfp==NULL)
   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=1; i<numzeros; i++)
    riem[i]=riem[i]/log(riem[i]);
  }
for (i=0; i<1000; i++)
   histo[i]=0;
count=0;
for (i=1; i<numzeros; i++) {
   temp=riem[i]-riem[i-1];
   if (temp<0.0)
      temp=-temp;
   if (temp<=delta) {
      count=count+1;
      j=(unsigned int)(temp*hisfac);
      histo[j]=histo[j]+1;
      if (outhis==0) {
	 printf("%d %f %f %f \n",i+1,riem[i],riem[i-1],temp);
	 fprintf(Outfp,"%d %f %f %f \n",i+1,riem[i],riem[i-1],temp);
	 }
      if (outhis==0) {
	 j=factor(i,tfact,tsize);
	 for (k=0; k<j; k++) {
	    printf("%d ",tfact[k]);
	    fprintf(Outfp,"%d ",tfact[k]);
	    }
	 printf("\n");
	 fprintf(Outfp,"\n");
	 j=factor(i+1,tfact,tsize);
	 for (k=0; k<j; k++) {
	    printf("%d ",tfact[k]);
	    fprintf(Outfp,"%d ",tfact[k]);
	    }
	 printf("\n");
	 fprintf(Outfp,"\n");
	 }
      }
   }
printf("count=%d \n",count);
fprintf(Outfp,"count=%d \n",count);
count=0;
if (compl!=0) {
  for (i=300; i>0; i--) {
    printf(" %d \n",histo[i]);
    fprintf(Outfp," %d \n",histo[i]);
    }
  }
for (i=0; i<300; i++) {
  printf(" %d \n",histo[i]);
  fprintf(Outfp," %d \n",histo[i]);
  if (histo[i]!=0)
     count=count+1;
  }
printf("count=%d \n",count);
fclose(Outfp);
 return;
}