nfunc.c

calc大数库

				FREEMPI(tmp); FREEMPI(tmp1);
				tmp = MULTI(V, V); tmp1 = MULTI(M, M);
				*Yptr = ADD0I(tmp, tmp1);
				FREEMPI(tmp); FREEMPI(tmp1);
				FREEMPI(X);
				FREEMPI(L); FREEMPI(M);
				FREEMPI(U); FREEMPI(V);
				FREEMPI(R); FREEMPI(S);
				FREEMPI(C); FREEMPI(B);
				FREEMPI(E); FREEMPI(F);
				return (MINUS_ONEI());
			}
		}
		FREEMPI(R); FREEMPI(S); FREEMPI(E); FREEMPI(F);
	} 
}

MPI  *PEL(MPI *D, MPI *Eptr, MPI **Xptr, MPI **Yptr)
/*
 * This is a program for finding the least solution of Pell's equation
 * x*x - D*y*y = +-1.
 * The algorithm is based on K. Rosen, Elementary number theory
 * and its applications, p382, B.A. Venkov, Elementary Number theory, p.62 
 * and D. Knuth, Art of computer programming, Vol.2, p359, with Pohst's trick 
 * of using half the period.
 * The norm of the least solution is returned.
 * The partial quotients are printed iff *E != 0.
 */
{
	unsigned int i;
	MPI *X, *B, *C, *G, *Y, *F, *E, *L, *M, *N, *U, *V, *tmp;
	MPI *tmp1, *tmp2, *K, *R, *S;
	unsigned int t;
	FILE *outfile;
	char buff[20];

	X = BIG_MTHROOT(D, 2);
	strcpy(buff, "pell.out");
	outfile = fopen(buff, "w");
	if (Eptr->S)
	{
	printf("A[0] = "); PRINTI(X); printf("\n");
	fprintf(outfile, "A[0] = "); FPRINTI(outfile, X); fprintf(outfile,"\n");
	}
	B = COPYI(X);
	C = ONEI();
	G = MULTI(X, X); tmp = ADD0_I(G, (USL)1); FREEMPI(G);
	t = EQUALI(D, tmp); FREEMPI(tmp);
	if (t) /* period 1, exceptional case */
	{
		*Xptr = X;
		*Yptr = ONEI();
		printf("period length 1\n");
		fprintf(outfile, "period length 1\n");
		fclose(outfile);
		FREEMPI(B); FREEMPI(C);
		return (MINUS_ONEI());
	}
	L = ZEROI(); K = ONEI(); M = ONEI(); N = ZEROI();
	for (i = 0; 1; i++)
	{
		tmp = ADDI(X, B); Y = INT(tmp, C); FREEMPI(tmp);
		if (i && Eptr->S)
		{
			printf("A[%u] = ", i); PRINTI(Y); printf("\n");
			fprintf(outfile, "A[%u] = ", i); FPRINTI(outfile, Y); fprintf(outfile,"\n");
		}
		F = COPYI(B);
		tmp = B; tmp1 = MULTI(Y, C); B = SUBI(tmp1, B); 
		FREEMPI(tmp); FREEMPI(tmp1);
		E = COPYI(C); tmp = C;
		tmp1 = MULTI(B, B); tmp2 = SUBI(D, tmp1); C = INT(tmp2, C);
		FREEMPI(tmp); FREEMPI(tmp1); FREEMPI(tmp2);
		if (i == 0)
		{
			FREEMPI(Y); Y = COPYI(X);
		}
		R = L; S = M;
		tmp = MULTI(K, Y); U = ADDI(tmp, L); FREEMPI(tmp);
		tmp = MULTI(N, Y); V = ADDI(tmp, M); FREEMPI(tmp);
		FREEMPI(Y);
		L = K; K = U;
		M = N; N = V;
/* U/V is the ith convergent to sqrt(D) */
		if (i)
		{
			if (EQUALI(B, F))
			/*\alpha_H=\alpha_{H+1}, even period 2H */
			{
				tmp = ADDI(U, R); tmp1 = MULTI(M, tmp);
				FREEMPI(tmp);
				if (i % 2 == 0)
					*Xptr = ADD0_I(tmp1, (USL)1);
				else
					*Xptr = SUB0_I(tmp1, (USL)1);
				FREEMPI(tmp1); 
				tmp = ADDI(V, S); *Yptr = MULTI(M, tmp);
				FREEMPI(tmp);
				FREEMPI(X);
				FREEMPI(L); FREEMPI(M);
				FREEMPI(U); FREEMPI(V);
				FREEMPI(R); FREEMPI(S);
				FREEMPI(C); FREEMPI(B);
				FREEMPI(E); FREEMPI(F);
				printf("even period length %u\n", 2*i);
				fprintf(outfile, "even period length %u\n", 2*i);
				fclose(outfile);
				return (ONEI());
			}
			if (EQUALI(C, E)) 
			/*\beta_H=\beta_{H-1}, odd period 2H-1 */
			{
				tmp = MULTI(U, V); tmp1 = MULTI(L, M);
				*Xptr = ADDI(tmp, tmp1);
				FREEMPI(tmp); FREEMPI(tmp1);
				tmp = MULTI(V, V); tmp1 = MULTI(M, M);
				*Yptr = ADD0I(tmp, tmp1);
				FREEMPI(tmp); FREEMPI(tmp1);
				FREEMPI(X);
				FREEMPI(L); FREEMPI(M);
				FREEMPI(U); FREEMPI(V);
				FREEMPI(R); FREEMPI(S);
				FREEMPI(C); FREEMPI(B);
				FREEMPI(E); FREEMPI(F);
				printf("odd period length %u\n", 2*i+1);
				fprintf(outfile, "odd period length %u\n", 2*i+1);
				fclose(outfile);
				return (MINUS_ONEI());
			}
		}
		FREEMPI(R); FREEMPI(S); FREEMPI(E); FREEMPI(F);
	} 
}

