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

/* i8.c */ 
#include  
#include "integer.h" 
#include "fun.h" 
 
MPI *APPROXN(MPI *Aptr, unsigned int m) 
/* 
 *  *Eptr = 2 ^ (1 + [i / m]), (where i + 1 is the number of binary digits in  
the positive  MPI *Aptr) is the initial approximation to the integer part of 
the m-th root of *Aptr.  (*Aptr)^(1/m) < *Eptr <= 2*(*Aptr)^(1/m). 
 */ 
{ 
	MPI *Eptr; 
	unsigned int i = 0; 
 
	i = BINARYB(Aptr) - 1; 
	Eptr = POWER_I(2L, 1 + (i / m)); 
	return (Eptr); 
} 
 
MPI *NEWTON(MPI *Aptr, unsigned int m, MPI *Fptr) 
/* 
 * Returns the integer part of the mth root of *Aptr, where m is 
 * an unsigned integer, 1 < m < R0.  
 * The method is Newton's, using the integer part function, with initial 
 * approximation *Fptr > (*Aptr)^(1/m). 
 */ 
{ 
	MPI *X, *Y, *Z, *Temp, *Temp1; 
 
        X = COPYI(Fptr); 
	while (1) 
	{ 
		Y = COPYI(X); 
		Z = POWERI(X, m - 1); 
		Temp = INT0(Aptr, Z); 
		FREEMPI(Z); 
		Z = MULT_I(X, m - 1); 
		Temp1 = ADD0I(Temp, Z); 
		FREEMPI(X); 
		FREEMPI(Temp); 
		X = INT0_(Temp1, m); 
		FREEMPI(Temp1); 
		FREEMPI(Z); 
		if (RSV(Y, X) != 1) 
		{ 
			FREEMPI(X); 
			return (Y); 
		} 
		FREEMPI(Y); 
	} 
} 
 
MPI *BABY_MTHROOT(MPI *Aptr, unsigned int m) 
/* 
 * *Eptr is the integer part of the mth root of *Aptr, where m is 
 * an unsigned integer, 1 < m < R0. 
 * The method is Newton's, using the integer part function, with initial 
 * approximation provided by APPROXN(), for "small" *Aptr. 
 */ 
{ 
	MPI *X, *Eptr; 
	 
	X = APPROXN(Aptr, m); 
	Eptr = NEWTON(Aptr, m, X); 
	FREEMPI(X); 
	return (Eptr); 
} 
 
MPI *BIG_MTHROOT(MPI *Aptr, unsigned int m) 
/* 
 * *Eptr, the integer part of the mth root of the positive MPI *Aptr, 
 * 1 < m < R0, is obtained by Newton's method, using the integer part function. 
 */ 
{ 
	MPI *a, *X, *Eptr, *Temp; 
	unsigned int i, n, r, t; 
	 
	n = Aptr->D; 
	if (n < m) 
	{ 
		Eptr = BABY_MTHROOT(Aptr, m); 
		return (Eptr); 
	} 
	else 
	{ 
		r = n / m; 
		a = BUILDMPI((n % m) + 1); /* a = [*Aptr/R0^(m*r)] */ 
		a->S = 1; 
		t = r * m; 
		for (i = 0; i <= a->D; i++)  
			a->V[i] = Aptr->V[i + t]; 
		X = BABY_MTHROOT(a, m); 
		FREEMPI(a); 
		Temp = X; 
		X = ADD0_I(X, (USL)1); 
		FREEMPI(Temp); 
		Temp = X; 
		X = BUILDMPI(Temp->D + r + 1); /* X = X* R0^r */ 
		X->S = 1; 
		for (i = X->D; i >= r; i--) 
			X->V[i] = Temp->V[i - r]; 
		FREEMPI(Temp); 
		INTSETUL(X->V, (USL)r, (USL)0); 	 
		/* sets the first r array elements to 0 */ 
 
		Eptr = NEWTON(Aptr, m, X); /* X is the initial approximation */ 
		FREEMPI(X); 
		return (Eptr); 
	} 
} 
 
