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

/* program:  "i2.c" */ 
/* 
 * A translation into C of the long division section of pascal program 
 * ARITHMETIC,  from "SCIENTIFIC PASCAL" by H. FLANDERS pp. 342-357. 
 */ 
 
#include  
#include  
#include "integer.h" 
#include "fun.h" 
	 
unsigned long l0, n0, bn, bn1; 
 
void INIT(MPI *Aptr, MPI *Bptr, MPI **A1ptr, MPI **B1ptr) 		 
/* 
 * see lemmas 2 and 3, pp. 349-350.  
 * input: Aptr, Bptr are pointers to MPI'S, Aptr->D >=  Bptr->D = n0 >= 1. 
 *  output: l0 = int[R0 / (Bptr->V[n0] + 1), *A1ptr = l0 * *Aptr, 
 * *B1ptr = l0 * *Bptr, bn = B1ptr->V[n0], bn1 = B1ptr->V[n0 - 1]. 
 */ 
{ 
	n0 = Bptr->D; 
	l0 = R0 / (Bptr->V[n0] + 1); 
	if (l0 > 1) 
	{ 
		*A1ptr = MULT_I(Aptr,(long)(l0)); 
		*B1ptr = MULT_I(Bptr,(long)(l0)); 
	} 
	else 
	{ 
		*A1ptr = COPYI(Aptr); 
		*B1ptr = COPYI(Bptr); 
	} 
	bn = (*B1ptr)->V[n0]; 
	bn1 = (*B1ptr)->V[n0 - 1]; 
	return; 
} 
 
unsigned long qp; 
void FIND(unsigned int k, MPI *A1ptr, MPI *B1ptr, unsigned long R[]) 
/* 
 * finds the k-th trial quotient digit qp = Q' using lemmas 1-4.  
 */ 
{ 
	unsigned long t, u, rp; 
 
	if (k == A1ptr->D - B1ptr->D) 
		R[A1ptr->D + 1] = 0; /* an extra 0 is added to the front of 
			the initial dividend to ensure qp < R0.(see p348) */ 
	 
	u = (R[k + B1ptr->D + 1] << T0) + R[k + B1ptr->D]; 
	t = u / bn; 
	if (t < R0) 
		qp = t; 
	else 
		qp = R0 - 1; 
 
	rp = u - qp * bn; /* Flanders erroneously has t instead of qp */  
	 
	if (rp >= R0) /* rp * R0 is then >= 2 ^ O32 > bn1 * qp. */ 
		return; 
	while (bn1 * qp > (rp << T0) + R[k + n0 - 1]) 
	{	/* see lemma 4(1) */  
		qp--; 
		rp = rp + bn; 
		if (rp >= R0) /* rp * R0 is then >= 2 ^ O32 > bn1 * qp. */ 
			return; 
	} 
	return; 
} 
 
void NEXT(int k, MPI *A1ptr, MPI *B1ptr, unsigned long R[]) 
/* 
 * the run of n0 + 1 "digits" in the current remainder, from 
 * the (k + n0 + 1)-th down to the k-th, is divided by *B1ptr. The quotient 
 * becomes the k-th "digit" of *Qptr in INT0() and the remainder R is updated. 
 */ 
{ 
	long l; 
	unsigned int j; 
	unsigned long t, ct = 0, cr = 1;	/* initialize "carries */ 
	 
	for (j = 0; j <= n0; j++) 
	{	/* multiply with carry */ 
		t = (B1ptr->V[j]) * qp + ct; 
		ct = t >> T0; 
		t = t & (R0 - 1); 	/* digit to subract with carry */ 
		t = R[k + j] - t + R0 + cr - 1; 
		R[k + j] = t & (R0 - 1); 
		cr = t >> T0; 
	} 
	if (k == A1ptr->D - B1ptr->D) 
		R[A1ptr->D + 1] = 0; /* avoids garbage value for R[k + n0 + 1] */ 
	l = R[k + n0 + 1] - ct + cr - 1; 
	if (l < 0) { /* qp is one too large */ 
		qp--; 
		cr = 0;  
		for (j = 0; j <= n0; j++) 
		{ /* add *B1ptr back */ 
			t = R[k + j] + B1ptr->V[j] + cr; 
			R[k + j] = t & (R0 - 1); 
			cr = t >> T0; 
		} 
		R[k + n0 + 1] = l + cr; 
	} 
	return; 
	 
} 
 
 
MPI *INT0_(MPI *Aptr, unsigned long b) 
/* 
 * input: Aptr is a pointer to an non-negative MPI, b is a positive integer, 
 * b < R0. output:  *Qptr = int(*Aptr / b). 
 */ 
{ 
	unsigned int k, n; 
	unsigned long t; 
	MPI *Qptr, *tmp; 
 
	if (b >= R0) 
	{ 
		fprintf(stderr, "in INT0_I, %lu >= R0\n", b); 
		exit(1); 
	} 
	Qptr = COPYI(Aptr); 
	n = Qptr->D; 
	t = Qptr->V[Qptr->D]; 
	if (Qptr->D != 0) 
	{ 
		for (k = Qptr->D; k >= 1; k--) 
		{ 
			Qptr->V[k] = t / b; 
			t  = Qptr->V[k - 1] + ((t % b) << T0); 
		} 
		if (Qptr->V[Qptr->D] == 0) 
		{ 
			 
			tmp = BANK_REALLOC(Qptr, n - 1); 
			FREEMPI(Qptr); 
			Qptr = tmp; 
			/* we remark that Qptr->S = 1 here. */ 
		} 
	} 	 
	Qptr->V[0] = t / b; 
	if (Qptr->D == 0)  
		Qptr->S = Qptr->V[0] ? 1 : 0; 
	return Qptr; 
} 
 
