log.c

calc大数库

			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);
	FREEMPI(E);
	fclose(outfile);
	return;
}

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

	MPI *X = stackPop(S);
	MPI *Y = stackPop(S);
	MPI *D = stackPop(S);
	MPI *R = stackPop(S);
	MPI *RR = stackPop(S);
	i = 0;
	j = 0;
	
	s = 2;
	ss = 2;
	r = CONVERTI(R);
	rr = CONVERTI(RR);
	FREEMPI(R);
	FREEMPI(RR);
	if(r >= rr){
		FREEMPI(X);
		FREEMPI(Y);
		FREEMPI(D);
		execerror("arg 4 >= arg5", "");
	}
	C = POWERI(D, (USI)r);
	A = MULTI(X, C);
	B = MULTI(Y, C);
	CC = POWERI(D, (USI)rr);
	AA = MULTI(X, CC);
	BB = MULTI(Y, CC);
	FREEMPI(D);
	FREEMPI(X);
	FREEMPI(Y);
	while (RSV(B, C) >= 0 && RSV(BB, CC) >= 0){
		while (RSV(A, B) >= 0){
			TMP1 = MULTI(A, C);
			TMP2 = A;
			A = INT(TMP1, B);
			FREEMPI(TMP1);
			FREEMPI(TMP2);
			i++;
		}
		n = i;
printf("n=%lu\n", n);
		s++;
		TMP1 = B;
		B = A;
		A = TMP1;
		i = 0; /* reset the partial quotient count */

printf("AA=");
PRINTI(AA);
printf("\n");
printf("BB=");
PRINTI(BB);
printf("\n");
		while (RSV(AA, BB) >= 0){
			TMP1 = MULTI(AA, CC);
			TMP2 = AA;
			AA = INT(TMP1, BB);
			FREEMPI(TMP1);
			FREEMPI(TMP2);
			j++;
printf("AA=");
PRINTI(AA);
printf("\n");
printf("BB=");
PRINTI(BB);
printf("\n");
		}
		nn = j;
printf("nn=%lu\n", nn);
		ss++;
		TMP1 = BB;
		BB = AA;
		AA = TMP1;
		j = 0; /* reset the partial quotient count */
		if (n != nn){
			printf("%lu != %lu\n", n, nn);
			printf("breaking\n");
			break;
		}
		printf("n[%lu]=%lu\n", s - 3, n);
	}
	FREEMPI(A);
	FREEMPI(B);
	FREEMPI(AA);
	FREEMPI(BB);
	FREEMPI(C);
	FREEMPI(CC);
	return;
}