void MTHROOT(MPI *Aptr, MPI *Bptr, unsigned int m, unsigned int r) 
/* 
 * *Aptr and *Bptr are positive MPI'S.  
 * the mthroot of *Aptr/(*Bptr) is computed to r decimal places, r >= 0;  
 * m, r are integers, 0 < m * r < R0 * R0.  
 */ 
{ 
	MPI *D, *E, *F, *G, *Y, *Temp; 
	unsigned int i, l; 
	 
	E = POWER_I(10L, r * m); 
	Temp = E; 
	E = MULTI(E, Aptr); 
	FREEMPI(Temp); 
	Temp = E; 
	E = INT0(E, Bptr); 
	FREEMPI(Temp); 
	if (E->S == 0) 
	{ 
		FREEMPI(E); 
		if (r == 0) 
		{ 
			printf("%uthroot of ", m); 
			PRINTI(Aptr); 
			printf("/"); 
			PRINTI(Bptr); 
			printf(" = 0\n"); 
			return; 
		} 
		else 
		{ 
			printf("%uthroot of ", m); 
			PRINTI(Aptr); 
			printf("/"); 
			PRINTI(Bptr); 
			printf(" = 0."); 
			for (i = 1; i <= r; i++) 
				printf("0"); 
			printf("\n"); 
			return; 
		} 
	} 
	Y = BIG_MTHROOT(E, m); 
	FREEMPI(E); 
	F = POWER_I(10L, r); 
	G = MOD0(Y, F); 
	D = INT0(Y, F); 
	FREEMPI(F); 
	FREEMPI(Y); 
	printf("%uthroot of ", m); 
	PRINTI(Aptr); 
	printf("/"); 
	PRINTI(Bptr); 
	printf(" = "); 
	PRINTI(D); 
	FREEMPI(D); 
	if (r == 0) 
	{ 
		FREEMPI(G); 
		printf("\n"); 
		return; 
	} 
	printf("."); 
	l = LENGTHI(G); /* the number of decimal digits in G */ 
	for (i = 1; i <= r - l; i++) 
		printf("0"); 
	PRINTI(G); 
	printf("\n"); 
	FREEMPI(G); 
	return; 
} 
 
MPR *MTHROOTR(MPR *Nptr, unsigned int m, unsigned int r) 
/* 
 * *Nptr is a positive MPR. 
 * the mthroot of *Nptr is computed to r decimal places, r >= 0 . 
 * m, r are integers, 0 < m * r < R0 * R0 . 
 */ 
{ 
	MPI *Tmp0I, *Tmp1I, *Tmp2I; 
	MPR *Tmp1R; 
	 
	Tmp0I = POWER_I(10L, r); 
	Tmp1I = POWERI(Tmp0I, m); 
	Tmp2I = MULTI(Tmp1I, Nptr->N); 
	Tmp1R = RATIOI(Tmp2I, Nptr->D); 
	FREEMPI(Tmp1I); 
	FREEMPI(Tmp2I); 
	Tmp1I = INTI(Tmp1R->N, Tmp1R->D); 
	FREEMPR(Tmp1R); 
	if (EQZEROI(Tmp1I)) 
	{ 
		FREEMPI(Tmp1I); 
		FREEMPI(Tmp0I); 
		return (ZEROR()); 
	} 
	Tmp2I = BIG_MTHROOT(Tmp1I, m); 
	FREEMPI(Tmp1I); 
	Tmp1R = RATIOI(Tmp2I, Tmp0I); 
	FREEMPI(Tmp0I); 
	FREEMPI(Tmp2I); 
	return (Tmp1R); 
} 
 
MPI *PI(unsigned int r) 
{ 
	MPI *PI, *X, *X1, *Y, *U, *V; 
	MPI *Tmp, *Tmp0I, *Tmp1I, *Tmp2I, *Tmp3I, *Tmp4I, *Tmp5I; 
 
	Tmp0I = POWER_I(10L, r); 
	Tmp1I = POWERI(Tmp0I, 2);/* 10^2r */ 
	Tmp2I = POWERI(Tmp1I, 2);/* 10^4r */ 
	Tmp3I = POWERI(Tmp2I, 2);/* 10^8r */ 
	Tmp = MULT_I(Tmp2I, 2L); 
	X = BIG_MTHROOT(Tmp, 2); 
	FREEMPI(Tmp); 
	Tmp4I = MULT_I(Tmp1I, 2L); 
	PI = ADD0I(Tmp4I, X); 
	FREEMPI(Tmp4I); 
	Tmp4I = MULT_I(Tmp3I, 2L); 
	Y = BIG_MTHROOT(Tmp4I, 4); 
	FREEMPI(Tmp3I); 
	FREEMPI(Tmp4I); 
	while (1) 
	{ 
		Tmp4I = ADD0I(X, Tmp1I); 
		Tmp = X; 
		X = MULTI(Tmp0I, Tmp4I); 
		FREEMPI(Tmp); 
		Tmp3I = BIG_MTHROOT(X, 2); 
		X1 = MULT_I(Tmp3I, 2L); 
		Tmp = X; 
		X = INT0(X, X1); 
		FREEMPI(Tmp); 
		FREEMPI(X1); 
 
		U = MULTI(X, Y); 
		Tmp = U; 
		U = ADD0I(U, Tmp2I); 
		printf(" U = "); 
		PRINTI(U); 
		printf("\n"); 
		FREEMPI(Tmp); 
		Tmp5I = ADD0I(Y, Tmp1I); 
		V = MULTI(Tmp3I, Tmp5I); 
		Tmp = V; 
		V = MULTI(V, Tmp0I); 
		FREEMPI(Tmp); 
		Tmp = Y; 
		Y = INT0(U, V); 
		printf(" V = "); 
		PRINTI(V); 
		printf("\n"); 
		FREEMPI(U); 
		FREEMPI(V); 
		printf(" Y = "); 
		PRINTI(Y); 
		printf("\n"); 
 
		if (EQUALI(Tmp1I, Y)) 
		{ 
			FREEMPI(Tmp0I); 
			FREEMPI(Tmp1I); 
			FREEMPI(Tmp2I); 
			FREEMPI(Tmp3I); 
			FREEMPI(Tmp4I); 
			FREEMPI(Tmp5I); 
			return (PI); 
		} 
		Tmp = PI; 
		PI = MULTI(PI, Tmp4I); 
		FREEMPI(Tmp); 
		Tmp = PI; 
		PI = INT0(PI, Tmp5I); 
		FREEMPI(Tmp4I); 
		FREEMPI(Tmp5I); 
	} 
} 
 
