www.pudn.com > calc.rar > davison.c

#include "integer.h" 
#include  
#include  
#include  
#include "fun.h" 
#include "stack.h" 
/*#include "unistd.h"*/ 
 
 
MPI *GLOBALA; 
MPI *GLOBALB; 
MPI *GLOBALC; 
MPI *GLOBALD; 
USL count, flagl, flagr, globalr, globall; 
char buff[20]; 
FILE *outfile; 
 
USL POSITIVITY(MPI *A, MPI *B, MPI *C, MPI *D){ 
	 
	if(A->S >= 0 && B->S >= 0 && C->S >= 0 && D->S >= 0){ 
		return(1); 
	}else{ 
		return(0); 
	} 
} 
 
 
/* NPROD(l,m) finds the least n such that A_n and D=A_0...A_n are  
 * non-negative. 
 * The matrix of global variables D=[GLOBALA,GLOBALB,GLOBALC,GLOBALD] is 
 * returned along with n. See Proposition 4.1, J.L. Davison, 'An algorithm for  
 * the continued fraction of e^{l/m}', Proceedings of the Eighth Manitoba  
 * Conference on Numerical Mathematics and Computing (Univ. Manitoba, Winnipeg, 
 * 1978), 169--179, Congress. Numer., XXII, Utilitas Math.  
 */ 
USL NPROD(USL l,USL m){ 
	USL k, s, temp; 
	MPI *TEMPM, *TEMPL, *T, *TEMP1, *TEMP2, *A1, *B1, *C1, *D1; 
	MPI *TEMP, *T1, *T2; 
	 
	TEMPM = CHANGE(m); 
	TEMPL = CHANGE(l); 
	temp = m + l; 
	GLOBALA =CHANGE(temp); 
	GLOBALB = CHANGE(m); 
	GLOBALC = CHANGE(m); 
	GLOBALD = SUBI(TEMPM, TEMPL); 
 
	for(k = 1; 1; k++){ 
		s = POSITIVITY(GLOBALA, GLOBALB, GLOBALC, GLOBALD); 
		if(s) 
			break; 
		T  = CHANGE((2 * k + 1) * m); 
		T1 = ADDI(GLOBALA, GLOBALB); 
		T2 = ADDI(GLOBALC, GLOBALD); 
		TEMP1 = MULTI(T, T1); 
		TEMP2 = MULTI(T, T2); 
		FREEMPI(T); 
		FREEMPI(T1); 
		FREEMPI(T2); 
		TEMP = MULTI(GLOBALA, TEMPL); 
		A1 = ADDI(TEMP1, TEMP); 
		FREEMPI(TEMP); 
		TEMP = MULTI(GLOBALB, TEMPL); 
		B1 = SUBI(TEMP1, TEMP); 
		FREEMPI(TEMP); 
		TEMP = MULTI(GLOBALC, TEMPL); 
		C1 = ADDI(TEMP2, TEMP); 
		FREEMPI(TEMP); 
		TEMP = MULTI(GLOBALD, TEMPL); 
		D1 = SUBI(TEMP2, TEMP); 
		FREEMPI(TEMP); 
		FREEMPI(TEMP1); 
		FREEMPI(TEMP2); 
		TEMP = GLOBALA; 
		GLOBALA = A1; 
		FREEMPI(TEMP); 
		TEMP = GLOBALB; 
		GLOBALB = B1; 
		FREEMPI(TEMP); 
		TEMP = GLOBALC; 
		GLOBALC = C1; 
		FREEMPI(TEMP); 
		TEMP = GLOBALD; 
		GLOBALD = D1; 
		FREEMPI(TEMP); 
	} 
	FREEMPI(TEMPL); 
	FREEMPI(TEMPM); 
	return(k - 1); 
} 
 
MPI *REFINE(MPI *A, MPI *B, MPI *C, MPI *D){ 
	MPI *G, *TEMP; 
 
	G = GCD(A, B); 
	TEMP = G; 
	G = GCD(G, C); 
	FREEMPI(TEMP); 
	TEMP = G; 
	G = GCD(G, D); 
	FREEMPI(TEMP); 
	return(G); 
} 
 
 
/* Raney factorization - G.N. Raney, Math, Annalen 206 (1973) 265-283. 
 * Input: a non-singular matrix A=[p,q;r,s], p,q,r,s>=0, A!=I_2, A!=[0,1;1,0]. 
 * With L=[1,0;1,1] and R=[1,1;0,1], we express A uniquely as 
 * a product of non-negative powers of L and R, (at least one is positive) 
 * followed by a row-balanced B. 
 * B=[a,b;c,d] is row-balanced if (ad) or (cb) and a,b,c>=0.  
 * The number k of powers of L and R is returned. 
 */ 
 
