log.c

calc大数库

#include "integer.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "fun.h"
#include "stack.h"
#ifdef _WIN32
#include "unistd_DOS.h"
#else
#include <unistd.h>
#endif

void INTLOG(Stack s)
/* This performs variant 0 of Shank's log_b(a),
 * working as in bc scale 2r, but doing it in integers.
 * See http://www.maths.uq.edu.au/~krm/log.pdf
 * We print out the partial quotients n[0],...,n[r-1] if a>b,
 * n[0],...,n[r] if a<b,
 * and, with c=10^{2r},  int[a[i]*c], where the a[i] are Shank's:
 * b[0]=[a[0]*c],,...,b[r+1]=[a[r+1]*c] if a[0]=a < a[1]=b,
 * b[0]=[a[1]*c],,...,b[r]=[a[r+1]*c] if a[0]=a > a[1]=b.
 * The output is not guaranteed to be correct, but seems certain to be so.
 * See D. Shanks, "A logarithm algorithm", Math. Tables and Other Aids to
 * Computation 8, (1954). 60--64.
 */
{
	USI j;
	USL i, r, rr, n;
	MPI *AA, *BB, *C, *TMP1, *TMP2;

	MPI *A = stackPop(s);
	MPI *B = stackPop(s);
	MPI *R = stackPop(s);
	MPI *N = stackPop(s);
	i = 0;
	j = 2;
	rr = CONVERTI(R);
	FREEMPI(R);
	n = CONVERTI(N);
	FREEMPI(N);
	if (RSV(A, B) == -1)
		rr++;
	r = rr + rr;
	C = POWER_I(10, r);
	AA = MULTI(A, C);
	if (n){
		printf("A[0] = "); 
		PRINTI(AA); 
		printf("\n");
	}
	BB = MULTI(B, C);
	if (n){
		printf("A[1] = "); 
		PRINTI(BB); 
		printf("\n");
	}
	while (rr){
		if (RSV(AA, BB) >= 0){
			TMP1 = MULTI(AA, C);
			TMP2 = AA;
			AA = INT(TMP1, BB);
			FREEMPI(TMP1);
			FREEMPI(TMP2);
			i++;
		}
		else{
			if (n){
				printf("n[%u]=%lu, A[%u]=", j - 2, i, j);
				PRINTI(AA);
			}
			else
				printf("n[%u]=%lu", j - 2, i);
			printf("\n");
			j++;
			TMP1 = BB;
			BB = AA;
			AA = TMP1;
			rr--;
			i = 0; /* reset the partial quotient count */

		}
	}
	FREEMPI(AA);
	FREEMPI(BB);
	FREEMPI(A);
	FREEMPI(B);
	FREEMPI(C);
	return;
}