void LOGGG(Stack S)
/* This performs variant 1 of Shank's log_b(a),
 * See http://www.maths.uq.edu.au/~krm/log.pdf
 */
{
	int tt;
	USL e, i, j, n, nn, r, rr, s, ss;
	MPI *A, *B, *C, *AA, *BB, *CC, *TMP1, *TMP2;
	MPI *PN, *QP, *QN, *PP, *X, *Y, *Z, *QNN, *PNN;
	char buff[20];
        FILE *outfile;

	MPI *XX = stackPop(S);
	MPI *YY = stackPop(S);
	MPI *D = stackPop(S);
	MPI *R = stackPop(S);
	MPI *RR = stackPop(S);
	MPI *E = stackPop(S);
	i = 0;
	j = 0;
	e = CONVERTI(E);
	FREEMPI(E);
	PN = ONEI();
	QP = ONEI();
	QN = ZEROI();
	PP = ZEROI();
	PNN = ONEI();
	QNN = ZEROI();

	s = 2;
	ss = 2;
	r = CONVERTI(R);
	rr = CONVERTI(RR);
	FREEMPI(R);
	FREEMPI(RR);
	if(r >= rr){
		FREEMPI(XX);
		FREEMPI(YY);
		FREEMPI(D);
		FREEMPI(E);
		execerror("arg 4 >= arg5", "");
	}
	C = POWERI(D, (USI)r);
	A = MULTI(XX, C);
	B = MULTI(YY, C);
	CC = POWERI(D, (USI)rr);
	AA = MULTI(XX, CC);
	BB = MULTI(YY, CC);
	FREEMPI(D);
	strcpy(buff, "log.out");
	outfile = fopen(buff, "w");

	while (RSV(B, C) > 0 && RSV(BB, CC) > 0){
		if (e){
			FREEMPI(PNN);
			FREEMPI(QNN);
			PNN = COPYI(PN);
			QNN = COPYI(QN);
		}
		while (RSV(A, B) >= 0){
			TMP1 = MULTI(A, C);
			TMP2 = A;
			A = INT(TMP1, B);
			FREEMPI(TMP1);
			FREEMPI(TMP2);
			i++;
			if (e){
				TMP1 = PN;
				PN = ADDI(PN, PP);
				FREEMPI(TMP1);
				TMP1 = QN;
				QN = ADDI(QN, QP);
				FREEMPI(TMP1);
			}
		}
		n = i;
		s++;
		TMP1 = B;
		B = A;
		A = TMP1;
		i = 0; /* reset the partial quotient count */

		while (RSV(AA, BB) >= 0){
			TMP1 = MULTI(AA, CC);
			TMP2 = AA;
			AA = INT(TMP1, BB);
			FREEMPI(TMP1);
			FREEMPI(TMP2);
			j++;
		}
		nn = j;
		ss++;
		TMP1 = BB;
		BB = AA;
		AA = TMP1;
		j = 0; /* reset the partial quotient count */
		if (n != nn){
		/*	printf("%lu != %lu\n", n, nn);
			printf("breaking\n");*/
			break;
		}
		if (e){
			printf("n[%lu]=%lu:", s - 3, n);
			fprintf(outfile, "n[%lu]=%lu:", s - 3, n);
		}
		else{
			printf("n[%lu]=%lu\n", s - 3, n);
			fprintf(outfile, "n[%lu]=%lu\n", s - 3, n);
		}
		if (e){
			TMP1 = PN;
			PN = PP;
			PP = TMP1;
			TMP1 = QN;
			QN = QP;
			QP = TMP1;
			if(EQUALI(B, C)){
     /* tricky-needed if B=C is reached and n !=nn is not encountered */
				FREEMPI(PNN);
				FREEMPI(QNN);
				PNN = COPYI(PN);
				QNN = COPYI(QN);
			}
			printf("p[%lu]/q[%lu]=", s - 3, s - 3);
			PRINTI(QP);
			printf("/");
			PRINTI(PP);
			fprintf(outfile, "p[%lu]/q[%lu]=", s - 3, s - 3);
			FPRINTI(outfile, QP);
			fprintf(outfile, "/");
			FPRINTI(outfile, PP);
			printf("\n");
			fprintf(outfile, "\n");
		}
	}
	if(e){
		printf("The log of "); 
		fprintf(outfile, "The log of "); 
		PRINTI(XX); 
		FPRINTI(outfile, XX); 
		printf(" to base "); 
		fprintf(outfile, " to base "); 
		PRINTI(YY); 
		FPRINTI(outfile, YY); 
		FREEMPI(XX);
		FREEMPI(YY);
		if (s == 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);
			FREEMPI(PNN);
			FREEMPI(QNN);
			fclose(outfile);
			return;
		}
		printf(" equals ");
		fprintf(outfile, " equals ");
		if (PN->S){
			fclose(outfile);
			tt = COMPARE_DIGITS(QP, PP, QNN, PNN, 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(CC);
				FREEMPI(PN);
				FREEMPI(PP);
				FREEMPI(QP);
				FREEMPI(QN);
				FREEMPI(QNN);
				FREEMPI(PNN);
				fclose(outfile);
				return;
			}
			else{
				printf(" truncated to %u decimal places\n", tt);
				fprintf(outfile, " truncated to %u decimal places\n", tt);
			}
		}
	}

	FREEMPI(A);
	FREEMPI(B);
	FREEMPI(AA);
	FREEMPI(BB);
	FREEMPI(C);
	FREEMPI(CC);
	FREEMPI(PN);
	FREEMPI(PP);
	FREEMPI(QP);
	FREEMPI(QN);
	FREEMPI(XX);
	FREEMPI(YY);
	FREEMPI(QNN);
	FREEMPI(PNN);
	fclose(outfile);

	return;
}

void SHANKSLOG(Stack s)
/* This performs Shank's log_b(a),
 * See http://www.maths.uq.edu.au/~krm/log.pdf
 * See D. Shanks, "A logarithm algorithm", Math. Tables and Other Aids to
 * Computation 8, (1954). 60--64.
 * We print the a[i] for i <= n, to d decimal places.
 */
{
	USL i, n;
	USI j, d;
	MPR *AA, *BB, *TMP1;

	MPI *A = stackPop(s);
	MPI *B = stackPop(s);
	MPI *D = stackPop(s);
	MPI *N = stackPop(s);
	i = 0;
	j = 2;
	d = (USI)CONVERTI(D);
	FREEMPI(D);
	n = CONVERTI(N);
	FREEMPI(N);
	printf("a[0] = "); 
	PRINTI(A); 
	printf("\n");
	printf("a[1] = "); 
	PRINTI(B); 
	printf("\n");
	AA = BUILDMPR();
	AA->N = COPYI(A);
	AA->D = ONEI();
	FREEMPI(A);
	BB = BUILDMPR();
	BB->N = COPYI(B);
	BB->D = ONEI();
	FREEMPI(B);
	while (j <= n){
			printf("Starting AA ");
			PRINTDR(d, AA);
			printf("; BB to ");
			PRINTDR(d, BB);
			printf("\n");
		while (COMPARER(AA, BB) >=0){
			TMP1 = AA;
			AA = RATIOR(AA, BB);
			printf("decreasing AA to ");
			PRINTDR(d, AA);
			printf("; BB=");
			PRINTDR(d, BB);
			printf("\n");
			FREEMPR(TMP1);
			i++;
		}
		printf("n[%u]=%lu, A[%u]=", j - 2, i, j);
		PRINTDR(d, AA);
		printf("\n");
		j++;
		TMP1 = BB;
		BB = AA;
		AA = TMP1;
		i = 0; /* reset the partial quotient count */
	}
	FREEMPR(AA);
	FREEMPR(BB);
	return;
}

