nfunc.c

calc大数库

/* nfunc.c */
#ifdef _WIN32
#include "unistd_DOS.h"
#else
#include <unistd.h>
#endif
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "integer.h"
#include "fun.h"

extern unsigned long PRIME[];
unsigned int MLLLVERBOSE;
unsigned int HERMITEVERBOSE;
unsigned int HERMITE1VERBOSE;
unsigned int GCDVERBOSE;
unsigned int HAVASFLAG = 0;
unsigned int FPRINTMATIFLAG = 0;
unsigned int GCD_MAX = 100000;
extern unsigned int GCDFLAG;

MPI *EUCLIDI(MPI *Pptr, MPI *Qptr, MPI **Hptr, MPI **Kptr)
/*
 * gcd(Pptr, Qptr) = Hptr * Pptr + Kptr * Qptr.
 */
{
	MPI *Q, *A, *B, *C, *H1, *H2, *K1, *K2, *L1, *L2, *Tmp1, *Tmp2;
	int s;

	if (Qptr->S == 0)
	{
		s = Pptr->S;
		if (Pptr->S != 0)
		{
			if (s == 1)
				*Hptr = ONEI();
			else
				*Hptr = MINUS_ONEI();
			*Kptr = ZEROI();
			return (ABSI(Pptr));
		}
		else
		{
			*Hptr = ZEROI();
			*Kptr = ZEROI();
			return (ZEROI());
		}
	}
	A = COPYI(Pptr);
	B = ABSI(Qptr);
	C = MOD(A, B);
	s = Qptr->S;
	if (C->S == 0)
	{
		if (s == 1)
			*Kptr = ONEI();
		else
			*Kptr = MINUS_ONEI();
		*Hptr = ZEROI();
		FREEMPI(A);
		FREEMPI(C);
		return (B);
	}
	L1 = ONEI();
	K1 = ZEROI();
	L2 = ZEROI();
	K2 = ONEI();
	while (C->S != 0)
	{
		Q = INT(A, B);
		FREEMPI(A);
		A = B;
		B = C;
		C = MOD0(A, B);
		Tmp1 = MULTI(Q, K1);
		Tmp2 = MULTI(Q, K2);
		FREEMPI(Q);
		H1 = SUBI(L1, Tmp1);
		H2 = SUBI(L2, Tmp2);
		FREEMPI(Tmp1);
		FREEMPI(Tmp2);
		FREEMPI(L1);
		L1 = K1;
		FREEMPI(L2);
		L2 = K2;
		K1 = H1;
		K2 = H2;
	}
	*Hptr = K1;
	if (s == -1)
		K2->S = -(K2->S);		
	*Kptr = K2;
	FREEMPI(L1);
	FREEMPI(L2);
	FREEMPI(A);
	FREEMPI(C);
	return (B);
}

void SERRET(MPI *P, MPI **Xptr, MPI **Yptr)
/*
 * This program finds positive integers X, Y such that "X*X+Y*Y=P, where P is a
 * prime, p=4n+1. The algorithm goes back to Serret and is from the book 
 * "Computational methods in number theory, part 1," edited by H.W.Lenstra and
 * R.Tijdeman.
 */
{
	int i, j;
	MPI *Q, *N, *Z, *U, *V, *W, *Tmp;

	j = 0;
	Q = SUB0_I(P, (USL)1);
	Z = BIG_MTHROOT(Q, 2);
	W = MULTI(Z, Z);
	if (EQUALI(W, Q))
	{
		PRINTI(P); printf(" = "); PRINTI(Z); printf("\n");
		*Xptr = Z;
		*Yptr = ONEI();
		FREEMPI(Q); FREEMPI(W);
		return;
	}
	N = INT_(Q, (USL)4);
	FREEMPI(W); FREEMPI(Q);
	for (i = 0; i <= Y0 - 1; i++)
	{
		printf("PRIME[%d] = %lu\n", i, PRIME[i]);
		U = MPOWER_((long)PRIME[i], N, P);
		V = MULTI(U, U);
		Tmp = V; V = MOD0(V, P); FREEMPI(Tmp);
		Tmp = V; V = ADD0_I(V, (USL)1); FREEMPI(Tmp);
		if (EQUALI(V, P))
		{
			j = 1;
			printf("U = "); PRINTI(U); printf("\n");
			/* U is a square-root of -1 mod P */
			FREEMPI(V); FREEMPI(N);
			break;
		}
		FREEMPI(U); FREEMPI(V);
	}
	if (j == 1)
	{
		CONT_FRAC(P, U, Z, Xptr, Yptr);
		printf("P = X^2 + Y^2, where\n");
		printf("P = "); PRINTI(P); printf("\n");
		printf("X = "); PRINTI(*Xptr); printf("\n");
		printf("Y = "); PRINTI(*Yptr); printf("\n");
		FREEMPI(U); FREEMPI(Z);
		return;
	}
	printf("I don't think that "); PRINTI(P);
	printf(" is a prime of the form 4n+1!\n");
}