MPIA A_SURD; /* used in REDUCED() */
MPIA U_SURD; /* used in REDUCED() */
MPIA V_SURD; /* used in REDUCED() */
unsigned int REDUCED(MPI *D, MPI *U, MPI *V, USI i)
/*
 * This is a function for finding the period of the continued fraction
 * expansion of reduced quadratic irrational a=(U+sqrt(D))/V.
 * Here D is non-square, 1<(U+sqrt(D))/V, -1<(U-sqrt(D))/V<0.
 * The algorithm also assumes that V divides d-U*U and is based on K. Rosen,
 * Elementary Number theory and its applications, p.379-381 and Knuth's The art
 * of computer programming, Vol. 2, p. 359. The period length is returned if
 * a is reduced.
 * variable i is created by SURD(D,T,U,V,P_ARRAY,Q_ARRAY) below and indexes 
 * the ith convergent of (U+T*sqrt(D))/V. 
*/
{
	MPI *A, *F, *R, *S, *tmp, *tmp1, *tmp2;
	unsigned int j;
	FILE *outfile;
	char buff[20];

	F = BIG_MTHROOT(D, 2);
	R = COPYI(U); S = COPYI(V);
	strcpy(buff, "surd.out");
	outfile = fopen(buff, "a");
	for(j = i; 1; j++)
	{
		tmp = ADD0I(F, R); A = INT0(tmp, S); FREEMPI(tmp);
		ADD_TO_MPIA(A_SURD, A, j);
		fprintf(outfile,", A[%u]=", j); FPRINTI(outfile, A);
		fprintf(outfile,", PERIOD\n");
		tmp = MULTI(A, S); tmp1 = R; R = SUB0I(tmp, R);
		 FREEMPI(tmp); FREEMPI(tmp1);
		tmp = S; tmp1 = MULTI(R, R);
		tmp2 = SUB0I(D, tmp1); S = INT0(tmp2, S);
		ADD_TO_MPIA(U_SURD, R, j + 1);
		fprintf(outfile,"P[%u]=", j + 1);
		FPRINTI(outfile, R);
		ADD_TO_MPIA(V_SURD, S, j + 1);
		fprintf(outfile,", Q[%u]=", j + 1);
		FPRINTI(outfile, S);
		FREEMPI(tmp); FREEMPI(tmp1); FREEMPI(tmp2);
		FREEMPI(A);
		if (EQUALI(U, R) && EQUALI(V, S))
		{
			fprintf(outfile,"\n");
			FREEMPI(R); FREEMPI(S); FREEMPI(F);
			fclose(outfile);
			return (j + 1 - i);
		}
		if(j == R0){
			FREEMPI(R); FREEMPI(S); FREEMPI(F);
			execerror("j = R0", "");
		}
	}
}