USL RANEY(MPI *P, MPI *Q, MPI *R, MPI *S){ 
	USL k, i, j; 
	MPI *PP, *QQ, *RR, *SS, *TEMP; 
 
	PP = COPYI(P); 
	QQ = COPYI(Q); 
	RR = COPYI(R); 
	SS = COPYI(S); 
	k = 0; 
	while(1){ 
                i = 0; 
		while(COMPAREI(PP, RR)>= 0 && COMPAREI(QQ, SS) >= 0){ 
			TEMP = PP; 
		        PP = SUBI(PP, RR); 
			FREEMPI(TEMP); 
			TEMP = QQ; 
		        QQ = SUBI(QQ, SS); 
			FREEMPI(TEMP); 
	                i++; 
		} 
		if(i > 0){ 
			globalr = globalr + i; 
			flagr = 1; 
			if(flagr * flagl){ 
				flagl = 0; 
				printf("a[%lu]:%lu\n", count, globall); 
				fprintf(outfile, "a[%lu]:%lu\n", count, globall); 
				globall = 0; 
				count++; 
			} 
			k++; 
		} 
                j = 0; 
		while(COMPAREI(RR, PP) >= 0 && COMPAREI(SS, QQ) >= 0){ 
			TEMP = RR; 
		        RR = SUBI(RR, PP); 
			FREEMPI(TEMP); 
			TEMP = SS; 
		        SS = SUBI(SS, QQ); 
			FREEMPI(TEMP); 
			j++; 
		} 
		if(j > 0){ 
			globall = globall + j; 
			flagl = 1; 
			if(flagr * flagl){ 
				flagr = 0; 
				printf("a[%lu]:%lu\n", count, globalr); 
				fprintf(outfile, "a[%lu]:%lu\n", count, globalr); 
				globalr = 0; 
				count++; 
			} 
			k++; 
		} 
		if((COMPAREI(PP, RR) < 0 && COMPAREI(QQ, SS) > 0) || (COMPAREI(PP, RR) > 0 && COMPAREI(QQ, SS) < 0)){ 
			break; 
		} 
	} 
	TEMP = GLOBALA; 
	GLOBALA = PP; 
	FREEMPI(TEMP); 
	TEMP = GLOBALB; 
	GLOBALB = QQ; 
	FREEMPI(TEMP); 
	TEMP = GLOBALC; 
	GLOBALC = RR; 
	FREEMPI(TEMP); 
	TEMP = GLOBALD; 
	GLOBALD = SS; 
	FREEMPI(TEMP); 
	/*fclose(outfile);*/ 
	return(k); 
} 
 
/* We perform the algorithm of J.L. Davison's paper. 
 * With n > =0, we first find the n* of Davison's Proposition 4.1 
 * and apply Raney's factorisation to A_0...A_k, for n*<=k<=n*+n. 
 * The number (count) of partial quotients of e^{l/m} found is returned. 
 * count becomes positive for all large n. 
 * We exit the program if count reaches 100000. 
 */ 