void PELL(MPI *Dptr, MPI *Eptr)
/*
 * This program finds the period of the regular continued
 * fraction expansion of square-root d, as well as the least
 * solution x,y of the Pellian equation x*x-d*y*y=+-1.
 * The algorithm is from Sierpinski's 'Theory of Numbers', p.296.
 * and Davenport's 'The Higher Arithmetic', p.109.
 * Here sqrt(d)=a[0]+1/a[1]+...+1/a[n-1]+1/2*a[0]+1/... ,
 * The length n of the period a[1],...,a[n-1],2*a[0] is printed.
 * The partial quotients are printed iff *Eptr != 0.
 */
{
	MPI *B, *C, *Y, *Z, *A0, *A1, *L, *K, *M, *N, *H, *P, *Tmp;
	unsigned int i, l, m;
	int j;

	Y = BIG_MTHROOT(Dptr, 2);
	A0 = COPYI(Y);
	B = COPYI(Y);
	Z = MULTI(Y, Y);
	if (EQUALI(Dptr, Z))
	{
		printf("You have inputted a perfect square\n");
		printf("The program is aborted\n");
		return;
	}
	C = SUB0I(Dptr, Z);
	A1 = COPYI(C);
	if (Eptr->S)
	{
		printf("\nA[0] = "); PRINTI(Y); printf("\n");
	}
	L = ONEI(); K = COPYI(A0); M = ZEROI(); N = ONEI();
	for (i = 1; 1; i++)
	{
		FREEMPI(Z); Z = ADD0I(B, A0); 
		FREEMPI(Y); Y = INT0(Z, C);
		if (Eptr->S)
		{
			printf("A[%u] = ", i); PRINTI(Y); printf("\n");
		}
		FREEMPI(Z); Z = MULTI(Y, C);
		Tmp = B; B = SUB0I(Z, B); FREEMPI(Tmp);
		FREEMPI(Z); Z = MULTI(B, B);
		Tmp = Z; Z = SUB0I(Dptr, Z); FREEMPI(Tmp);
		Tmp = C; C = INT0(Z, C); FREEMPI(Tmp);
		FREEMPI(Z); Z = MULTI(K, Y);
		H = ADDI(Z, L);
		FREEMPI(Z); Z = MULTI(N, Y);
		P = ADDI(Z, M); FREEMPI(M);
		FREEMPI(L); L = K;
		M = N; K = H; N = P;
		if (EQUALI(B, A0))
		{
			if (EQUALI(C, A1))
			{
				printf("The continued fraction for sqrt(");
				PRINTI(Dptr);
				printf(") has period length equal to %u.\n", i);
				printf("Also the least solution (x, y) of x*x - ");
				PRINTI(Dptr);
				j = i % 2 ? -1 : 1;
				printf("y*y = %d\n",  j);
				l = LENGTHI(L); m = LENGTHI(M);
				printf(" x has %u digits;\n", l);
				printf(" y has  %u digits;\n", m);
				if (l>500)
					return;
				if (m>500)
					return;
				printf(" and \n");
				printf("x = "); PRINTI(L); printf("\n");
				printf("y = "); PRINTI(M); printf("\n");
				FREEMPI(Z); FREEMPI(L); FREEMPI(M); FREEMPI(K);
				FREEMPI(N); FREEMPI(A0); FREEMPI(A1);
				FREEMPI(B); FREEMPI(C); FREEMPI(Y);
				return;
			}
		}
		
	}
}

