/******************************************************************************/
/*									      */
/*  COMPUTE UPPER BOUNDS OF MINIMA					      */
/*  04/19/10 (dkc)							      */
/*									      */
/*  This C program computes upper bounds of minima where n is prime.	      */
/*									      */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
#include "table2.h"

void addn(unsigned int *a, unsigned int *b, unsigned int n);
void subn(unsigned int *c, unsigned int *d, unsigned int n);
void shiftn(unsigned int *a, unsigned int *b, unsigned int shift, unsigned int n);
void copyn(unsigned int *a, unsigned int *b, unsigned int n);
void setn(unsigned int *a, unsigned int b, unsigned int n);
int main () {
unsigned int flag=0;		// form flag (0, 1, 2, or 3)
unsigned int n=2048;		// number of words
unsigned int i,j,k,a,b,count;
unsigned int O[2048],X[2048],Z[2048];
unsigned int histo[100];	 // histogram
double x,l2,l3;
int m,sum,lastp,lastm,delta,savep,iters,oldcount;
FILE *Outfp;
Outfp = fopen("out14g.dat","w");
for (i=0; i<100; i++)		// clear histogram
   histo[i]=0;
sum=0;
iters=0;			// period count
lastp=0;			// last prime
lastm=0;			// last minimum
savep=0;			// prime from last period
count=0;			// prime count
oldcount=1;			// old prime count
l2=log(2);			// log(2)
l3=log(3);			// log(3)
a=0;				// clear a
b=0;				// clear b
setn(O, 0, n);
O[0]=0x20000000;		// O=1
setn(X, 0, n);
X[0]=0x10000000;		// x=1/2
for (i=0; i<41346; i++) {	// starts losing precision at i=41346
   copyn(O, Z, n);
   subn(X, Z, n);		// x-1
   if ((Z[0]&0x80000000)==0) {	// x>1
      copyn(X, Z, n);
      addn(X, X, n);
      addn(Z, X, n);
      shiftn(X, X, 2, n);	// x=(3/4)x
      a=a+1;
      }
   else {			// x<1
      copyn(X, Z, n);
      addn(X, X, n);
      addn(Z, X, n);
      shiftn(X, X, 1, n);	// x=(3/2)x
      b=b+1;
      }
   x=(double)(a+1)/((double)(2*a+b+1)*l2-(double)(a+b)*l3);
   m=(int)x;			// upper bound of minimum
   j=i+1;			// n
   if (j>table[5131]) {
      printf("prime look-up table too small \n");
      goto zskip;
      }
   for (k=0; k<5132; k++) {	// check if n is prime
      if (table[k]==j)
	 goto askip;
      if (table[k]>j)
	 goto bskip;
      }
   goto bskip;
askip:
   if (flag==0) {		// check if 12 divides n+1
      if ((j+1)!=((j+1)/12)*12)
	 continue;
      else
	 count=count+1;
      }
   if (flag==1) {		// check if 12 divides n+5
      if ((j+5)!=((j+5)/12)*12)
	 continue;
      else
	 count=count+1;
      }
   if (flag==2) {		// check if 12 divides n+7
      if ((j+7)!=((j+7)/12)*12)
	 continue;
      else
	 count=count+1;
      }
   if (flag==3) {		// check if 12 divides n+11
      if ((j+11)!=((j+11)/12)*12)
	 continue;
      else
	 count=count+1;
      }
   if ((lastm>0)&&(m<0)) {	// check for new period
      delta=lastp-savep;	// difference in n values
      savep=lastp;		// save n value
      if (iters!=0) {
	 sum=sum+delta; 	// sum of differences
	 if ((delta/12)<100)	// histogram differences
	    histo[delta/12]+=1;
	 else {
	    printf("histogram too small \n");
	    goto zskip;
	    }
	 printf("period=%d, n=%d, minimum=%d, delta=%d \n",count-oldcount,lastp,lastm,delta/12);
	 }
      oldcount=count;		// save prime count
      iters=iters+1;		// increment period count
      }
   lastp=j;			// save n
   lastm=m;			// save minimum
   fprintf(Outfp," %d %d %d \n",count,j,m);
bskip:
   j=0;
   }
printf("\n");
printf("sum=%d, count=%d, average=%f \n",sum,iters-1,(float)sum/(float)(iters-1));
//
// histogram of differences
//
printf("\n");
for (i=20; i<80; i++)
   printf("i=%d, h[i]=%d \n",i,histo[i]);
zskip:
fclose(Outfp);
return(0);
}