unsigned int SURD(MPI *D, MPI *T, MPI *U, MPI *V, MPIA *AA_SURD, MPIA *UU_SURD, MPIA *VV_SURD, MPIA *P_SURD, MPIA *Q_SURD, USI surd_flag)
/*
 * This function uses the continued fraction algorithm expansion in K. Rosen,
 * Elementary Number theory and its applications,p.379-381 and Knuth's
 * The art of computer programming, Vol.2, p. 359. It locates the first complete
 * quotient that is reduced and then uses the function REDUCED(D,U,V,i) above
 * to locate and return the period of the continued fraction expansion
 * of (U+T*sqrt(D))/V. 
 * AA_SURD is the sequence of partial quotients up to the end of the period;
 * UU_SURD and VV_SDURD are the sequences of U[i] and V[i], where the i-th
 * complete quotient is (U[i]+sqrt(D))/V[i];
 * P_SURD/Q_SURD give the convergents.
 * Output is sent to surd.out.
 * Regarding surd_flag (added 29th May 2000). This is needed in PATZ.
 * If surd_flag = 1 and the period length is odd, we do an extra period.
 */
{
	MPI *A, *DD, *F, *W, *tmp, *tmp1, *tmp2, *UU, *VV, *X;
	unsigned int i, j, l;
	int z, t;
	FILE *outfile;
	char buff[20];

	A_SURD = BUILDMPIA();
	U_SURD = BUILDMPIA();
	V_SURD = BUILDMPIA();
	F = BIG_MTHROOT(D, 2);
	UU = COPYI(U); VV = COPYI(V);
	z = T->S;
	UU->S = (U->S) * z;
	VV->S = (V->S) * z;
	DD = MULTI3(D, T, T);
	W = ABSI(VV);
	tmp1 = MULTI(UU, UU); tmp2 = SUBI(DD, tmp1); FREEMPI(tmp1);
	tmp1 = MOD(tmp2, W); FREEMPI(tmp2);
	if (tmp1->S)
	{
		tmp2 = DD; DD = MULTI3(DD, VV, VV); FREEMPI(tmp2);
		tmp2 = UU; UU = MULTI(UU, W); FREEMPI(tmp2);
		tmp2 = VV; VV = MULTI(VV, W); FREEMPI(tmp2);
	}
	FREEMPI(W);
	FREEMPI(tmp1);
	FREEMPI(F); 
	F = BIG_MTHROOT(DD, 2);
	strcpy(buff, "surd.out");
	outfile = fopen(buff, "w");
	if (V->S == -1)
		fprintf(outfile, "-");
	if (!EQONEI(V))
		fprintf(outfile, "(");
	if(U->S)
		FPRINTI(outfile, U); 
	if(T->S >0 && U->S)
		fprintf(outfile," + ");
	if(!EQONEI(T) && !EQMINUSONEI(T))
	{
		FPRINTI(outfile, T); 
		fprintf(outfile,"*");
	}
	if (EQMINUSONEI(T))
        	fprintf(outfile,"-");
        fprintf(outfile,"sqrt(");
	FPRINTI(outfile, D); fprintf(outfile,")");
	if(!EQONEI(V))
		fprintf(outfile,")");
	if (!EQONEI(V) && !EQMINUSONEI(V))
	{
		fprintf(outfile,"/");
		X = ABSI(V);
		FPRINTI(outfile, X); 
		FREEMPI(X);
	}
	fprintf(outfile,":\n");
	for (j = 0; 1; j++)
	{
		if(j == R0){
			FREEMPIA(A_SURD);
			FREEMPIA(U_SURD);
			FREEMPIA(V_SURD);
			FREEMPI(DD);
			FREEMPI(F);
			FREEMPI(UU);
			FREEMPI(VV);
			fclose(outfile);
			execerror("j = R0", "");
		}
		ADD_TO_MPIA(U_SURD, UU, j);
		fprintf(outfile,"P[%u]=", j);
		FPRINTI(outfile, UU);
		ADD_TO_MPIA(V_SURD, VV, j);
		fprintf(outfile,", Q[%u]=", j);
		FPRINTI(outfile, VV);
		if (VV->S > 0 && UU->S > 0 && (RSV(UU, F) <= 0))
		{
			tmp1 = ADD0I(UU, VV);
			z = RSV(tmp1, F);
			FREEMPI(tmp1);
			tmp1 = SUBI(VV, UU);
			t = RSV(tmp1, F);
			FREEMPI(tmp1);
			if (z == 1 && t <= 0)/* (U+sqrt(D))/V is reduced */
			{
				fclose(outfile);
				l = REDUCED(DD, UU, VV, j);
				if(surd_flag && l % 2)
					l = REDUCED(DD, UU, VV, j + l);
				outfile = fopen(buff, "a");
				fprintf(outfile,"period length =  %u\n", l);
				CONVERGENTS(A_SURD, P_SURD, Q_SURD);
				FREEMPI(DD);
				FREEMPI(F);
				FREEMPI(UU);
				FREEMPI(VV);
				fprintf(outfile,"convergents:\n");
				for(i = 0; i < A_SURD->size; i++)
				{
					fprintf(outfile,"A[%u]/B[%u]=", i, i);
					FPRINTI(outfile, (*(P_SURD))->A[i]);
					fprintf(outfile,"/");
					FPRINTI(outfile, (*(Q_SURD))->A[i]);
					fprintf(outfile,"\n");
				}
/* Actually U_SURD and V_SURD are one unit longer than A_SURD */
				
				*(AA_SURD) = BUILDMPIA();
				for(i = 0; i < A_SURD->size; i++)
					ADD_TO_MPIA(*(AA_SURD), A_SURD->A[i], i);
				*(UU_SURD) = BUILDMPIA();
				*(VV_SURD) = BUILDMPIA();
				for(i = 0; i < A_SURD->size; i++)
				{
					ADD_TO_MPIA(*(UU_SURD), U_SURD->A[i], i);
					ADD_TO_MPIA(*(VV_SURD), V_SURD->A[i], i);
				}
				FREEMPIA(A_SURD);
				FREEMPIA(U_SURD);
				FREEMPIA(V_SURD);
				fclose(outfile);
				return (l); /* l is the period length */
			}
		}
		tmp1 = ADDI(F, UU);
		if (VV->S > 0)
			 A = INT(tmp1, VV);
		else /* See Knuth p. 359 */
		{
			tmp = ONEI();
			tmp2 = ADDI(tmp1, tmp); A = INTI(tmp2, VV);
			FREEMPI(tmp); FREEMPI(tmp2);
		}
		FREEMPI(tmp1); 
		fprintf(outfile,", A[%u]=", j);
		FPRINTI(outfile, A);
		fprintf(outfile,"\n");
		ADD_TO_MPIA(A_SURD, A, j);
		tmp2 = MULTI(A, VV); tmp1 = UU; UU = SUBI(tmp2, UU);
	 	FREEMPI(A); FREEMPI(tmp2); FREEMPI(tmp1);
		tmp1 = MULTI(UU, UU); tmp2 = SUBI(DD, tmp1); FREEMPI(tmp1);
		tmp1 = VV; VV = INTI(tmp2, VV); FREEMPI(tmp1); FREEMPI(tmp2);
	}
}
 
