i5i.c

calc大数库

MPI *GCD(MPI *Aptr, MPI *Bptr)
/*
 * GCD(|*Aptr|, |*Bptr|). 
 */
{
	MPI *A, *B, *Cptr, *Temp;

	A = ABSI(Aptr);
	if (Bptr->S == 0)
		return A;
	B = ABSI(Bptr);
	Cptr = MOD0(A, B);
	while (Cptr->S == 1)
	{
		Temp = A;
		A = B;
		FREEMPI(Temp);
		B = Cptr;
		Cptr = MOD0(A, B);
	}
	Temp = Cptr;
	Cptr = B;
	FREEMPI(Temp);
	FREEMPI(A);
	return Cptr;
}

MPI *LCM(MPI *Aptr, MPI *Bptr)
/*
 * *Cptr = LCM(*Aptr, *Bptr). 
 */
{
	MPI *C, *N, *Cptr, *Temp;

	if (Aptr->S == 0 || Bptr->S == 0)
		return ZEROI();
	else
	{
		N = GCD(Aptr, Bptr);
		C = MULTI(Aptr, Bptr);
		Temp = C;
		C = ABSI(Temp);
		FREEMPI(Temp);
		Cptr = INT0(C, N);
		FREEMPI(C);
		FREEMPI(N);
		return Cptr;
	}
}

MPI *LCM_ARRAY(MPIA M)
/*
 * *Aptr = LCM(M[0], ..., M[n - 1]).
 */
{
	MPI *L, *Aptr;
	unsigned int i;
	USL n = M->size;
	if (n == 0) 
	  return ZEROI();
	
	Aptr = COPYI(M->A[0]);
	if (n == 1)
		return Aptr;
	else
	{
		for ( i = 1; i < n; i++)
		{
			L = Aptr;
			Aptr = LCM(L, M->A[i]);
			FREEMPI(L);
		}
		return Aptr;
	}
}

void CONT_FRAC(MPI *Aptr, MPI *Bptr, MPI *Zptr, MPI **Xptr, MPI **Yptr)
/*
 * Euclid's algorithm for use in 2SQUARES. 
 */
{
	MPI *A, *B, *C, *Temp;

	A = COPYI(Aptr);
	B = COPYI(Bptr);
	C = MOD0(A, B);
	while (C->S)
	{
		if (RSV(C, Zptr) <= 0)
		{
			*Xptr = COPYI(C);
			Temp = A;
			A = B;
			FREEMPI(Temp);
			B = C;
			*Yptr = MOD0(A, B);
			FREEMPI(A);
			FREEMPI(B);
			return;
		}
		Temp = A;
		A = B;
		FREEMPI(Temp);
		B = C;
		C = MOD0(A, B);
	}
}

void INTSETUI(unsigned int *Aptr, unsigned long n, unsigned int m)
/*
 * this converts the n array elements *Aptr,....*(Aptr + n - 1) to m. 
 */
{
	for( ; n; n--)
		*(Aptr++) = m;
	return;
}

void INTSETUL(unsigned long *Aptr, unsigned long n, unsigned long m)
/*
 * this converts the n array elements *Aptr,....*(Aptr + n - 1) to m. 
 */
{
	for( ; n; n--)
		*(Aptr++) = m;
	return;
}

void INTSETI(int *Aptr, unsigned long n, int m)
/*
 * this converts the n array elements *Aptr,....*(Aptr + n - 1) to m. 
 */
{
	for( ; n; n--)
		*(Aptr++) = m;
	return;
}

void INTSETL(long *Aptr, unsigned long n, long m)
/*
 * this converts the n array elements *Aptr,....*(Aptr + n - 1) to m. 
 */
{
	for( ; n; n--)
		*(Aptr++) = m;
	return;
}

unsigned long LENGTHI(MPI *Mptr)
/*
 * returns the number of decimal digits in the MPI *Mptr,
 * increased by 1 if *Mptr is negative.
 */
{
	MPI *F, *Temp;
	unsigned long i = 0, s;
	
	if (Mptr->S >= 0)
		s = 0;
	else
		s = 1;
	F = ABSI(Mptr);
	do
	{
		Temp = F;
		F = INT0_(Temp, (USL)10);
		FREEMPI(Temp);
		i++;
	} while (F->S == 1);
	FREEMPI(F);
	return (i + s);
}

unsigned long CONVERTI(MPI *N)
/*
 * Returns |N| as an unsigned long, providing |N| < 2^32.
 * If N's too big, you get an overflow error.
 */
{
	if (N->D > 1)
	{
		fprintf(stderr, "overflow in CONVERTI\n");
		exit(1);
	}
	if (N->S == 0)
		return (USL)0;
	if (N->D == 0)
		return N->V[0];
	else
		return ((N->V[1]) * R0 + (N->V[0]));
}