USL DAVISON(USL l, USL m, USL n){ 
	MPI *G, *ONE, *TEMP, *T, *T1, *T2, *TEMP1, *TEMP2; 
	MPI *A1, *B1, *C1, *D1, *TEMPL; 
	USL i, j, k, t; 
 
	printf("l,m,n=%lu,%lu,%lu\n",l,m,n); 
	strcpy(buff, "davison.out"); 
	outfile = fopen(buff, "w"); 
	count = 0; 
	flagr= 0; 
	flagl= 0; 
	globalr = 0; 
	globall= 0; 
	k = NPROD(l, m); 
	G = REFINE(GLOBALA, GLOBALB, GLOBALC, GLOBALD); 
	ONE = ONEI(); 
	if(RSV(G, ONE)>0){ 
		TEMP = GLOBALA; 
		GLOBALA = INT(GLOBALA, G); 
		FREEMPI(TEMP); 
		TEMP = GLOBALB; 
		GLOBALB = INT(GLOBALB, G); 
		FREEMPI(TEMP); 
		TEMP = GLOBALC; 
		GLOBALC = INT(GLOBALC, G); 
		FREEMPI(TEMP); 
		TEMP = GLOBALD; 
		GLOBALD = INT(GLOBALD, G); 
		FREEMPI(TEMP); 
	} 
	FREEMPI(G); 
	i = k; 
	j =k + n; 
	TEMPL = CHANGE(l); 
	while(i <= j){ 
		if(i > k){ 
			T  = CHANGE((2 * i + 1) * m); 
			T1 = ADDI(GLOBALA, GLOBALB); 
			T2 = ADDI(GLOBALC, GLOBALD); 
			TEMP1 = MULTI(T, T1); 
			TEMP2 = MULTI(T, T2); 
			FREEMPI(T); 
			FREEMPI(T1); 
			FREEMPI(T2); 
			TEMP = MULTI(GLOBALA, TEMPL); 
			A1 = ADDI(TEMP1, TEMP); 
			FREEMPI(TEMP); 
			TEMP = MULTI(GLOBALB, TEMPL); 
			B1 = SUBI(TEMP1, TEMP); 
			FREEMPI(TEMP); 
			TEMP = MULTI(GLOBALC, TEMPL); 
			C1 = ADDI(TEMP2, TEMP); 
			FREEMPI(TEMP); 
			TEMP = MULTI(GLOBALD, TEMPL); 
			D1 = SUBI(TEMP2, TEMP); 
			FREEMPI(TEMP); 
			FREEMPI(TEMP1); 
			FREEMPI(TEMP2); 
			TEMP = GLOBALA; 
			GLOBALA = A1; 
			FREEMPI(TEMP); 
			TEMP = GLOBALB; 
			GLOBALB = B1; 
			FREEMPI(TEMP); 
			TEMP = GLOBALC; 
			GLOBALC = C1; 
			FREEMPI(TEMP); 
			TEMP = GLOBALD; 
			GLOBALD = D1; 
			FREEMPI(TEMP); 
		} 
		i++; 
		t = RANEY(GLOBALA, GLOBALB, GLOBALC, GLOBALD); 
		if(count == 1000000){ 
			printf("1000000 partial quotients found\n"); 
			i = j + 1; 
		} 
                /* the input matrix will not be I_2 or [0,1;1,0] here. */ 
		G = REFINE(GLOBALA, GLOBALB, GLOBALC, GLOBALD); 
		if(RSV(G, ONE)>0){ 
			TEMP = GLOBALA; 
			GLOBALA = INT(GLOBALA, G); 
			FREEMPI(TEMP); 
			TEMP = GLOBALB; 
			GLOBALB = INT(GLOBALB, G); 
			FREEMPI(TEMP); 
			TEMP = GLOBALC; 
			GLOBALC = INT(GLOBALC, G); 
			FREEMPI(TEMP); 
			TEMP = GLOBALD; 
			GLOBALD = INT(GLOBALD, G); 
			FREEMPI(TEMP); 
		} 
		FREEMPI(G); 
	} 
	FREEMPI(TEMPL); 
	FREEMPI(ONE); 
	FREEMPI(GLOBALA); 
	FREEMPI(GLOBALB); 
	FREEMPI(GLOBALC); 
	FREEMPI(GLOBALD); 
	fflush(stdout); 
	printf("n* = %lu, n = %lu\n", k, n);  
	fprintf(outfile, "n* = %lu, n = %lu\n", k, n);  
	printf("The number of partial quotients found for e^(%lu/%lu) is %lu\n", l, m, count); 
	fprintf(outfile, "The number of partial quotients found for e^(%lu/%lu) is %lu\n", l, m, count); 
	fclose(outfile); 
	return(count); 
} 
 