MPI * JUGGLERT(MPI *Dptr)
{
	MPI *N, *M;
	unsigned long int t;

	t = MOD0_(Dptr, 3);
	printf("res class mod %lu\n", t);
	if (t == 0)
		N = COPYI(Dptr);
	else if (t == 1)
		N = POWERI(Dptr, 2);
	else
		N = POWERI(Dptr, 13);
	M = BIG_MTHROOT(N, 3);
	FREEMPI(N);
	return (M);
}

void JUGGLER(MPI *Dptr, MPI *Iptr)
{
	unsigned long i, k;
	MPI *X, *Y, *Temp;

	Y = BUILDMPI(1);
	Y->S = 1;
	Y->V[0] = 2;
	k = CONVERTI(Iptr);
	X = COPYI(Dptr);
	for(i = 0; i <= k; i++)
	{
		printf("i = %lu: size = %u\n", i, X->D);
		if (EQUALI(X, Y))
		{
		  printf("starting number = ");PRINTI(Dptr);printf("\n");
		  printf("the number of iterations taken to reach 2 is %lu\n", i);
			break;
		}
		else if (EQONEI(X))
		{
		  printf("starting number = ");PRINTI(Dptr);printf("\n");
		  printf("the number of iterations taken to reach 1 is %lu\n", i);
			break;
		}
		Temp = X;
		/*printf("i = %lu, LENGTHI(X[%lu])=%lu\n", i, i, LENGTHI(X));*/
		X = JUGGLERT(Temp);
		FREEMPI(Temp);
	/*	PRINTI(X);printf("\n"); */
	}
	FREEMPI(X);
	FREEMPI(Y);
	return;
}