void LOG_(Stack s)
/* This performs Shank's log a to base b algorithm, working in effect with 
 * scale 2r, but doing it in integers.
 * We work with scale = 2r and believe that this is sufficient to correctly
 * output partial quotients a[0],...,a[s], where s=r-1 if n>b, but r if n<b.
 * the decimal expansion of log_b{n} is printed, truncated correct to
 * as many decimal places as possible: the r-1th convergent is compared with
 * rth convergent and decimal expansion is truncated where they differ.
 * e=0 prints only the partial quotients, convergents and decimal expansion,
 * while e !=0 is verbose.
 * Initially written in by Alan Offer and improved by Sean Seefried, Dec 1999.
 * Modified by Sean Seefried on 17th January 2000.
 * Calc version by Keith Matthews, 13th March 2000.
 * See D. Shanks, "A logarithm algorithm", Math. Tables and Other Aids to
 * Computation 8, (1954). 60--64.
 */
{
	USI j;
	int tt;
	USL i, r, rr, n;
	MPI *AA, *BB, *C, *TMP1, *TMP2;
	MPI *PN, *QP, *QN, *PP, *X, *Y, *Z;
	char buff[20];
        FILE *outfile;

	MPI *A = stackPop(s);
	MPI *B = stackPop(s);
	MPI *R = stackPop(s);
	MPI *N = stackPop(s);
	i = 0;
	j = 2;
	PN = ONEI();
	QP = ONEI();
	QN = ZEROI();
	PP = ZEROI();
	rr = CONVERTI(R);
	FREEMPI(R);
	n = CONVERTI(N);
	FREEMPI(N);
	if (RSV(A, B) == -1)
		rr++;
	r = rr + rr;
	C = POWER_I(10, r);
	AA = MULTI(A, C);
	printf("Partial Quotient, Convergent\n");
	printf("----------------------------\n");
	if (n){
		printf("A[0] = "); 
		PRINTI(AA); 
		printf("\n");
	}
	strcpy(buff, "log.out");
	if(access(buff, R_OK) == 0)
		  unlink(buff);
	outfile = fopen(buff, "w");
	fprintf(outfile, "A[0] = "); 
	FPRINTI(outfile, AA); 
	fprintf(outfile, "\n");
	BB = MULTI(B, C);
	if (n){
		printf("A[1] = "); 
		PRINTI(BB); 
		printf("\n");
	}
	fprintf(outfile, "A[1] = "); 
	FPRINTI(outfile, BB); 
	fprintf(outfile, "\n");
	while (rr){
		if (RSV(AA, BB) >= 0){
			TMP1 = MULTI(AA, C);
			TMP2 = AA;
			AA = INT(TMP1, BB);
			FREEMPI(TMP1);
			FREEMPI(TMP2);
			i++;
			TMP1 = PN;
			PN = ADDI(PN, PP);
			FREEMPI(TMP1);
			TMP1 = QN;
			QN = ADDI(QN, QP);
			FREEMPI(TMP1);
			if(EQUALI(BB, C)){
				PRINTI(A);
				FPRINTI(outfile, A);
				printf("^");
				fprintf(outfile, "^");
				PRINTI(PP);
				FPRINTI(outfile, PP);
				printf("=");
				fprintf(outfile, "=");
				PRINTI(B);
				FPRINTI(outfile, B);
				printf("^");
				fprintf(outfile, "^");
				PRINTI(QP);
				FPRINTI(outfile, QP);
				printf("\n");
				fprintf(outfile, "\n");
				FREEMPI(AA);
				FREEMPI(BB);
				FREEMPI(A);
				FREEMPI(B);
				FREEMPI(C);
				FREEMPI(PN);
				FREEMPI(PP);
				FREEMPI(QP);
				FREEMPI(QN);
				return;
			}	
		}
		else{
			if (n){
				printf("n[%u]=%lu, A[%u]=", j - 2, i, j);
				PRINTI(AA);
			}
			else
				printf("n[%u]=%lu", j - 2, i);
			fprintf(outfile, "n[%u]=%lu, A[%u]=", j - 2, i, j);
			FPRINTI(outfile, AA);
			j++;
			TMP1 = BB;
			BB = AA;
			AA = TMP1;
			TMP1 = PN;
			PN = PP;
			PP = TMP1;
			TMP1 = QN;
			QN = QP;
			QP = TMP1;
			printf(", p[%u]/q[%u]=", j - 3, j - 3); 
			PRINTI(QP); 
			printf("/");
			PRINTI(PP);
			fprintf(outfile, ", p[%u]/q[%u]=", j - 3, j - 3); 
			FPRINTI(outfile, QP); 
			fprintf(outfile, "/");
			FPRINTI(outfile, PP);
			printf("\n");
			fprintf(outfile, "\n");
			rr--;
			i = 0; /* reset the partial quotient count */
			
		}
	}
	printf("The log of "); 
	fprintf(outfile, "The log of "); 
	PRINTI(A); 
	FPRINTI(outfile, A); 
	printf(" to base "); 
	fprintf(outfile, " to base "); 
	PRINTI(B); 
	FPRINTI(outfile, B); 
	if (j == 3){
		printf(" has integer part ");
		fprintf(outfile, " has integer part ");
		PRINTI(QP);
		FPRINTI(outfile, QP);
		printf("\n");
		fprintf(outfile, "\n");
		FREEMPI(AA);
		FREEMPI(BB);
		FREEMPI(A);
		FREEMPI(B);
		FREEMPI(C);
		FREEMPI(PN);
		FREEMPI(PP);
		FREEMPI(QP);
		FREEMPI(QN);
		fclose(outfile);
		return;
	}
	printf(" equals ");
	fprintf(outfile, " equals ");
	if (PN->S){
		fclose(outfile);
		tt = COMPARE_DIGITS(QP, PP, QN, PN, 10);
		outfile = fopen(buff, "a");
		if (tt == -1){
			X = MULTI(QP, PN);
			Y = MULTI(PP, QN);
			printf(" has integer part ");
			fprintf(outfile, " has integer part ");
			if (RSV(X, Y)< 0)
				Z = INTI(QP, PP);
			else
				Z = INTI(QN, PN);
			PRINTI(Z);
			FPRINTI(outfile, Z);
			FREEMPI(Z);
			FREEMPI(X);
			FREEMPI(Y);
			printf("\n");
			fprintf(outfile, "\n");
			FREEMPI(AA);
			FREEMPI(BB);
			FREEMPI(A);
			FREEMPI(B);
			FREEMPI(C);
			FREEMPI(PN);
			FREEMPI(PP);
			FREEMPI(QP);
			FREEMPI(QN);
			fclose(outfile);
			return;
		}
		else{
			printf(" truncated to %u decimal places\n", tt);
			fprintf(outfile, " truncated to %u decimal places\n", tt);
		}
	}
	FREEMPI(AA);
	FREEMPI(BB);
	FREEMPI(A);
	FREEMPI(B);
	FREEMPI(C);
	FREEMPI(PN);
	FREEMPI(PP);
	FREEMPI(QP);
	FREEMPI(QN);
	fclose(outfile);
	return;
}