unsigned long MOD0_(MPI *Aptr, unsigned long b) 
/* 
 * input: Aptr is a pointer to a non-negative MPI, b is a positive integer, 
 * b < R0.  output: *Aptr mod b. 
 */ 
{ 
	MPI *R; 
	unsigned int k; 
	unsigned long t; 
	if (b >= R0) 
	{ 
		fprintf(stderr, "in MOD0_, %lu >= R0\n", b); 
		exit(1); 
	} 
	R = COPYI(Aptr); 
	t = R->V[R->D]; 
	if (R->D != 0) 
	{ 
		for (k = R->D; k >= 1; k--) 
		{ 
			R->V[k] = t / b; 
			t  = R->V[k - 1] + ((t % b) << T0); 
		} 
	} 	 
 
	FREEMPI(R); 
	return (t % b); 
} 
 
MPI *INT_(MPI *Aptr, unsigned long b) 
/* 
 * input: Aptr is a pointer to an MPI, b is a positive integer, b < R0. 
 * output:  *Qptr = int(*Aptr, b). 
 */ 
{ 
	MPI *Q, *Qptr; 
	unsigned long r; 
 
	if (b >= R0) 
	{ 
		fprintf(stderr, "in INT_, %lu >= R0\n", b); 
		exit(1); 
	} 
	if (Aptr->S >= 0) 
		return (INT0_(Aptr, b)); 
	else 
	{ 
		Aptr->S =  1; 
		Q = INT0_(Aptr, b); 
		r = MOD0_(Aptr, b); 
		Aptr->S = -1; 
		/* -*Aptr = *Q * b + r, 0 <= r < b, 
				*Aptr = (-*Q - 1) * b + b - r */ 
		if (r) 
		{/* add 1 to |*Q|, but not changing the sign of *Q */ 
			Qptr = ADD0_I(Q, (USL)1); 
			Qptr->S = -1; 
			FREEMPI(Q); 
			return Qptr; 
		} 
		else 
		{					/* r = 0 here */ 
			Q->S = -1; 
			Qptr = COPYI(Q); 
			FREEMPI(Q); 
			return Qptr;				 
		} 
	} 
} 
 
unsigned long MOD_(MPI *Aptr, unsigned long b) 
/* 
 * input: Aptr is a pointer to an MPI, b is a positive integer, b < R0. 
 * output: *Aptr mod b. 
 */ 
{ 
	unsigned long r; 
 
	if (b >= R0) 
	{ 
		fprintf(stderr, "in MOD_, %lu >= R0\n", b); 
		exit(1); 
	} 
	if (Aptr->S >= 0) 
		return (MOD0_(Aptr, b)); 
	else 
	{ 
		Aptr->S =  1; 
		r = MOD0_(Aptr, b); 
		Aptr->S = -1; 
		/* -*Aptr = *Q * b + r, 0 <= r < b, 
				*Aptr = (-*Q - 1) * b + b - r */ 
		if (r)  
			return (b - r); 
		else 
			return (0);				 
	} 
} 
 