void HERMITE()
{
	USI *alpha, nz, answer;
	MPMATI *MATI1, *MATI2, *MATI3, *MATI4;
	FILE *outfile;
	char buff[20];

	printf("Do you wish to use an existing matrix from a file? (Y/N)\n");
	answer = GetYN();
	if (answer)
	{
		MATI1 = FINPUTMATFILEI_I();
		if (MATI1 == NULL)
			exit (0);
	}
	else
	{
		printf("enter the matrix of integers (first row non-zero) :\n");
		MATI1 = INPUTMATI();
	}
	printf("The matrix entered is\n\n");
	PRINTMATI(0,MATI1->R-1,0,MATI1->C-1,MATI1);
	printf("Do you want the unimodular transformation matrix P? (Y/N)\n");
	answer = GetYN();
	if (!answer)
	{
		alpha = KB_ROW(MATI1, &nz);
		MATI2 = HERMITE1(MATI1, alpha, nz);
		printf("The Hermite normal form = \n");
		PRINTMATI(0,MATI2->R-1,0,MATI2->C-1,MATI2);
		strcpy(buff, "hermite.out");
		outfile = fopen(buff, "w");
		FPRINTMATI(outfile,0,MATI2->R-1,0,MATI2->C-1,MATI2);
		fclose(outfile);
	}
	else
	{
		alpha = KB_ROWP(MATI1, &MATI3, &nz);
		MATI2 = HERMITE1P(MATI1, MATI3, &MATI4, alpha, nz);
		FREEMATI(MATI3);
		printf("The Hermite normal form = \n");
		PRINTMATI(0,MATI2->R-1,0,MATI2->C-1,MATI2);
		strcpy(buff, "hermite.out");
		outfile = fopen(buff, "w");
		FPRINTMATI(outfile,0,MATI2->R-1,0,MATI2->C-1,MATI2);
		fclose(outfile);
		printf("The unimodular transformation matrix P = \n");
		PRINTMATI(0,MATI4->R-1,0,MATI4->C-1,MATI4);
		strcpy(buff, "hermitep.out");
		outfile = fopen(buff, "w");
		FPRINTMATI(outfile,0,MATI4->R-1,0,MATI4->C-1,MATI4);
		fprintf(outfile, "%u", nz);
		fclose(outfile);
		FREEMATI(MATI4);
	}
	FREEMATI(MATI2);
	FREEMATI(MATI1);
	ffree((char *)alpha, (MATI1->C) * sizeof(USI));
	return;
}