MPI *INPUTSI(char **ptr, unsigned int *uptr)
/*
 * Converts decimal input from a stream into an MPI, leaving *ptr to point
 * to the character following the last decimal digit.
 * Ignores the combination '\' followed by '\n'.
 * If a rubbish character is met before decimal input, Mptr is set to 0
 * and *uptr = 0 is returned.
 * If a rubbish character is met immediately after decimal input, 0 is
 * returned,  Otherwise *uptr = 1 is returned.
 * In any case *Mptr is set equal to any inputted decimal.
 */
{
	unsigned long n;
	int t;
	MPI *Mptr, *Temp;

	Mptr = ZEROI();
	while (**ptr == ' ' || **ptr == '\n')
		(*ptr)++;
	if (**ptr != '-' && (**ptr < '0' || **ptr > '9')) 
	{
		printf("in INPUTSI1, illegal character %c entered:\n", **ptr);
		*uptr = 0;
		return (Mptr);
	}
	else if (**ptr == '0')
	{
		(*ptr)++;
		if (**ptr >=0 || **ptr <= 9)
			t = 1;
		else if (**ptr != ' ' && **ptr != '\n' && **ptr != 'i' && **ptr != '/' && **ptr != '+' && **ptr != '-' && **ptr != '\0')
		{
			printf("in INPUTSI2, illegal character %c entered:\n", **ptr);
			*uptr = 0;
			return (Mptr);
		}
		else
		{
			*uptr = 1;
			return (Mptr);
		}
	}
	else if (**ptr == '-') 
	{
		do
			(*ptr)++;
		while (**ptr == ' ');
		if (**ptr <= '0' || **ptr > '9')
		{
			printf("in INPUTSI3 illegal character %c entered:\n", **ptr);
			*uptr = 0;
			return (Mptr);
		}
		t = -1;
		while (**ptr == ' ')
			(*ptr)++;
	}
        else
		t = 1;
	while ((**ptr == '\\') || (**ptr >= '0' && **ptr <= '9')) 
	{
		if (**ptr == '\\')
		/* for use with bc output, where '\' is followed by '\n'*/
			(*ptr)++;
		else 
		{
			Temp = Mptr;
			Mptr = MULT_I(Temp, 10L);
			FREEMPI(Temp);
			n = (unsigned long)(**ptr - '0');
			Temp = Mptr;
			Mptr = ADD0_I(Temp, n);
			FREEMPI(Temp);
		}
		(*ptr)++;
	}
	if (t == -1)
		Mptr->S = -1;
	if (**ptr != '.' && **ptr != ' ' && **ptr != '\n' && **ptr != 'i' && **ptr != '/' && **ptr != '+' && **ptr != '-' && **ptr != '\0')
	{
		printf("in INPUTSI4 illegal character %c entered:\n", **ptr);
		*uptr = 0;
		return (Mptr);
	}
	*uptr = 1;
	return (Mptr);
}

void SPRINTI(char *buffer, MPI *Mptr)
/*
 * The decimal digits of the MPI *Mptr are placed in the string buffer.
 * no new-line is incorporated. 
 */
{
	MPI *F, *Temp;
	unsigned long *A, l ;
	int i = 0;
	
	if (Mptr->S == 0) 
	{
		sprintf(buffer, "%lu", (USL)0);
		return;
	}
	if (Mptr->S == -1) 
	{
		sprintf(buffer, "-");
		buffer++;
	}
	F = ABSI(Mptr);
	l = LENGTHI(F);
	A = (unsigned long *)mmalloc((USL)(l * sizeof(unsigned long)));
	do
        {
		A[i] = MOD0_(F, (USL)10);
		Temp = F;
		F = INT0_(Temp, (USL)10);
		FREEMPI(Temp);
		i++;
	}
 	while (F->S);
	i--;
	while (i >= 0)
        {
		sprintf(buffer, "%lu", A[i]);
		buffer++;
		i--;
	}
	ffree((char *)A, l * sizeof(unsigned long)); 
	FREEMPI(F);
	return;
}

MPI *NEAREST_INTI(MPI *Aptr, MPI *Bptr)
/*
 * Y is the nearest integer to *Aptr/(*Bptr).
 */
{
	MPI *X, *Y, *Z, *Tmp1, *Tmp2;

	Y = INTI(Aptr, Bptr);
	X = MULTI(Bptr, Y);
	Z = SUBI(Aptr, X);
	FREEMPI(X);
	Tmp1 = Z;
	Z = MULT_I(Z, 2L);
	FREEMPI(Tmp1);
	if (RSV(Z, Bptr) == 1)
	{
		Tmp2 = ONEI();
		Tmp1 = Y;
		Y = ADDI(Y, Tmp2);
		FREEMPI(Tmp2);
		FREEMPI(Tmp1);
	}
	FREEMPI(Z);


/*
	printf("Aptr = ");PRINTI(Aptr);printf("\n");
	printf("Bptr = ");PRINTI(Bptr);printf("\n");
	printf("Y = ");PRINTI(Y);printf("\n");

*/

	return (Y);
}