MPI *INT0(MPI *Aptr, MPI *Bptr) 
/* 
 * input: Aptr, Bptr, pointers to non-negative MPI'S, *Bptr > 0. 
 * output: *Qptr = int(*Aptr / *Bptr), *Rptr = *Aptr mod *Bptr. 
 */ 
{ 
	MPI *A1, *B1, *Qptr, *tmp; 
	int k; 
	unsigned int e, f, degree; 
	unsigned long *R; 
 
	if (Bptr->D == 0) 
		return (INT0_(Aptr, Bptr->V[0])); 
	if (Aptr->D < Bptr->D) 
		return ZEROI(); 
	INIT(Aptr, Bptr, &A1, &B1); 
	f = 2 + A1->D; 
	R = (unsigned long *)mmalloc((USL)(f * sizeof(unsigned long))); 
	/* we malloc an extra place needed in FIND(). */ 
	for (k = 0; k <= A1->D; k++) 
		R[k] = A1->V[k]; 
	e = A1->D - B1->D; 
	Qptr = BUILDMPI(e + 1); 
	for (k = e; k >= 0; k--) 
	{ 
		FIND((unsigned int) k, A1, B1, R); 
		NEXT(k, A1, B1, R); 
		Qptr->V[k] = qp; 
	} 
	ffree((char *)R, f * sizeof(unsigned long)); 
	Qptr->S = 0;			/* finding Qptr->S and Qptr->D */ 
	degree = 0; 
	k = e; 
	FREEMPI(A1); 
	FREEMPI(B1); 
	while (k >= 0) 
	{ 		   	 
		if (Qptr->V[k] != 0) 
		{ 
			Qptr->S = 1; 
			degree = k; 
			k = -1; 
		} 
		else 
			k--; 
	} 							 
	tmp = Qptr; 
	Qptr = BANK_REALLOC(Qptr, degree); 
	FREEMPI(tmp); 
	return Qptr; 
} 
 
MPI *MOD0(MPI *Aptr, MPI *Bptr) 
/* 
 * input: Aptr, Bptr, pointers to MPI'S, *Aptr >= 0, *Bptr > 0. 
 * output: *Rptr = *Aptr mod *Bptr. 
 */ 
{ 
	MPI *A1, *B1, *Rptr, *Temp; 
	int k, s;  
	unsigned int d, e, f; 
	unsigned long *R; 
 
	if (Bptr->D == 0) 
	{ 
		Rptr = BUILDMPI(1); 
		Rptr->V[0] = MOD0_(Aptr, Bptr->V[0]); 
		Rptr->S = Rptr->V[0] ? 1 : 0; 
 
		return Rptr; 
	} 
	if (Aptr->D < Bptr->D) 
		return COPYI(Aptr); 
	INIT(Aptr, Bptr, &A1, &B1); 
	f = 2 + A1->D; 
	R = (unsigned long *)mmalloc((USL)(f * sizeof(unsigned long))); 
	/* we malloc an extra place needed in FIND(). */ 
	for (k = 0; k <= A1->D; k++) 
		R[k] = A1->V[k]; 
	e = A1->D - B1->D; 
	for (k = e; k >= 0; k--) 
	{ 
		FIND((unsigned int) k, A1, B1, R); 
		NEXT(k, A1, B1, R); 
	} 
	FREEMPI(A1); 
	FREEMPI(B1); 
 
	/* finding Rptr->S and Rptr->D */ 
 
	s = 0; 
	d = 0; 
	k = n0; 
	while (k >= 0) 
	{ 
		if (R[k] != 0) 
		{ 
			s = 1; 
			d = k; 
			k = -1; 
		} 
		else 
			k--; 
	} 
	Rptr = BUILDMPI(1 + d); 
	Rptr->S = s; 
	for (k = 0; k <= Rptr->D; k++) 
		Rptr->V[k] =  R[k]; 
	ffree((char *)R, f*sizeof(unsigned long)); 
	Temp = Rptr; 
	Rptr = INT0_(Temp, l0);			 
	FREEMPI(Temp);	 
 
	return Rptr; 
} 
 
