﻿ /* COMPUTE GENERALIZED CONTINUED
```/* 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;
}
```