int COMPARE_DIGITS(MPI *M, MPI *N, MPI *U, MPI *V, USL b)
/* m,n,u,v are positive integers, b is the base.
 * We compare the base b expansions. If the integer parts
 * of m/n and u/v are equal, we go on to compare the expansions
 * after the decimal point. We print out the t common decimals and
 * the number t.
 */
{
	USI i, o;
	MPI *R, *S, *T, *P, *RR, *PP;
	MPI *MM, *UU, *TMP1;
	char buff[20];
        FILE *outfile;

	R = INTI(M, N);
	RR = INTI(U, V);
	strcpy(buff, "log.out");
	outfile = fopen(buff, "a");
	if (EQUALI(R, RR)){
		/* integer parts are equal */
		PRINTI(R);
		FPRINTI(outfile, R);
	}
	else{
		FREEMPI(R);
		FREEMPI(RR);
		fclose(outfile);
		return(-1);
	}
	TMP1 = MULTI(R, N);
	MM = SUBI(M, TMP1);
	FREEMPI(TMP1);
	TMP1 = MULTI(RR, V);
	UU = SUBI(U, TMP1);
	FREEMPI(TMP1);
	i = 0;
	P = MULT_I(MM, b);
	FREEMPI(MM);
	PP = MULT_I(UU, b);
	FREEMPI(UU);
	S = INT(P, N); /* the first decimal digit of {m/n} */
	T = INT(PP, V);/* the first decimal digit of {u/v} */
	
	FREEMPI(R);
	FREEMPI(RR);
	R = MOD(P, N);
	RR = MOD(PP, V);
	FREEMPI(P);
	FREEMPI(PP);
	printf(".");
	fprintf(outfile, ".");
	o = EQUALI(S, T);
	while(o){
		PRINTI(S);
		FPRINTI(outfile, S);
		FREEMPI(S);
		FREEMPI(T);
		P = MULT_I(R, b);
		PP = MULT_I(RR, b);
		FREEMPI(R);
		FREEMPI(RR);
		R = MOD(P, N);
		RR = MOD(PP, V);
		S = INT(P, N); /* the first decimal digit of {m/n} */
		T = INT(PP, V);/* the first decimal digit of {u/v} */
		o = EQUALI(S, T);
		FREEMPI(P);
		FREEMPI(PP);
		i++;
	}
	FREEMPI(S);
	FREEMPI(T);
	FREEMPI(R);
	FREEMPI(RR);
	printf("\n");
	fprintf(outfile, "\n");
	fclose(outfile);
	return (i); 
	/* the no. of common decimal digits after the decimal point */
}