MPI *CONGR(MPI *A, MPI *B, MPI *M, MPI **N)
/*
 * Returns the least solution (mod N) of the congruence AX=B(mod M), where
 * N = M / gcd(A, M), otherwise returns the null pointer.
 */
{
	MPI *D, *E, *U, *V, *X, *tmp1, *tmp2;
	int t;

	D = EUCLIDI(A, M, &U, &V);
	FREEMPI(V);
	*N = INT(M, D);
	E = MOD(B, D);
	t = E->S;
	FREEMPI(E);
	if (t)
	{
		FREEMPI(U);
		FREEMPI(D);
		return ((MPI *)NULL);
	}
	tmp1 = MULTI(U, B);
	FREEMPI(U);
	tmp2 = INT(tmp1, D);
	X = MOD(tmp2, *N);
	FREEMPI(tmp1);
	FREEMPI(tmp2);
	FREEMPI(D);
	return (X);
}

MPI *CHINESE(MPI *A, MPI *B, MPI *M, MPI *N, MPI **Mptr)
/*
 * Returns the solution mod *Mptr=lcm[M,N] of the congruences X = A (mod M)
 * and X = B (mod N), if soluble; otherwise returns NULL.
 */
{
	MPI * D, *E, *F, *R, *S, *T, *U, *V;
	int t;

	D = EUCLIDI(M, N, &U, &V);
	S = SUBI(A, B);
	R = MOD(S, D);
	t = R->S;
	FREEMPI(S);
	FREEMPI(R);
	*Mptr = LCM(M, N);
	if (t)
	{
		FREEMPI(D);
		FREEMPI(U);
		FREEMPI(V);
		return ((MPI *)NULL);
	}
	S = MULTI3(B, U, M);
	R = MULTI3(A, V, N);
	FREEMPI(U);
	FREEMPI(V);
	T = ADDI(S, R);
	FREEMPI(S);
	FREEMPI(R);
	E = INT(T, D);
	FREEMPI(D);
	FREEMPI(T);
	F = MOD(E, *Mptr);
	FREEMPI(E);
	return (F);
}

MPI *CHINESE_ARRAY(MPI *A[], MPI *M[], MPI **Mptr, USI n)
/*
 * Returns the solution mod *Mptr=lcm[M[0],...,M[n-1] of the congruences
 * X = A[i] (mod M[i]),0<=i<n, if soluble; otherwise returns NULL.
 */
{
	MPI *D, *MM, *Z, *tmpMM, *tmpZ, *S, *T, *tmp;
	unsigned int i, j;
	int t;

	for (i = 0; i < n - 1; i++)
		for (j = i+1; j < n; j++)
		{
			D = GCD(M[i], M[j]);
			S = SUBI(A[i], A[j]);
			T = MOD(S, D);
			t = T->S;
			FREEMPI(D);
			FREEMPI(S);
			FREEMPI(T);
			if (t)
			{
				*Mptr = ZEROI();
				return ((MPI *)NULL);
			}
		}
	MM = COPYI(M[0]);
	Z = COPYI(A[0]);
	for (i = 1; i < n; i++)
	{
		tmpMM = MM;
		tmpZ = Z;
		Z = CHINESE(A[i], Z, M[i], MM, &tmp);
		MM = tmp;
		FREEMPI(tmpMM);
		FREEMPI(tmpZ);
	}
	*Mptr = MM;
	return (Z);
}

MPI *COLLATZT(MPI *Dptr)
{
	MPI *Eptr, *one, *Tmp;

	if (Dptr->V[0] & 1)
	{
		one = ONEI();
		Eptr = MULT_I(Dptr, 3L);
		Tmp = Eptr;
		Eptr = ADDI(Tmp, one);
		FREEMPI(Tmp);
		Tmp = Eptr;
		Eptr = INT_(Tmp, (USL)2);
		FREEMPI(Tmp);
		FREEMPI(one);
	}
	else
		Eptr = INT_(Dptr, (USL)2);
	return (Eptr);
}