void MLLL()
{
	USI answer, m1, n1;
	MPMATI *MATI1, *MATI2, *MATI3;
	char buff[20];
	FILE *outfile;

	printf("Do you wish to use an existing matrix from a file? (Y/N)\n");
	answer = GetYN();
	if (answer)
	{
		MATI1 = FINPUTMATFILEI_I();
		if (MATI1 == NULL)
			exit (0);
	}
	else
	{
		printf("enter the matrix of integers (first row non-zero) :\n");
		MATI1 = INPUTMATI();
	}
	printf("The matrix entered is\n\n");
	PRINTMATI(0,MATI1->R-1,0,MATI1->C-1,MATI1);
	MLLLVERBOSE = 0;
	printf("VERBOSE? (Y/N)\n");
	answer = GetYN();
	if (answer)
		MLLLVERBOSE = 1;
	MATI3 = IDENTITYI(MATI1->R);
	printf("enter the parameters m1 and n1 (normally 1 and 1) :");
	scanf("%u %u", &m1, &n1);
	Flush();
	MATI2 = BASIS_REDUCTION(MATI1, &MATI3, 0, m1, n1);
	printf("\n\n");
	printf("The corresponding transformation matrix is\n");
	PRINTMATI(0,MATI3->R-1,0,MATI3->C-1,MATI3);
	strcpy(buff, "mllltran.out");
	outfile = fopen(buff, "w");
	FPRINTMATI(outfile,0,MATI3->R-1,0,MATI3->C-1,MATI3);
	fclose(outfile);
	printf("The corresponding reduced basis is\n");
	PRINTMATI(0,MATI2->R-1,0,MATI2->C-1,MATI2);
	strcpy(buff, "mlllbas.out");
	outfile = fopen(buff, "w");
	FPRINTMATI(outfile,0,MATI2->R-1,0,MATI2->C-1,MATI2);
	fclose(outfile);
	FREEMATI(MATI1);
	FREEMATI(MATI2);
	FREEMATI(MATI3);
	return;
}

void SMITH()
{
	MPI **M, *MAX_ENTRY;
	unsigned int i, j, u, r, z, answer, m1, n1;
	MPMATI *MATI1, *MATI2, *MATI3, *Temp, *Temp1;
	char buff[20];
	FILE *outfile;

	printf("This program takes as input a rectangular matrix A of integers and computes\n");
	printf("unimodular integer matrices P, Q such that PAQ=D, where D is a diagonal matrix\n");
	printf("with positive integer diagonal elements d[1],...,d[r].\n");
 	printf("Here d[i] divides d[i+1], 1<=i<=r-1.\n");
	printf("d[1],...,d[r] are called the invariant factors of A.\n\n");
	printf("Do you wish to use an existing matrix from a file? (Y/N)\n\n");
	answer = GetYN();
	if (answer)
	{
		MATI1 = FINPUTMATFILEI_I();
		if (MATI1 == NULL)
			exit (0);
	}
	else
	{
		printf("enter the matrix of integers:\n");
		MATI1 = INPUTMATI();
	}