MPR *PII(unsigned int r) 
/* 
 * From " A very rapidly convergent product expansion for pi," Bit 23 (1983)  
 * 538-540, by J.M. Borwein and P.B. Borwein. The algorithm delivers pi  
 * to 2^r decimals. 
 */ 
{ 
	MPR * Tmp1, * Tmp2, * Tmp3, * Tmp4, * Tmp5, * Tmp6, * Tmp7, *Tmp8; 
	MPR *X, *Y, *PI, *ONE, *TWO, *Tmp, *TEST, *SS, *TEN, *Tmp9; 
	MPI *S, *TmpI; 
	unsigned int s; 
 
	S = POWER_I(2L, r); 
	s = CONVERTI(S); 
	FREEMPI(S); 
	SS = BUILDMPR(); 
	SS->N = POWER_I(10L, s + 1); 
	SS->D = ONEI(); 
	ONE = ONER(); 
	TWO = BUILDMPR(); 
	TWO->N = CHANGE((USL)2); 
	TWO->D = ONEI(); 
	TEN = BUILDMPR(); 
	TEN->N = CHANGE((USL)10); 
	TEN->D = ONEI(); 
	X = MTHROOTR(TWO, 2, s + 1); 
	PI = ADDR(TWO, X); 
	Tmp1 = MTHROOTR(X, 2, s + 1); 
	Tmp2 = INVERSER(Tmp1); 
	Tmp3 = ADDR(Tmp1, Tmp2); 
	Tmp = X; 
	X = RATIOR(Tmp3, TWO); 
	printf("X = "); 
	PRINTDR(s, X); 
	printf("\n"); 
	FREEMPR(Tmp); 
	FREEMPR(Tmp1); 
	FREEMPR(Tmp2); 
	FREEMPR(Tmp3); 
	Y = MTHROOTR(TWO, 4, s + 1); 
	while (1) 
	{ 
		Tmp1 = MTHROOTR(X, 2, s + 1); 
		Tmp2 = INVERSER(Tmp1); 
		Tmp3 = ADDR(Tmp1, Tmp2); 
		Tmp4 = RATIOR(Tmp3, TWO); 
		Tmp5 = MULTR(Y, Tmp1); 
		Tmp6 = ADDR(Tmp5, Tmp2); 
		Tmp7 = ADDR(Y, ONE); 
		Tmp8 = ADDR(X, ONE); 
		Tmp = PI; 
		PI = MULTR(Tmp8, PI); 
		FREEMPR(Tmp); 
		Tmp = PI; 
		PI = RATIOR(PI, Tmp7); 
		FREEMPR(Tmp); 
		Tmp = X; 
		X = Tmp4; 
		FREEMPR(Tmp); 
		printf("X = "); 
		PRINTDR(s, X); 
		printf("\n"); 
		Tmp = Y; 
		Y = RATIOR(Tmp6, Tmp7); 
		FREEMPR(Tmp); 
		printf("Y = "); 
		PRINTDR(s, Y); 
		printf("\n"); 
		TEST = SUBR(Y, ONE); 
		Tmp = TEST; 
		Tmp9 = MULTR(SS, TEN); 
		TEST = MULTR(Tmp9, TEST); 
		FREEMPR(Tmp); 
		TmpI = INT0(TEST->N, TEST->D); 
		FREEMPR(TEST); 
		if (TmpI->S == 0) 
		{ 
			printf("Y - 1 = "); 
			PRINTDR(s, Y); 
			printf("\n"); 
			FREEMPR(Tmp1); 
			FREEMPR(Tmp2); 
			FREEMPR(Tmp3); 
			FREEMPR(Tmp5); 
			FREEMPR(Tmp6); 
			FREEMPR(Tmp7); 
			FREEMPR(Tmp8); 
			FREEMPR(Tmp9); 
			FREEMPI(TmpI); 
			FREEMPR(ONE); 
			FREEMPR(TWO); 
			FREEMPR(TEN); 
			FREEMPR(SS); 
			FREEMPR(X); 
			FREEMPR(Y); 
			return (PI); 
		} 
		FREEMPR(Tmp1); 
		FREEMPR(Tmp2); 
		FREEMPR(Tmp3); 
		FREEMPR(Tmp5); 
		FREEMPR(Tmp6); 
		FREEMPR(Tmp7); 
		FREEMPR(Tmp8); 
		FREEMPR(Tmp9); 
		FREEMPI(TmpI); 
	} 
}