void INTLOG1(Stack s)
/* This performs variant 1 of Shank's log_b(a),
 * See http://www.maths.uq.edu.au/~krm/log.pdf
 */
{
	USI j;
	USL i, r, rr, n;
	MPI *AA, *BB, *C, *TMP1, *TMP2, *E, *F;

	MPI *A = stackPop(s);
	MPI *B = stackPop(s);
	MPI *R = stackPop(s);
	MPI *N = stackPop(s);
	i = 0;
	j = 2;
	rr = CONVERTI(R);
	FREEMPI(R);
	n = CONVERTI(N);
	FREEMPI(N);
	if (RSV(A, B) == -1)
		rr++;
	r = rr + rr;
	C = POWER_I(10, r);
	F = POWER_I(10, rr + 2);
	E = ADD0I(C, F);
	FREEMPI(F);
	AA = MULTI(A, C);
	if (n){
		printf("A[0] = "); PRINTI(AA); printf("\n");
	}
	BB = MULTI(B, C);
	if (n){
		printf("A[1] = "); PRINTI(BB); printf("\n");
	}
	while (RSV(BB, E)>= 0){
		if (RSV(AA, BB) >= 0){
			TMP1 = MULTI(AA, C);
			TMP2 = AA;
			AA = INT(TMP1, BB);
			FREEMPI(TMP1);
			FREEMPI(TMP2);
			i++;
		}
		else{
			if (n){
				printf("n[%u]=%lu, A[%u]=", j - 2, i, j);
				PRINTI(AA);
			}
			else
				printf("n[%u]=%lu", j - 2, i);
			printf("\n");
			j++;
			TMP1 = BB;
			BB = AA;
			AA = TMP1;
			/*rr--;*/
			i = 0; /* reset the partial quotient count */

		}
	}
	FREEMPI(AA);
	FREEMPI(BB);
	FREEMPI(A);
	FREEMPI(B);
	FREEMPI(C);
	FREEMPI(E);
	return;
}

void INTLOG2(Stack s)
/* This performs variant 2 of Shank's log_b(a),
 * See http://www.maths.uq.edu.au/~krm/log.pdf
 */
{
	USI j;
	USL d, i, r, rr, n;
	MPI *AA, *BB, *C, *TMP1, *TMP2, *E, *F;
	MPI *U, *X, *Y, *Z;

	MPI *A = stackPop(s);
	MPI *B = stackPop(s);
	MPI *R = stackPop(s);
	MPI *N = stackPop(s);
	i = 0;
	j = 2;
	rr = CONVERTI(R);
	FREEMPI(R);
	n = CONVERTI(N);
	FREEMPI(N);
	if (RSV(A, B) == -1)
		rr++;
	r = rr + rr;
	C = POWER_I(10, r);
	F = POWER_I(10, rr + 2);
	E = ADD0I(C, F);
	FREEMPI(F);
	AA = MULTI(A, C);
	if (n){
		printf("A[0] = "); PRINTI(AA); printf("\n");
	}
	BB = MULTI(B, C);
	if (n){
		printf("A[1] = "); PRINTI(BB); printf("\n");
	}
	while (RSV(BB, E)>= 0){
		Y = ONEI();
		X = COPYI(AA);
		while(1){
			U = INT0(X, Y);
			TMP1 = X;
			TMP2 = Y;
			X = MULTI(C, X);
			Y = MULTI(BB, Y);
			FREEMPI(TMP1);
			FREEMPI(TMP2);
			Z = MULTI(C, Y);
			d = RSV(Z, X);
			FREEMPI(Z);
			if (d == 1)
				break;
			i++;
			FREEMPI(U);
		}
		FREEMPI(X);
		FREEMPI(Y);
		FREEMPI(AA);
		AA = COPYI(U);
		FREEMPI(U);
		if (n){
			printf("n[%u]=%lu, A[%u]=", j - 2, i, j);
			PRINTI(AA);
		}
		else
			printf("n[%u]=%lu", j - 2, i);
		printf("\n");
		j++;
		TMP1 = BB;
		BB = AA;
		AA = TMP1;
		i = 0; /* reset the partial quotient count */
	}
	FREEMPI(AA);
	FREEMPI(BB);
	FREEMPI(A);
	FREEMPI(B);
	FREEMPI(C);
	FREEMPI(E);
	return;
}