MPI *INT(MPI *Aptr, MPI *Bptr) 
/* 
 * input: Aptr, Bptr, pointers to MPI'S, *Bptr > 0, 
 * output: MPI *Qptr = int(*Aptr, *Bptr). 
 */ 
{ 
	MPI *Q, *R, *Qptr; 
 
	if (Bptr->D == 0) 
		return	INT_(Aptr, Bptr->V[0]); 
	if (Aptr->S >= 0) 
		return INT0(Aptr, Bptr); 
	else 
	{ 
		Aptr->S =  1; 
		Q = INT0(Aptr, Bptr); 
		R = MOD0(Aptr, Bptr); 
		Aptr->S = -1; 
		/* -*Aptr = *Q * *Bptr + *R, 0 <= *R < *Bptr, 
				*Aptr = (-*Q - 1) * *Bptr + *Bptr - R */ 
		if (R->S == 1) 
		{/* add 1 to |*Q|, but not changing the sign of *Q */ 
			Qptr = ADD0_I(Q, (USL)1); 
			Qptr->S = -1; 
		} 
		else 
		{					/* *R = 0 here */ 
			Q->S = -1; 
			Qptr = COPYI(Q); 
		} 
		FREEMPI(Q); 
		FREEMPI(R); 
		return Qptr; 
	} 
} 
 
MPI *INTI(MPI *Aptr, MPI *Bptr) 
/* 
 * input: Aptr, Bptr, pointers to MPI'S, *Bptr != 0, 
 * output: MPI *Qptr = int(*Aptr, *Bptr). 
 */ 
{ 
	MPI *A, *B, *Qptr; 
 
	if (Bptr->S > 0) 
		return INT(Aptr, Bptr); 
	else 
	{ 
			A = COPYI(Aptr); 
		B = COPYI(Bptr); 
		A->S = -(A->S); 
		B->S = -(B->S); 
		Qptr = INT(A, B); 
		FREEMPI(A); 
		FREEMPI(B); 
		return Qptr; 
		} 
} 
 
MPI *MOD(MPI *Aptr, MPI *Bptr) 
/* 
 * input: Aptr, Bptr, pointers to MPI'S, *Bptr > 0. 
 * output: MPI *Rptr = *Aptr mod *Bptr. 
 */ 
{ 
	MPI *R, *Rptr; 
 
	if (Bptr->D == 0) 
		return CHANGE(MOD_(Aptr, Bptr->V[0])); 
	if (Aptr->S >= 0) 
		return 	MOD0(Aptr, Bptr); 
	else 
	{ 
		Aptr->S = 1; 
		R = MOD0(Aptr, Bptr); 
		Aptr->S = -1; 
		/* -*Aptr = *Q * *Bptr + *R, 0 <= *R < *Bptr, 
				*Aptr = (-*Q - 1) * *Bptr + *Bptr - *R */ 
		if (R->S == 1) 
			Rptr = SUB0I(Bptr, R); 
		else 
			Rptr = ZEROI(); 
		FREEMPI(R); 
		return Rptr;				 
	} 
} 
 
unsigned long MODINT0_(MPI *Aptr, unsigned long b, MPI **Qptr) 
/* 
 * input: Aptr is a pointer to an non-negative MPI, b is a positive integer, 
 * b < R0. output: *Aptr (mod b) and *Qptr = int(*Aptr / b). 
 */ 
{ 
	unsigned int k, n; 
	unsigned long t; 
	MPI *tmp; 
 
	if (b >= R0) 
	{ 
		fprintf(stderr, "in MODINT0_, %lu >= R0\n", b); 
		exit(1); 
	} 
	*Qptr = COPYI(Aptr); 
	n = (*Qptr)->D; 
	t = (*Qptr)->V[(*Qptr)->D]; 
	if ((*Qptr)->D != 0) 
	{ 
		for (k = (*Qptr)->D; k >= 1; k--) 
		{ 
			(*Qptr)->V[k] = t / b; 
			t  = (*Qptr)->V[k - 1] + ((t % b) << T0); 
		} 
		if ((*Qptr)->V[(*Qptr)->D] == 0) 
		{ 
			 
			tmp = BANK_REALLOC(*Qptr, n - 1); 
			FREEMPI(*Qptr); 
			*Qptr = tmp; 
			/* we remark that (*Qptr)->S = 1 here. */ 
		} 
	} 	 
	(*Qptr)->V[0] = t / b; 
	if ((*Qptr)->D == 0)  
		(*Qptr)->S = (*Qptr)->V[0] ? 1 : 0; 
	return (t % b); 
} 
 
MPI *HALFMOD(MPI *Aptr, MPI *Pptr) 
/* Returns R=Aptr(mod Pptr) if R <= Pptr/2, otherwise R-Pptr. */ 
{ 
	MPI *T, *U, *W; 
	int t; 
 
	T = MOD(Aptr, Pptr); 
	W = MULT_II(T, 2); 
	t = RSV(W, Pptr); 
	FREEMPI(W); 
	if (t == 1) 
		U = SUBI(T, Pptr); 
	else 
		U = COPYI(T); 
	FREEMPI(T); 
	return (U); 
} 
 