MPI *DAVISONX(MPI *L, MPI *M, MPI *N){ 
	USL l, m, n, t; 
	MPI *G, *TEMP, *TEMPL, *TEMPM, *ONE; 
 
	if(L->S <= 0){ 
		printf("l<=0\n"); 
		return(NULL); 
	} 
	if(M->S <= 0){ 
		printf("m<=0\n"); 
		return(NULL); 
	} 
	if(N->S < 0){ 
		printf("n<0\n"); 
		return(NULL); 
	} 
	TEMP = CHANGEI(1000); 
	if(RSV(L, TEMP) > 0){ 
		printf("l>1000\n"); 
		FREEMPI(TEMP); 
		return(NULL); 
	} 
	if(RSV(M, TEMP) > 0){ 
		printf("m>1000\n"); 
		FREEMPI(TEMP); 
		return(NULL); 
	} 
	FREEMPI(TEMP); 
	TEMP = CHANGEI(100000); 
	if(RSV(N, TEMP) > 0){ 
		printf("n>=10^6\n"); 
		FREEMPI(TEMP); 
		return(NULL); 
	} 
	FREEMPI(TEMP); 
	TEMPL = COPYI(L); 
	TEMPM = COPYI(M); 
	ONE = ONEI(); 
	G = GCD(L, M); 
	if(RSV(G, ONE)>0){ 
		TEMP = TEMPL; 
		TEMPL = INT0(L, G); 
		FREEMPI(TEMPL); 
		TEMP = TEMPM; 
		TEMPM = INT0(M, G); 
		FREEMPI(TEMP); 
	} 
	FREEMPI(G); 
	FREEMPI(ONE); 
	l = CONVERTI(TEMPL); 
	m = CONVERTI(TEMPM); 
	FREEMPI(TEMPL); 
	FREEMPI(TEMPM); 
	n = CONVERTI(N); 
	t = DAVISON(l, m, n); 
	return (CHANGE(t)); 
} 
 
/* Raney factorization - G.N. Raney, Math, Annalen 206 (1973) 265-283. 
 * Input: a non-singular matrix A=[p,q;r,s], p,q,r,s>=0, A!=I_2, A!=[0,1;1,0]. 
 * With L=[1,0;1,1] and R=[1,1;0,1], we express A uniquely as 
 * a product of non-negative powers of L and R, (at least one is positive) 
 * followed by a row-balanced B. 
 * B=[a,b;c,d] is row-balanced if (ad) or (cb) and a,b,c>=0.  
 * The number k of powers of L and R is returned. 
 */ 
 
USL RANEY1(MPI *P, MPI *Q, MPI *R, MPI *S){ 
	USL k, i, j; 
	MPI *PP, *QQ, *RR, *SS, *TEMP; 
	char buff1[20]; 
	FILE *outfile1; 
 
	strcpy(buff1, "raney.out"); 
	outfile1 = fopen(buff1, "w"); 
	printf("["); 
	fprintf(outfile1, "["); 
	PRINTI(P); 
	FPRINTI(outfile1, P); 
	printf(","); 
	fprintf(outfile1, ","); 
	PRINTI(Q); 
	FPRINTI(outfile1, Q); 
	printf(";"); 
	fprintf(outfile1, ";"); 
	PRINTI(R); 
	FPRINTI(outfile1, R); 
	printf(","); 
	fprintf(outfile1, ","); 
	PRINTI(S); 
	FPRINTI(outfile1, S); 
	printf("]="); 
	fprintf(outfile1, "]="); 
	PP = COPYI(P); 
	QQ = COPYI(Q); 
	RR = COPYI(R); 
	SS = COPYI(S); 
	k = 0; 
	while(1){ 
                i = 0; 
		while(COMPAREI(PP, RR)>= 0 && COMPAREI(QQ, SS) >= 0){ 
			TEMP = PP; 
		        PP = SUBI(PP, RR); 
			FREEMPI(TEMP); 
			TEMP = QQ; 
		        QQ = SUBI(QQ, SS); 
			FREEMPI(TEMP); 
	                i++; 
		} 
		if(i > 0){ 
			k++; 
			if(i > 1){ 
				printf("R^%lu", i);	 
				fprintf(outfile1, "R^%lu", i);	 
			}else{ 
				printf("R");	 
				fprintf(outfile1, "R");	 
			} 
		} 
                j = 0; 
		while(COMPAREI(RR, PP) >= 0 && COMPAREI(SS, QQ) >= 0){ 
			TEMP = RR; 
		        RR = SUBI(RR, PP); 
			FREEMPI(TEMP); 
			TEMP = SS; 
		        SS = SUBI(SS, QQ); 
			FREEMPI(TEMP); 
			j++; 
		} 
		if(j > 0){ 
			if(j > 1){ 
				printf("L^%lu", j);	 
				fprintf(outfile1, "L^%lu", j);	 
			}else{ 
				printf("L");	 
				fprintf(outfile1, "L");	 
			} 
			k++; 
 
		} 
		if((COMPAREI(PP, RR) < 0 && COMPAREI(QQ, SS) > 0) || (COMPAREI(PP, RR) > 0 && COMPAREI(QQ, SS) < 0)){ 
			break; 
		} 
	} 
	printf("D, where D=["); 
	fprintf(outfile1, "D, where D=["); 
	PRINTI(PP); 
	FPRINTI(outfile1, PP); 
	printf(","); 
	fprintf(outfile1, ","); 
	PRINTI(QQ); 
	FPRINTI(outfile1, QQ); 
	printf(";"); 
	fprintf(outfile1, ";"); 
	PRINTI(RR); 
	FPRINTI(outfile1, RR); 
	printf(","); 
	fprintf(outfile1, ","); 
	PRINTI(SS); 
	FPRINTI(outfile1, SS); 
	printf("] is row-balanced\n"); 
	fprintf(outfile1, "] is row-balanced\n"); 
	FREEMPI(PP); 
	FREEMPI(QQ); 
	FREEMPI(RR); 
	FREEMPI(SS); 
	fclose(outfile1); 
	return(k); 
} 
 