void COLLATZ(MPI *Dptr, MPI *Eptr)
/* The iterates are printed iff *Eptr != 0. */
{
	unsigned int i;
	MPI *X, *Y, *Z, *Temp;

	Y = BUILDMPI(1);
	Y->S = -1;
	Y->V[0] = 5;
	Z = BUILDMPI(1);
	Z->S = -1;
	Z->V[0] = 17;

	X = COPYI(Dptr);
	if(EQZEROI(X))
	{
		printf("starting number = ");PRINTI(Dptr);printf("\n");
		printf("the number of iterations taken to reach 0 is %u\n", 0);
		FREEMPI(X);
		FREEMPI(Y);
		FREEMPI(Z);
		return;
	}
	for(i = 0; 1 ; i++)
	{
		if(EQONEI(X))
		{
			printf("starting number = ");PRINTI(Dptr);printf("\n");
			printf("the number of iterations taken to reach 1 is %u\n", i);
			break;
		}
		if(EQMINUSONEI(X))
		{
			printf("starting number = ");PRINTI(Dptr);printf("\n");
			printf("the number of iterations taken to reach -1 is %u\n", i);
			break;
		}
		if(EQUALI(X, Y))
		{
			printf("starting number = ");PRINTI(Dptr);printf("\n");
			printf("the number of iterations taken to reach -5 is %u\n", i);
			break;
		}
		if(EQUALI(X, Z))
		{
			printf("starting number = ");PRINTI(Dptr);printf("\n");
			printf("the number of iterations taken to reach -17 is %u\n", i);
			break;
		}
		Temp = X;
		X = COLLATZT(Temp);
		FREEMPI(Temp);
		if (Eptr->S)
		{
			PRINTI(X);printf("\n");
		}
	}
	FREEMPI(X);
	FREEMPI(Y);
	FREEMPI(Z);
	return;
}

MPI  *FUND_UNIT(MPI *D, MPI **Xptr, MPI **Yptr)
/*
 * This is a program for finding the fundamental unit of Q(sqrt(D)).
 * 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.
 * w=(1+sqrt(D))/2 if D=1 (mod 4), w=sqrt(D) otherwise.
 * The norm of the fundamental unit (*Xptr) + (*Yptr)*w is returned.
 */
{
	unsigned int i;
	MPI *X, *B, *C, *H, *T, *G, *Y, *F, *E, *L, *M, *N, *U, *V, *tmp;
	MPI *tmp1, *tmp2, *K, *R, *S;
	unsigned int s, t;
	
	if ((D->D == 0) && (D->V[0]) == 5)
	{
		*Xptr = ZEROI();
		*Yptr = ONEI();
		return (MINUS_ONEI());
	}
	X = BIG_MTHROOT(D, 2);
	B = COPYI(X);
	C = ONEI();
	tmp = SUB0_I(X, (USL)1); H = INT0_(tmp, (USL)2); FREEMPI(tmp);
	tmp = MULT_I(H, 2L); T = ADD0_I(tmp, (USL)1); FREEMPI(tmp);
	s = (D->V[0]) % 4;
	if (s == 1)
	{
		FREEMPI(B); B = COPYI(T);
		tmp = C; C = ADD0_I(C, (USL)1); FREEMPI(tmp); /* C = 2 */
	}
	G = MULTI(X, X); tmp = ADD0_I(G, (USL)1); FREEMPI(G);
	t = EQUALI(D, tmp); FREEMPI(tmp);
	if (t) /* period 1, exceptional case */
	{
		if (s == 1)
		{
			*Xptr = SUB0_I(X, (USL)1);
			*Yptr = TWOI();
			FREEMPI(X);
		}
		else
		{
			*Xptr = X;
			*Yptr = ONEI();
		}
		FREEMPI(B); FREEMPI(C); FREEMPI(T); FREEMPI(H);
		return (MINUS_ONEI());
	}
	tmp = MULTI(T, T); tmp1 = ADD0_I(tmp, (USL)4); FREEMPI(tmp);
	t = EQUALI(D, tmp1); FREEMPI(tmp1);
	if (t) /* period 1, exceptional case */
	{
		*Xptr = H;
		*Yptr = ONEI();
		FREEMPI(B); FREEMPI(C); FREEMPI(T); FREEMPI(X);
		return (MINUS_ONEI());
	}
	FREEMPI(T);
	L = ZEROI(); K = ONEI(); M = ONEI(); N = ZEROI();
	for (i = 0; 1; i++)
	{
		tmp = ADDI(X, B); Y = INT(tmp, C); FREEMPI(tmp);
		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);
			if ((D->V[0]) % 4 == 1)
				Y = COPYI(H);
			else
				Y = COPYI(X);
			FREEMPI(H);
		}
		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) or (sqrt(D)-1)/2 */
		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);
				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);