void LOG(MPI *A, MPI *B, MPI *D, MPI *R, MPIA *M, MPI **L)
/* Returns an array M[] of L positive integers that are hopefully
 * partial quotients of log(A)/log(B), using C=D^R.
 * Here A > B > 1, D > 1, R >= 1
 * Uses an algorithm in http://www.numbertheory.org/pdfs/log.pdf
 */
{
	USL i, ii, r, s;
	USI l, flag=1;
	MPI *AA, *BB, *AAA, *BBB, *C, *TMP1, *TMP2, *I;
	MPIA P, Q;
	char buff[20];
        FILE *outfile;
	P= NULL;
	Q= NULL;

	s = 0;
	r = CONVERTI(R);
	C = POWERI(D, r);
	AA = MULTI(A, C);
	BB = MULTI(B, C);
	AAA = COPYI(AA);
	BBB = COPYI(BB);
	*M = BUILDMPIA();
	i = ii = 0;
	strcpy(buff, "log.out");
	outfile = fopen(buff, "w");
	fprintf(outfile, "Output from log(");
	FPRINTI(outfile, A);
	fprintf(outfile, ",");
	FPRINTI(outfile, B);
	fprintf(outfile, ",");
	FPRINTI(outfile, D);
	fprintf(outfile, ",");
	FPRINTI(outfile, R);
	fprintf(outfile, ")\n");
       
	while (RSV(BB, C) > 0 && RSV(BBB, C) > 0){
		i = 0;
		while (RSV(AA, BB) >= 0){
			TMP1 = MULTI(AA, C);
			TMP2 = AA;
			AA = INT(TMP1, BB);
			FREEMPI(TMP1);
			FREEMPI(TMP2);
			i++;
                }
		TMP1 = BB;
		BB = AA;
		AA = TMP1;
		ii = 0;
		while (RSV(AAA, BBB) >= 0){
			TMP1 = MULTI(AAA, C);
			TMP2 = AAA;
			AAA = CEILINGI(TMP1, BBB);
			FREEMPI(TMP1);
			FREEMPI(TMP2);
			ii++;
                }
		TMP1 = BBB;
		BBB = AAA;
		AAA = TMP1;
		if(!EQUALI(BB, BBB))
			flag = 0;
		if (i == ii){
			I = CHANGEI(i);
			ADD_TO_MPIA(*M, I, s);
		/*	fprintf(outfile, "m[%lu]=%lu\n", s, i);*/
			fprintf(outfile, "%lu", i); /*for use in Latexing*/
	if((s+1)%20!=0)
			fprintf(outfile, "&");
	else
			fprintf(outfile, "\\\\\n");
		
			FREEMPI(I);
		}
		else{
			printf("M[%lu]=%lu != %lu=M'[%lu]\n", s,i,ii, s);
			fprintf(outfile, "M[%lu]=%lu != %lu=M'[%lu]\n", s,i,ii, s);
			flag = 0;
			break;
		}
		s++;
	}
	*L = CHANGE(s);
	if(flag && EQUALI(BB, C) && EQUALI(BBB, C)){
		printf("log("); PRINTI(A); printf(")");
                printf("/log("); PRINTI(B); printf(")"); 
                printf(" is likely to be the rational number ");
		l = (*M)->size;
		CONVERGENTS(*M, &P, &Q);
		PRINTI(P->A[l-1]); printf("/"); PRINTI(Q->A[l-1]);
		printf("\n");
		FREEMPIA(P);
		FREEMPIA(Q);
	}
	FREEMPI(C);
	FREEMPI(AA);
	FREEMPI(AAA);
	FREEMPI(BB);
	FREEMPI(BBB);
	fclose(outfile);
}