MPI *RANEY1X(MPI *P, MPI *Q, MPI *R, MPI *S){ 
	MPI *TEMP1, *TEMP2, *TEMP; 
	USI t; 
	USL k; 
 
		if(P->S < 0){ 
			printf("P < 0\n"); 
			return(NULL); 
		} 
		if(Q->S < 0){ 
			printf("Q < 0\n"); 
			return(NULL); 
		} 
		if(R->S < 0){ 
			printf("R < 0\n"); 
			return(NULL); 
		} 
		if(S->S < 0){ 
			printf("S < 0\n"); 
			return(NULL); 
		} 
		if(EQONEI(P) && EQZEROI(Q) && EQZEROI(R) && EQONEI(S)){ 
			printf("A = I_2\n"); 
			return(NULL); 
		} 
		if(EQZEROI(P) && EQONEI(Q) && EQONEI(R) && EQZEROI(S)){ 
			printf("A = J = [0,1;1,0]\n"); 
			return(NULL); 
		} 
		TEMP1 = MULTI(P, S); 
		TEMP = CHANGEI(1000); 
		TEMP2 = MULTI(Q, R); 
		t = EQUALI(TEMP1, TEMP2); 
		FREEMPI(TEMP1); 
		FREEMPI(TEMP2); 
		if(t){ 
			printf("P*S-Q*R=0\n"); 
			return(NULL); 
		} 
		k = RANEY1(P, Q, R, S); 
		return(CHANGE(k)); 
} 
 
/* 
 * This algorithm expresses a unimodular matrix A !=I_2 or U=[0 1] 
 *                                                           [1 0] 
 * with non-negative coefficients, as a product of one of the 
 * following forms: 
 * P, UP, PU, or UPU, where P is a product of matrices of the form 
 * [a 1], a>0. 
 * [1 0] 
 * The representation is unique. See Kjell Kolden, 'Continued fractions 
 * and linear substitutions', Arch. Math. Naturvid. 50 (1949), 141-196. 
 * The number n of matrices in the product [d[0]1]***[d[n-1]1] 
 * is returned.                            [  1 0]   [  1   0] 
 */ 
 
