/* COMPUTE GENERALIZED CONTINUED-FRACTION CONVERGENTS OF LOG(3)/LOG(2)	     */
/* 03/09/13 (dkc)							     */
/*									     */
/* This C program approximates given real numbers with rational numbers using*/
/* the continued fraction algorithm.					     */
/*									     */
/*****************************************************************************/
#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 copyn(unsigned int *a, unsigned int *b, unsigned int n);
void setn(unsigned int *a, unsigned int b, unsigned int n);

unsigned int L2[4]={	// log(2)
 0x162e42fe,
 0xfa39ef35,
 0x793c7673,
 0x007e5ed5};
unsigned int L3[4]={   // log(3)
 0x2327d4f5,
 0x5a06152e,
 0xd48331aa,
 0xa0a76f96};

void main() {
unsigned int m=4;
unsigned int T[4],U[4];
unsigned int i,j;
double X,Y,B,MAX;
unsigned int I,J,k,temph,tempk;
unsigned int A[51],AP[51],H[51],K[51];
FILE *Outfp;
Outfp = fopen("expand3.dat","w");
//
// COMPUTE THE PARTIAL QUOTIENTS
//
I=1;
A[1]=1;
printf("accurate partial quotients \n");
printf("j=%d \n",A[1]);
copyn(L2, T, m);
subn(L3, L2, m);  // B=X-(double)A[1]
copyn(T, L3, m);  // X=1.0/B
//printf("L3=%#010x, L2=%#010x \n",L3[0],L2[0]);
for (I=2; I<=20; I++) {
   j=0;
   copyn(L2, T, m);
   U[0]=0;
   while ((U[0]&0x80000000)!=0x80000000) {  // A[i]=(unsigned int)X
      j=j+1;
      addn(L2, T, m);	// L2*(j+1)
      copyn(T, U, m);
      subn(L3, U, m);
      }
   A[I]=j;
   printf("j=%d \n",j);
   copyn(L2, U, m);
   subn(T, L2, m);   // L2*j
   subn(L3, L2, m);  // B=X-(double)A[1]
   copyn(U, L3, m);  // X=1.0/B
// printf("L3=%#010x, L2=%#010x \n",L3[0],L2[0]);
   }
i=20;
printf("\n");
//
      printf("CONTINUED FRACTION EXPANSION OF A REAL NUMBER \n");
      MAX=65536.0;
      MAX=MAX*32768.0-1.0;
//
// REAL NUMBER
//
      X=log(3)/log(2);
      printf("X=%f \n",X);
//
// COMPUTE THE PARTIAL QUOTIENTS
//
      I=1;
      AP[1]=(unsigned int)X;
      printf("partial quotients computed using floating point arithmetic \n");
      printf("A=%d \n",AP[1]);
      if((X<0.0)&&(X!=(double)AP[1])) AP[1]=AP[1]-1;
      B=X-(double)AP[1];
      if(B<=0.0) goto L300;
      X=1.0/B;
      for (I=2; I<=20; I++) {
	 if(X>MAX) goto L250;
	 AP[I]=(unsigned int)X;
	 B=X-(double)AP[I];
//	 printf("A=%d, X=%e, B=%e \n",AP[I],X,B);
	 printf("A=%d \n",AP[I]);
	 if(B<=0.0) goto L300;
	 X=1.0/B;
	 }
      I=21;
L250: I=I-1;
//
// COMPUTE THE CONVERGENTS
//
L300: printf("\n");
      J=i;
      H[1]=A[1];
      K[1]=1;
      Y=(double)H[1];
      Y=Y/(double)K[1];
      I=1;
      printf(" I=%d, H[I]=%d, K[I]=%d, Y=%e \n",I,H[I],K[I],Y);
      fprintf(Outfp," %d, %d, \n",H[I],K[I]);
      if(J==1) goto L500;
//
      H[2]=A[2]*H[1]+1;
      K[2]=A[2]*K[1];
      Y=(double)H[2];
      Y=Y/(double)K[2];
      I=2;
      printf(" I=%d, H[I]=%d, K[I]=%d, Y=%e \n",I,H[I],K[I],Y);
      fprintf(Outfp," %d, %d, \n",H[I],K[I]);
      if(J==2) goto L500;
//
      for (I=3; I<=J; I++) {
	 for (k=1; k<=A[I]; k++) {
	    H[I]=k*H[I-1]+H[I-2];
	    K[I]=k*K[I-1]+K[I-2];
	    Y=(double)H[I];
	    Y=Y/(double)K[I];
	    printf(" I=%d, H[I]=%d, K[I]=%d, Y=%e \n",I,H[I],K[I],Y);
	    fprintf(Outfp," %d, %d, \n",H[I],K[I]);
	    fprintf(Outfp," %d, %d, \n",H[I-1]*(k+1),K[I-1]*(k+1));
	    }
	 temph=(A[I]+1)*H[I-1]+H[I-2];
	 tempk=(A[I]+1)*K[I-1]+K[I-2];
	 Y=(double)temph;
	 Y=Y/(double)tempk;
//	 printf(" I=%d, H[I]=%d, K[I]=%d, Y=%e \n",I,temph,tempk,Y);
	 printf("\n");
	 }
L500:
fclose(Outfp);
      return;
      }