MPI *LOGX(MPI *A, MPI *B, MPI *D, MPI *R, MPIA *M, MPI **L)
{
	MPI *TMP = ONEI();
	int s;

	s = COMPAREI(B, TMP);
	if (s <= 0){
		printf("B <= 1\n");
		FREEMPI(TMP);
		M = NULL;
		return NULL;
	}
	s = COMPAREI(D, TMP);
	if (s <= 0){
		printf("D <= 1\n");
		FREEMPI(TMP);
		M = NULL;
		return NULL;
	}
	s = COMPAREI(A, B);
	if (s <= 0){
		printf("A <= B\n");
		FREEMPI(TMP);
		M = NULL;
		return NULL;
	}
	s = COMPAREI(R, TMP);
	if (s < 0){
		printf("R < 1\n");
		FREEMPI(TMP);
		M = NULL;
		return NULL;
	}
	FREEMPI(TMP);
	
        LOG(A, B, D, R, M, L);
	return ONEI();
}

void TESTLOG1(MPI *A, MPI *B, MPI *D, MPI *R)
/* This is similar to LOG, but is for use in TESTLOG2 */
{
	USL i, ii, r, s;
	MPI *AA, *BB, *AAA, *BBB, *C, *TMP1, *TMP2;
	char buff[20];
        FILE *outfile;

	s = 0;
	r = CONVERTI(R);
	C = POWERI(D, r);
	AA = MULTI(A, C);
	BB = MULTI(B, C);
	AAA = COPYI(AA);
	BBB = COPYI(BB);
	i = ii = 0;
	strcpy(buff, "testlog.out");
	outfile = fopen(buff, "a");
	while (RSV(BB, C) > 0 && RSV(BBB, C) > 0){
		i = 0;
		while (RSV(AA, BB) >= 0){
			TMP1 = MULTI(AA, C);
			TMP2 = AA;
			AA = INT(TMP1, BB);
			FREEMPI(TMP1);
			FREEMPI(TMP2);
			i++;
                }
		TMP1 = BB;
		BB = AA;
		AA = TMP1;
		ii = 0;
		while (RSV(AAA, BBB) >= 0){
			TMP1 = MULTI(AAA, C);
			TMP2 = AAA;
			AAA = CEILINGI(TMP1, BBB);
			FREEMPI(TMP1);
			FREEMPI(TMP2);
			ii++;
                }
		TMP1 = BBB;
		BBB = AAA;
		AAA = TMP1;
		if (i == ii){
			printf("%lu,", i);
			fflush(stdout);
			fprintf(outfile, "%lu,", i);
		}
		else
			break;
		s++;
	}
	printf("\n");
	fprintf(outfile,"\n");
	FREEMPI(C);
	FREEMPI(AA);
	FREEMPI(AAA);
	FREEMPI(BB);
	FREEMPI(BBB);
	fclose(outfile);
}

void TESTLOG(MPI *A, MPI *B, MPI *D, MPI *M, MPI *N)
/* runs TESTLOG1(A,B,D,R) for R=M,...,N. */
{
	USL r, m, n;
	MPI *R;
	char buff[20];
        FILE *outfile;

	strcpy(buff, "testlog.out");
	outfile = fopen(buff, "w");
	m = CONVERTI(M);
	n = CONVERTI(N);
	for (r = m; r <= n; r++){
		R = CHANGE(r);
		TESTLOG1(A, B, D, R);
		FREEMPI(R);
	}
	fclose(outfile);
}

MPI *TESTLOGX(MPI *A, MPI *B, MPI *D, MPI *M, MPI *N)
{
	MPI *TMP = ONEI();
	int s;

	s = COMPAREI(B, TMP);
	if (s <= 0){
		printf("B <= 1\n");
		FREEMPI(TMP);
		M = NULL;
		return NULL;
	}
	s = COMPAREI(D, TMP);
	if (s <= 0){
		printf("D <= 1\n");
		FREEMPI(TMP);
		M = NULL;
		return NULL;
	}
	s = COMPAREI(A, B);
	if (s <= 0){
		printf("A <= B\n");
		FREEMPI(TMP);
		M = NULL;
		return NULL;
	}
	s = COMPAREI(M, TMP);
	if (s < 0){
		printf("M < 1\n");
		FREEMPI(TMP);
		M = NULL;
		return NULL;
	}
	s = COMPAREI(M, TMP);
	if (s < 0){
		printf("M < 1\n");
		FREEMPI(TMP);
		M = NULL;
		return NULL;
	}
	s = COMPAREI(N, M);
	if (s < 0){
		printf("N < M\n");
		FREEMPI(TMP);
		M = NULL;
		return NULL;
	}
	FREEMPI(TMP);
	
        TESTLOG(A, B, D, M, N);
	return ONEI();
}