/******************************************************************************/
/*									      */
/*  COMPUTE A AND B VALUES						      */
/*  04/03/10 (dkc)							      */
/*									      */
/*  This C program computes a and b values.				      */
/*									      */
/******************************************************************************/
#include <stdio.h>
#include <math.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 n=128;			  // number of words
unsigned int i,j,k,a,b,length;
unsigned int M[256],O[256],X[256],W[256],Y[256],Z[256];
double w,x,y,z;
unsigned int v[178]={0,3,6,8,11,	  // I	(invalid lengths)
		    13,16,19,21,24,	  // I
		    26,29,32,34,37,	  // I
		    39,42,44,47,50,	  // J
		    52,55,57,60,63,	  // J
		    65,68,70,73,75,	  // K
		    78,81,83,86,88,	  // K
		    91,94,96,99,101,	  // K
		    104,106,109,112,114,  // L
		    117,119,122,125,127,  // L
		    130,132,135,138,140,  // L
		    143,145,148,150,153,  // M
		    156,158,161,163,166,  // M
		    171,174,176,179,	  // M
		    184,187,189,192,194,  // I'
		    197,200,202,205,207,  // I'
		    210,212,215,218,220,  // J'
		    223,225,228,231,233,  // J'
		    236,238,241,243,246,  // K'
		    249,251,254,256,259,  // K'
		    262,264,267,269,272,  // K'
		    275,277,280,282,285,  // K'
		    287,290,293,295,298,  // L'
		    300,303,306,308,311,  // L'
		    313,316,318,321,324,  // M'
		    326,329,331,334,337,  // M'
		    339,342,344,347,349,  // N
		    352,355,357,360,362,  // N
		    365,368,370,373,375,  // N
		    378,380,383,386,388,  // O
		    391,393,396,399,401,  // O
		    404,406,409,412,414,  // O
		    417,419,422,424,427,  // P
		    430,432,435,437,440,  // P
		    443,445,448,450,453,  // P
		    458,461,463,466};	  // I
FILE *Outfp;
Outfp = fopen("out13.dat","w");
a=0;
b=0;
setn(Y, 0, n);	    // Y=0
setn(O, 0, n);
O[0]=0x20000000;    // O=1
setn(X, 0, n);
X[0]=0x10000000;    // x=1/2
setn(M, 0, n);
M[0]=0x40000000;    // M=large number
j=0;
for (i=0; i<2582;  i++) {   // starts losing precision at i=2582
   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;
      }
   copyn(X, W, n);
   subn(O, W, n);		// 1-x
   if ((W[0]&0x80000000)!=0)
      subn(Y, W, n);		// abs(1-x)
   copyn(W, Z, n);
   subn(M, Z, n);		// minimum - abs(1-x)
   if ((Z[0]&0x80000000)==0)	// check if negative
      copyn(W, M, n);		// set minimum to abs(1-x)
   length=3*a+2*b;
   for (k=0; k<178; k++) {
      if (length==v[k]) {
	 printf("invalid length=%d, i=%d \n",length,i);
	 fprintf(Outfp,"invalid length=%d, i=%d \n",length,i);
	 goto zskip;
	 }
      }
   w=(double)a;
   w=w/(double)b;		// a/b
   x=(double)X[0];		// approximate x
   x=x/8192.0;
   x=x/65536.0;
   y=(double)W[0];		// approximate abs(1-x)
   y=y/8192.0;
   y=y/65536.0;
   printf("a=%d, b=%d, l=%d, x=%e, delta=%e, ratio=%e \n",a,b,length,x,y,w);
// fprintf(Outfp,"a=%d, b=%d, l=%d, x=%e, delta=%e, ratio=%e \n",a,b,length,x,y,w);
   fprintf(Outfp," %d, %d,",a,b);
// if (i>2570)
//    fprintf(Outfp,"  %#10x %#10x %#10x %#10x %#10x %#10x %#10x %#10x \n", X[120],X[121],X[122],X[123],X[124],X[125],X[126],X[127]);
   j=j+1;
   if (j==5) {
      fprintf(Outfp,"\n");
      j=0;
      }
   }
z=(double)M[0]; 		// approximate minimum
z=z/8192.0;
z=z/65536.0;
printf("min=%e \n",z);
//fprintf(Outfp,"min=%e \n",z);
zskip:
fclose(Outfp);
return(0);
}