MPI *HALFMODX(MPI *Aptr, MPI *Pptr) 
/* Returns R=Aptr(mod Pptr) if R <= Pptr/2, otherwise R-Pptr. */ 
{ 
	if(Pptr->S <= 0){ 
		printf("N <= 0\n"); 
		return(NULL); 
	} 
	else 
		return(HALFMOD(Aptr, Pptr)); 
} 
 
MPI *CEILINGI(MPI *A, MPI *B) 
/* Returns the least integer not less than A/B. */ 
{ 
	MPI *X, *TMP, *TMP1; 
	USI t; 
 
	X = INTI(A, B); 
	TMP = MULTI(X, B); 
	t = EQUALI(TMP, A); 
	FREEMPI(TMP); 
	if (t) 
		return (X); 
	else{ 
		TMP = ONEI(); 
		TMP1 = ADDI(X, TMP); 
		FREEMPI(TMP); 
		FREEMPI(X); 
		return(TMP1); 
	} 
} 
 
MPI *CEILINGIX(MPI *A, MPI *B) 
{ 
	if(B->S == 0){ 
		printf("B = 0\n"); 
		return(NULL); 
	} 
	else 
		return(CEILINGI(A, B)); 
} 
 
MPI *NEARINT(MPI *M, MPI *N){ 
/* Returns Z, the nearest integer to M/N, 
 * with Z=T if M/N=1/2+T. 
 */ 
	MPI *TEMP1, *TEMP2, *R; 
	 
	R = HALFMOD(M, N); 
	TEMP1 = SUBI(M, R); 
	TEMP2 = INT(TEMP1, N); 
	FREEMPI(R); 
	FREEMPI(TEMP1); 
	return(TEMP2); 
} 
 
MPI *NEARINTX(MPI *M, MPI *N){ 
/* Returns Z, the nearest integer to M/N, 
 * with Z=T if M/N=1/2+T. 
 */ 
	if(N->S <= 0){ 
		printf("N <= 0\n"); 
		return(NULL); 
	} 
	else 
		return(NEARINT(M, N)); 
} 
 
void POWERD(MPI *A, MPI *B, MPI *D, MPI *N, MPI **AA, MPI **BB){ 
/* *A + *B\sqrt(D)=(A + \sqrt(D))^N  */ 
    MPI *Y, *TEMP, *TEMP1, *TEMP2, *TEMP3, *TEMP4, *TEMP5; 
    MPI *X1, *X2, *T2; 
 
   X1 = COPYI(A); 
   X2 = COPYI(B); 
   Y = COPYI(N); 
   *AA = ONEI(); 
   *BB = ZEROI(); 
   while(Y->S){ 
      while((Y->V[0])%2 == 0){ 
          TEMP = Y; 
          Y = INT_(Y, 2); 
          FREEMPI(TEMP); 
          TEMP = X1; 
          TEMP1 = MULTI(X2, X2); 
          TEMP2 = MULTI(X1, X1); 
          TEMP3 = MULTI(D, TEMP1); 
          FREEMPI(TEMP1); 
          X1 = ADDI(TEMP2, TEMP3); 
          FREEMPI(TEMP2); 
          FREEMPI(TEMP3); 
          TEMP4 = MULTI(TEMP, X2); 
          FREEMPI(TEMP); 
          T2 = X2; 
          X2 = MULT_I(TEMP4, 2); 
          FREEMPI(TEMP4); 
          FREEMPI(T2); 
      } 
      TEMP = Y; 
      Y = SUB0_I(Y, 1); 
      FREEMPI(TEMP); 
      TEMP = *AA; 
      TEMP1 = MULTI(*BB, X2); 
      TEMP2 = MULTI(*AA, X1); 
      TEMP3 = MULTI(D, TEMP1); 
      FREEMPI(TEMP1); 
      *AA = ADDI(TEMP2, TEMP3); 
      FREEMPI(TEMP2); 
      FREEMPI(TEMP3); 
      TEMP4 = MULTI(TEMP, X2); 
      FREEMPI(TEMP); 
      TEMP5 = MULTI(*BB, X1); 
      T2 = *BB; 
      *BB = ADDI(TEMP4, TEMP5); 
      FREEMPI(T2); 
      FREEMPI(TEMP4); 
      FREEMPI(TEMP5); 
   } 
   FREEMPI(Y); 
   FREEMPI(X1); 
   FREEMPI(X2); 
   return; 
}