void LOG1(Stack s)
/* This performs variant 1 of Shank's log a to base b algorithm, working in 
 * effect with scale 2r, but doing it in integers.
 */
{
	USI j;
	int tt;
	USL i, r, rr, n;
	MPI *AA, *BB, *C, *TMP1, *TMP2;
	MPI *PN, *QP, *QN, *PP, *X, *Y, *Z;
	MPI *E, *F;
	char buff[20];
        FILE *outfile;

	MPI *A = stackPop(s);
	MPI *B = stackPop(s);
	MPI *D = stackPop(s);
	MPI *R = stackPop(s);
	MPI *N = stackPop(s);
	i = 0;
	j = 2;
	PN = ONEI();
	QP = ONEI();
	QN = ZEROI();
	PP = ZEROI();
	rr = CONVERTI(R);
	FREEMPI(R);
	n = CONVERTI(N);
	FREEMPI(N);
	if (RSV(A, B) == -1)
		rr++;
	r = rr + rr;
	C = POWERI(D, r);
	F = POWERI(D, rr + 2);
	FREEMPI(D);
        E = ADD0I(C, F);
	FREEMPI(F);

	AA = MULTI(A, C);
	printf("Partial Quotient, Convergent\n");
	printf("----------------------------\n");
	if (n){
		printf("A[0] = "); 
		PRINTI(AA); 
		printf("\n");
	}
	strcpy(buff, "log1.out");
	if(access(buff, R_OK) == 0)
		  unlink(buff);

	outfile = fopen(buff, "w");
	fprintf(outfile, "A[0] = "); 
	FPRINTI(outfile, AA); 
	fprintf(outfile, "\n");
	BB = MULTI(B, C);
	if (n){
		printf("A[1] = "); 
		PRINTI(BB); 
		printf("\n");
	}
	fprintf(outfile, "A[1] = "); 
	FPRINTI(outfile, BB); 
	fprintf(outfile, "\n");
	while (RSV(BB, E) >= 0){
		if (RSV(AA, BB) >= 0){
			TMP1 = MULTI(AA, C);
			TMP2 = AA;
			AA = INT(TMP1, BB);
			FREEMPI(TMP1);
			FREEMPI(TMP2);
			i++;
			TMP1 = PN;
			PN = ADDI(PN, PP);
			FREEMPI(TMP1);
			TMP1 = QN;
			QN = ADDI(QN, QP);
			FREEMPI(TMP1);
			if(EQUALI(BB, C)){
				PRINTI(A);
				FPRINTI(outfile, A);
				printf("^");
				fprintf(outfile, "^");
				PRINTI(PP);
				FPRINTI(outfile, PP);
				printf("=");
				fprintf(outfile, "=");
				PRINTI(B);
				FPRINTI(outfile, B);
				printf("^");
				fprintf(outfile, "^");
				PRINTI(QP);
				FPRINTI(outfile, QP);
				printf("\n");
				fprintf(outfile, "\n");
				FREEMPI(AA);
				FREEMPI(BB);
				FREEMPI(A);
				FREEMPI(B);
				FREEMPI(C);
				FREEMPI(PN);
				FREEMPI(PP);
				FREEMPI(QP);
				FREEMPI(QN);
				return;
			}	
		}
		else{
			if (n){
				printf("n[%u]=%lu, A[%u]=", j - 2, i, j);
				PRINTI(AA);
			}
			else
				printf("n[%u]=%lu", j - 2, i);
			fprintf(outfile, "n[%u]=%lu, A[%u]=", j - 2, i, j);
			FPRINTI(outfile, AA);
			j++;
			TMP1 = BB;
			BB = AA;
			AA = TMP1;
			TMP1 = PN;
			PN = PP;
			PP = TMP1;
			TMP1 = QN;
			QN = QP;
			QP = TMP1;
			printf(", p[%u]/q[%u]=", j - 3, j - 3); 
			PRINTI(QP); 
			printf("/");
			PRINTI(PP);
			fprintf(outfile, ", p[%u]/q[%u]=", j - 3, j - 3); 
			FPRINTI(outfile, QP); 
			fprintf(outfile, "/");
			FPRINTI(outfile, PP);
			printf("\n");
			fprintf(outfile, "\n");
		/*	rr--; */
			i = 0; /* reset the partial quotient count */
			
		}
	}
	printf("The log of "); 
	fprintf(outfile, "The log of "); 
	PRINTI(A); 
	FPRINTI(outfile, A); 
	printf(" to base "); 
	fprintf(outfile, " to base "); 
	PRINTI(B); 
	FPRINTI(outfile, B); 
	if (j == 3){
		printf(" has integer part ");
		fprintf(outfile, " has integer part ");
		PRINTI(QP);
		FPRINTI(outfile, QP);
		printf("\n");
		fprintf(outfile, "\n");
		FREEMPI(AA);
		FREEMPI(BB);
		FREEMPI(A);
		FREEMPI(B);
		FREEMPI(C);
		FREEMPI(PN);
		FREEMPI(PP);
		FREEMPI(QP);
		FREEMPI(QN);
		fclose(outfile);
		return;
	}
	printf(" equals ");
	fprintf(outfile, " equals ");
	if (PN->S){
		fclose(outfile);
		tt = COMPARE_DIGITS(QP, PP, QN, PN, 10);
		outfile = fopen(buff, "a");
		if (tt == -1){