USL UNIMODULAR(MPI *P, MPI *Q, MPI *R, MPI *S){ 
	USL i; 
	MPI *D, *D1, *D2, *TEMP, *TEMPP, *TEMPQ, *PP, *QQ, *RR, *SS; 
	char buff1[20]; 
	FILE *outfile1; 
 
	i = 0;	 
	strcpy(buff1, "unimodular.out"); 
	outfile1 = fopen(buff1, "w"); 
	printf("["); 
	fprintf(outfile1, "["); 
	PRINTI(P); 
	FPRINTI(outfile1, P); 
	printf(","); 
	fprintf(outfile1, ","); 
	PRINTI(Q); 
	FPRINTI(outfile1, Q); 
	printf(";"); 
	fprintf(outfile1, ";"); 
	PRINTI(R); 
	FPRINTI(outfile1, R); 
	printf(","); 
	fprintf(outfile1, ","); 
	PRINTI(S); 
	FPRINTI(outfile1, S); 
	printf("]="); 
	fprintf(outfile1, "]="); 
	PP = COPYI(P); 
	QQ = COPYI(Q); 
	RR = COPYI(R); 
	SS = COPYI(S); 
	while(1){ 
		if(SS->S == 0){ 
			printf("U["); 
			fprintf(outfile1, "U["); 
			PRINTI(PP); 
			FPRINTI(outfile1, PP); 
			printf("]"); 
			fprintf(outfile1, "]"); 
		        fflush(stdout); 
			i++; 
			break; 
		} 
		if(RR->S == 0){ 
			printf("U["); 
			fprintf(outfile1, "U["); 
			PRINTI(QQ); 
			FPRINTI(outfile1, QQ); 
			printf("]"); 
			fprintf(outfile1, "]"); 
		        fflush(stdout); 
			i++; 
			printf("U[0]"); 
			fprintf(outfile1, "U[0]"); 
		        fflush(stdout); 
			i++; 
			break; 
		} 
		D1 = INT0(QQ, SS); 
		D2 = INT0(PP, RR); 
		D = MINMPI(D1, D2); 
		FREEMPI(D1); 
		FREEMPI(D2); 
		TEMPP = PP; 
		PP = RR; 
		TEMP = MULTI(D, RR); 
		RR = SUBI(TEMPP, TEMP); 
		FREEMPI(TEMPP); 
		FREEMPI(TEMP); 
		TEMPQ = QQ; 
		QQ = SS; 
		TEMP = MULTI(D, SS); 
		SS = SUBI(TEMPQ, TEMP); 
		FREEMPI(TEMPQ); 
		FREEMPI(TEMP); 
		printf("U["); 
		fprintf(outfile1, "U["); 
		PRINTI(D); 
		FPRINTI(outfile1, D); 
		printf("]"); 
		fprintf(outfile1, "]"); 
		FREEMPI(D); 
		fflush(stdout); 
		i++; 
	} 
	printf("\n"); 
	fprintf(outfile1, "\n"); 
	FREEMPI(PP); 
	FREEMPI(QQ); 
	FREEMPI(RR); 
	FREEMPI(SS); 
	fclose(outfile1); 
	return(i); 
} 
 
MPI *UNIMODULARX(MPI *P, MPI *Q, MPI *R, MPI *S){ 
	USI t1, t2; 
	USL i; 
	MPI *DET, *TEMP1, *TEMP2; 
 
	if(P->S < 0){ 
		printf("P < 0\n"); 
		return(NULL); 
	} 
	if(Q->S < 0){ 
		printf("Q < 0\n"); 
		return(NULL); 
	} 
	if(R->S < 0){ 
		printf("R < 0\n"); 
		return(NULL); 
	} 
	if(S->S < 0){ 
		printf("S < 0\n"); 
		return(NULL); 
	} 
	TEMP1 = MULTI(P, S); 
	TEMP2 = MULTI(Q, R); 
	DET = SUBI(TEMP1, TEMP2); 
	FREEMPI(TEMP1); 
	FREEMPI(TEMP2); 
	t1 = EQONEI(DET); 
	t2 = EQMINUSONEI(DET); 
	FREEMPI(DET); 
	if(EQONEI(P) && EQZEROI(Q) &&EQZEROI(R) &&EQONEI(S)){ 
		printf("A=I_2\n"); 
		return(NULL); 
	} 
	if(EQZEROI(P) && EQONEI(Q) &&EQONEI(R) &&EQZEROI(S)){ 
		printf("A = U = [0,1;1,0]\n"); 
		return(NULL); 
	} 
	if(t1 == 0 && t2 == 0){ 
		printf("matrix is not unimodular\n"); 
		return(NULL); 
	} 
	i = UNIMODULAR(P, Q, R, S); 
	return(CHANGE(i)); 
}