/******************************************************************************/
/*									      */
/*  COMPUTE UPPER BOUNDS OF MINIMA					      */
/*  04/16/10 (dkc)							      */
/*									      */
/*  This C program computes upper bounds of minima where n and n+2 are prime. */
/*									      */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
#include "table1.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=1;		// form flag
unsigned int n=1024;		// number of words
unsigned int i,j,k,a,b;
unsigned int O[1024],X[1024],Z[1024];
double x,l2,l3;
int m;
FILE *Outfp;
Outfp = fopen("out14f.dat","w");
l2=log(2);			// log(2)
l3=log(3);			// log(3)
printf("log2=%e, log3=%e \n",l2,l3);
a=0;
b=0;
setn(O, 0, n);
O[0]=0x20000000;		// O=1
setn(X, 0, n);
X[0]=0x10000000;		// x=1/2
for (i=0; i<20671;  i++) {	// starts losing precision at i=20671
   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[2760]) {
      printf("prime look-up table too small \n");
      goto zskip;
      }
   for (k=0; k<2761; k++) {	// check if n, n+2 are primes
      if ((table[k]==j)&&(table[k+1]==(j+2)))
	 goto askip;
      if (table[k]>j)
	 goto bskip;
      }
   goto bskip;
askip:
   if (flag==1) {		// check if 12 divides n+7
      if ((j+7)!=((j+7)/12)*12)
	 continue;
      }
   else {
      if ((j+7)==((j+7)/12)*12)
	 continue;
      }
   fprintf(Outfp," %d, %d,\n",j,m);
   printf("n=%d, x=%d \n",j,m);
bskip:
   j=0;
   }
zskip:
fclose(Outfp);
return(0);
}