MPI *GCD_ARRAY(MPIA M)
/*
 * *Aptr = GCD(M[0], ..., M[n - 1]).
 */
{
	MPI *L, *Aptr;
	unsigned int i;
	unsigned n = M->size;
	Aptr = COPYI(M->A[0]);
	if (n == 1)
		return Aptr;
	else
	{
		for (i = 1; i < n; i++)
		{
			L = Aptr;
			Aptr = GCD(L, M->A[i]);
			FREEMPI(L);
		}
		return Aptr;
	}
}

void PRINTIA(MPIA Mptr)
/* 
 * Prints an array Mptr[0],...,Mptr[n], n=Aptr, consisting of *Aptr MPI's.
 */
{

	USI i;
	for (i = 0; i < Mptr->size; i++)
	{
		printf("[%u]:" , i);
		PRINTI(Mptr->A[i]);
		printf("\n");
	}
	return;
}

/*void EUCLID(MPI *Aptr, MPI *Bptr, MPI **C[], MPI **Dptr)*/
 /* Returns the partial quotients C[0],...C[n] of the continued fraction 
 * expansion of Aptr/Bptr. Here *Dptr = n.
 */
/*{
	MPI *A, *B, *Cptr, *Temp;
	USL n, i;

	n = 2 * BINARYB(Bptr) + 1;
	*C = (MPI **)mmalloc((USL)(n * sizeof(MPI *)));
	A = COPYI(Aptr);
	B = COPYI(Bptr);
	(*C)[0] = INT0(A, B);
	Cptr = MOD0(A, B);
	i = 1;
	while (Cptr->S == 1)
	{
		Temp = A;
		A = B;
		FREEMPI(Temp);
		B = Cptr;
		(*C)[i] = INT0(A, B);
		Cptr = MOD0(A, B);
		i++;
	}
	*C = (MPI **)rrealloc(*C, (USL)(i * sizeof(MPI *)), -((long)(n - i) * sizeof(MPI *)));
	*Dptr = CHANGE(i - 1);
	FREEMPI(A);
	FREEMPI(B);
	FREEMPI(Cptr);
	return;
}
*/

void CONVERGENTS(MPIA A, MPIA *P, MPIA *Q)
/*
 * P(0)=A[0],P[1)=A[0]*A[1]+1,P[i+1]=A[i+1]*P[i]+P[i-1] if i >= 1.
 * Q(0)=1,Q[1)=A[1],Q[i+1]=A[i+1]*Q[i]+Q[i-1] if i >= 1.
 * The arrays P[0],...,P[N] and Q[0],..., Q[N] are returned.
 */
{
        USL i, n;
	MPI *X1, *X2, *Y1, *Y2, *Z1, *Z2, *temp=NULL;
	MPI *ONE;
	char buff[20];
	FILE *outfile;
        strcpy(buff, "convergents.out");
	if(access(buff, R_OK) == 0)
	  unlink(buff);
	outfile = fopen(buff, "w");

	ONE = ONEI();
	n = A->size-1;
	Z1=(MPI *)NULL;
	Z2=(MPI *)NULL;

	*P = BUILDMPIA();
	*Q = BUILDMPIA();
	ADD_TO_MPIA(*P, A->A[0], 0);
	ADD_TO_MPIA(*Q, ONE, 0);

	if (n == 0)
	{
		fprintf(outfile, "P[0]=");
                FPRINTI(outfile, (*P)->A[0]);
		fprintf(outfile, ";Q[0]=");
                FPRINTI(outfile, (*Q)->A[0]);
		fprintf(outfile, "\n");
		fclose(outfile);
		FREEMPI(ONE);
		return;
	}
	temp = MULTABC(ONE, A->A[0], A->A[1]);
	ADD_TO_MPIA(*P, temp, 1);
	FREEMPI(temp);
	ADD_TO_MPIA(*Q, A->A[1], 1);
	if (n == 1) 
	{
		for (i = 0; i <= n;i++)
		{
			fprintf(outfile, "P[%lu]=", i);
			FPRINTI(outfile, (*P)->A[i]);
			fprintf(outfile, ";Q[%lu]=", i);
			FPRINTI(outfile, (*Q)->A[i]);
			fprintf(outfile, "\n");
		}
		fclose(outfile);
		FREEMPI(ONE);
		return;
	}
	X1 = COPYI(A->A[0]);
	Y1 = MULTABC(ONE, A->A[0], A->A[1]);
	X2 = ONEI();
	Y2 = COPYI(A->A[1]);
	for (i = 2; i <= n; i++)
	{
		Z1 = MULTABC(X1, A->A[i], Y1);
		Z2 = MULTABC(X2, A->A[i], Y2);
		ADD_TO_MPIA(*P, Z1, i);
		ADD_TO_MPIA(*Q, Z2, i);
		FREEMPI(X1);
		FREEMPI(X2);
		X1 = Y1;
		Y1 = Z1;
		X2 = Y2;
		Y2 = Z2;
	}
	FREEMPI(X1);
	FREEMPI(X2);
	FREEMPI(Z1);
	FREEMPI(Z2);
	for (i = 0; i <= n;i++)
	{