primes1.c

calc大数库

		if (x == R0)
		{
			if (VERBOSE)
				printf("SELFRIDGE'S TEST IS INCONCLUSIVE:\n");
			FREEMPI(M);
			return (0); 
		}
	}
	FREEMPI(M);
	if (*uptr == 0)
	{
		if (VERBOSE)
		{
			printf("SELFRIDGE'S TEST IS FINISHED:\n");
			PRINTI(Nptr);
			printf(" is prime\n");
			printf("--\n");
		}
		return (1);
	}
	else
	{
		if (VERBOSE)
		{
			printf("SELFRIDGE'S TEST IS FINISHED:\n");
			PRINTI(Nptr);
			printf(" retains its q-prime status\n");
			printf("--\n");
		}
		return (1);
	}
}

unsigned int scount; /* the number of prime factors of N < 1000 */
unsigned int s_count; /* the number of q-prime factors of N */

unsigned int FACTOR(MPI *Nptr)
/*
 * a factorization of *Nptr into prime and q-prime factors is first obtained.
 * Selfridge's primality test is then applied to any q-prime factors; the test 
 * is applied repeatedly until either a q-prime is found to be composite or
 * likely to be composite (in which case the initial factorization is doubtful 
 * and an extra base should be used in Miller's test) or we run out of q-primes,
 * in which case every q-prime factor of *Nptr is certified as prime.
 */
{
	unsigned int c, i, t, u, r, count;

	j = j_ = ltotal = 0;
	c = PRIME_FACTORS(Nptr);
	scount = j;
	s_count = j_;
	if (c == 0)
	{
		if (VERBOSE)
		{
			printf("NO Q-PRIMES:\n\n");
			PRINTI(Nptr);
			printf(" has the following factorization:\n\n");
			for (i = 0; i < scount; i++)
				printf("prime factor: %lu exponent: %u\n", Q[i], K[i]);
			for (i = 0; i < s_count; i++)
			{
				printf("prime factor: ");
				PRINTI(Q_[i]);
				printf(" exponent: %u\n", K_[i]);
			}
		}
		count = scount + s_count;
	}
	else
	{
		if (VERBOSE)
		{
			printf("TESTING Q-PRIMES FOR PRIMALITY\n");
			printf("--\n");
		}
		i = 0;
		while (c)
		{
			t = SELFRIDGE(Q_PRIME[i], &u);
			FREEMPI(Q_PRIME[i]);
			if (t == 0)
			{
				printf("do FACTOR() again with an extra base\n");
				for (r = i + 1; r < ltotal; r++)
					FREEMPI(Q_PRIME[r]);
				for (i = 0; i < j_; i++)
					FREEMPI(Q_[i]);
				return (0);
			}		
			c = c + u;
			i++;
			c--;
		}
		if (VERBOSE)
		{
			printf("ALL Q-PRIMES ARE PRIMES:\n\n");
			PRINTI(Nptr);
			printf(" has the following factorization:\n\n");
			for (i = 0; i < scount; i++)
				printf("prime factor: %lu exponent: %u\n", Q[i], K[i]);
			for (i = 0; i < s_count; i++)
			{
				printf("prime factor: ");
				PRINTI(Q_[i]);
				printf(" exponent: %u\n", K_[i]);
			}
		}
		count = scount + s_count;
	}
	/*printf("s_count = %u, j_ = %u\n", s_count, j_);*/
	for (i = s_count; i < j_; i++)
		FREEMPI(Q_[i]);
	if (FREE_PRIMES)
	{
		for (i = 0; i < s_count; i++)
			FREEMPI(Q_[i]);
	}
	return (count);
}

MPI *DIVISOR(MPI *N)
/*
 *  Returns the number of divisors of N,
 *  returns NULL if factorization unsuccessful .
 */
{
	MPI *U, *Tmp;
	unsigned int i, s;

	if (EQONEI(N))
		return (ONEI());
	U = ONEI();
	FREE_PRIMES = 1; /* This frees the Q_[i] in FACTOR()*/
	s = FACTOR(N);
	FREE_PRIMES = 0; 
	if (s == 0)
	{
		FREEMPI(U);
		return ((MPI *)NULL);
	}
	for (i = 0; i < scount; i++)
	{
		Tmp = U;
		U = MULT_I(U, 1 + (USL)K[i]);
		FREEMPI(Tmp);
	}
	for (i = 0; i < s_count; i++)
	{
		Tmp = U;
		U = MULT_I(U, 1 + (USL)K_[i]);
		FREEMPI(Tmp);
	}
	return (U);
}

MPI *MOBIUS(MPI *N)
/*
 *  Returns the Mobius function mu(N),
 *  returns NULL if factorization unsuccessful .
 */
{
	MPI *S;
	unsigned int i, s;

	if (EQONEI(N))
		return (ONEI());
	FREE_PRIMES = 1; /* This frees the Q_[i] in FACTOR()*/
	s = FACTOR(N);
	FREE_PRIMES = 0; 
	if (s == 0)
		return ((MPI *)NULL);
	for (i = 0; i < scount; i++)
		if (K[i] >= 2)
			return (ZEROI());
	for (i = 0; i < s_count; i++)
	{
		if (K_[i] >= 2)
			return (ZEROI());
	}
	if (s % 2)
		return (MINUS_ONEI());
	else
		return (ONEI());
	return (S);
}

MPI *EULER(MPI *N)
/*
 *  Returns Euler's function  phi(N),
 *  returns NULL if factorization unsuccessful .
 */
{
	MPI *U, *V, *W, *Tmp;
	unsigned int i, s;

	if (EQONEI(N))
		return (ONEI());
	U = ONEI();
	/* FREE_PRIMES = 0, so we can free the Q_[i] later */ 
	s = FACTOR(N);
	if (s == 0)
	{
		FREEMPI(U);
		return ((MPI *)NULL);
	}
	for (i = 0; i < scount; i++)
	{
		if (K[i] == 1)
			V = ONEI();
		else
			V = POWER_I((long)(Q[i]), K[i] - 1);
		W = CHANGE(Q[i] - 1);
		Tmp = U;
		U = MULTI3(U, V, W);
		FREEMPI(Tmp);
		FREEMPI(V);
		FREEMPI(W);
	}
	for (i = 0; i < s_count; i++)
	{
		if (K_[i] == 1)
			V = ONEI();
		else
			V = POWERI(Q_[i], K_[i] - 1);
		W = SUB0_I(Q_[i], 1);
		Tmp = U;
		U = MULTI3(U, V, W);
		FREEMPI(Tmp);
		FREEMPI(V);
		FREEMPI(W);
	}
	for (i = 0; i < s_count; i++)
		FREEMPI(Q_[i]);
	return (U);
}

MPI *SIGMA(MPI *N)
/*
 *  Returns sigma(N), the sum of the divisors of N,
 *  returns NULL if factorization unsuccessful .
 */
{
	MPI *U, *V, *W, *Tmp;
	unsigned int i, s;

	if (EQONEI(N))
		return (ONEI());
	U = ONEI();
	/* FREE_PRIMES = 0, so we can free the Q_[i] later */ 
	s = FACTOR(N);
	if (s == 0)
	{
		FREEMPI(U);
		return ((MPI *)NULL);
	}
	for (i = 0; i < scount; i++)
	{
		V = POWER_I((long)(Q[i]), K[i] + 1);
		Tmp = V;
		V = SUB0_I(V, 1);
		FREEMPI(Tmp);
		Tmp = V;
		V = INT0_(V, Q[i] - 1);
		FREEMPI(Tmp);
		Tmp = U;
		U = MULTI(U, V);
		FREEMPI(Tmp);
		FREEMPI(V);
	}
	for (i = 0; i < s_count; i++)
	{	V = POWERI(Q_[i], K_[i] + 1);
		Tmp = V;
		V = SUB0_I(V, 1);
		FREEMPI(Tmp);
		Tmp = V;
		W = SUB0_I(Q_[i], 1);
		V = INT0(V, W);
		FREEMPI(W);
		FREEMPI(Tmp);
		Tmp = U;
		U = MULTI(U, V);
		FREEMPI(Tmp);
		FREEMPI(V);
	}
	for (i = 0; i < s_count; i++)
		FREEMPI(Q_[i]);
	return (U);
}

MPI *LPRIMROOT(MPI *P)
/*
 *  Returns  the least primitive root mod P, an odd prime;
 *  returns NULL if factorization of P - 1 is unsuccessful .
 */
{
	MPI *Tmp, *QQ, *R, *RR, *T;
	unsigned int s, u, success = 1;
	long i;
	int r;

	T = ZEROI();
	QQ = SUB0_I(P, 1);
	s = FACTOR(QQ);
	if (s == 0)
		return (T);
	for (i = 2; i >= 2; i++)
	{
		RR = CHANGEI(i);
		r = JACOBIB(RR, P);
		FREEMPI(RR);
		if (r == -1)
		{
			for (u = 0; u < scount; u++)
			{
				Tmp = INT0_(QQ, Q[u]);
				R = MPOWER_(i, Tmp, P);
				FREEMPI(Tmp);
				if (EQONEI(R))
				{
					success = 0;
					FREEMPI(R);
					break;
				}
				FREEMPI(R);
			}	
			if (success)
			{
				for (u = 0; u < s_count; u++)
				{
					Tmp = INT0(QQ, Q_[u]);
					R = MPOWER_(i, Tmp, P);
					FREEMPI(Tmp);
					if (EQONEI(R))
					{
						success = 0;
						FREEMPI(R);
						break;
					}
					FREEMPI(R);
				}	
			}
			if (success)
			{
				FREEMPI(QQ);
				for (u = 0; u < s_count; u++)
					FREEMPI(Q_[u]);
				FREEMPI(T);
				T = CHANGE(i);
				break;
			}
			else
				success  = 1;
		}
	}
	return (T);
}

MPI *ORDERP(MPI *A, MPI *P)
/*
 * Returns the order of A mod P, a prime.
 */
{
	MPI *AA, *Tmp, *QQ, *M;
	unsigned int s, u, i;

	AA = MOD(A, P);
	Tmp = SUB0_I(AA, 1);
	if (Tmp->S == 0) 
	{
		FREEMPI(Tmp);
		FREEMPI(AA);
		return (ONEI());
	}
	FREEMPI(Tmp);
	Tmp = ADD0_I(AA, 1);
	if (Tmp->S == 0) 
	{
		FREEMPI(Tmp);
		FREEMPI(AA);
		return (CHANGE(2));
	}
	FREEMPI(Tmp);
	QQ = SUB0_I(P, 1);
	/* FREE_PRIMES = 0, so we can free the Q_[i] later */ 
	s = FACTOR(QQ);
	for (u = 0; u < scount; u++)
	{
		for (i = 1; i <= K[u]; i++)
		{
			Tmp = INT0_(QQ, Q[u]);
			M = MPOWER(AA, Tmp, P);
			FREEMPI(Tmp);
			if (!EQONEI(M))
			{
				FREEMPI(M);
				break;
			}
			FREEMPI(M);
			Tmp = QQ;
			QQ = INT0_(QQ, Q[u]);
			FREEMPI(Tmp);
		}
	}
	for (u = 0; u < s_count; u++)
	{
		for (i = 1; i <= K_[u]; i++)
		{
			Tmp = INT0(QQ, Q_[u]);
			M = MPOWER(AA, Tmp, P);
			FREEMPI(Tmp);
			if (!EQONEI(M))
			{
				FREEMPI(M);
				break;
			}
			FREEMPI(M);
			Tmp = QQ;
			QQ = INT0_(QQ, Q[u]);
			FREEMPI(Tmp);
		}
	}
	FREEMPI(AA);
	for (i = 0; i < s_count; i++)
		FREEMPI(Q_[i]);
	return (QQ);
}

MPI *ORDERQ(MPI *A, MPI *P, unsigned int n)
/*
 * Returns the order of A mod P^n, P a prime.
 */
{
	MPI *D, *E, *Tmp, *Tmp1, *Q;
	unsigned int h;

	if (EQONEI(A))
		return (ONEI());
	if (EQMINUSONEI(A))
		return (CHANGE(2));
	if (n == 1)
		return (ORDERP(A, P));
	if (P->D == 0 && P->V[0] == 2)
	{
		if ((A->V[0] - 1) % 4 == 0)
			D = ONEI();
		/*if ((A->V[0] + 1) % 4 == 0)*/
		else
			D = CHANGE(2);
	}
	else
		D = ORDERP(A, P);
	Q = COPYI(P);
	E = ZEROI();
	h = 0;
	while (E->S == 0)
	{
		Tmp = Q;
		Q = MULTI(Q, P);
		FREEMPI(Tmp);
		Tmp = E;
		E = MPOWER(A, D, Q);	
		FREEMPI(Tmp);
		Tmp = E;
		E = SUB0_I(E, 1);
		FREEMPI(Tmp);
		h++;
	}
	FREEMPI(E);
	FREEMPI(Q);
	if (n <= h)
		return (D);
	else
	{
		Tmp = POWERI(P, n - h);
		Tmp1 = MULTI(Tmp, D);
		FREEMPI(Tmp);
		FREEMPI(D);
		return (Tmp1);
	}
}

MPI *ORDERM(MPI *A, MPI *M)
/*
 * Returns the order of A mod M.
 */
{
	MPI *O, **Y, *Tmp, *Tmp1;
	unsigned int i, s, *x;

	if (EQONEI(A))
		return (ONEI());
	if (EQMINUSONEI(A))
	{
		if (M->D == 0 && M->V[0] == 2)
			return (ONEI());
		else
			return (CHANGE(2));
	}
	/* FREE_PRIMES = 0, so we can free the Q_[i] later */ 
	s = FACTOR(M);
	x = (USI *)mmalloc(s * sizeof(USI));
	Y = (MPI **)mmalloc(s * sizeof(MPI *));
	for (i = 0; i < scount; i++)
	{
		x[i] = K[i];
		Y[i] = CHANGE(Q[i]);
	}
	for (i = 0; i < s_count; i++)/* bugfix here - ANU,27/9/94. */ 
	{
		x[i + scount] = K_[i];
		Y[i + scount] = COPYI(Q_[i]);
		FREEMPI(Q_[i]);
	}
	O = ORDERQ(A, Y[0], x[0]);
	for (i = 1; i < s; i++)
	{
		Tmp = O;
		Tmp1 = ORDERQ(A, Y[i], x[i]);
		O = LCM(O, Tmp1);
		FREEMPI(Tmp);
		FREEMPI(Tmp1);
	}
	ffree(x, s * sizeof(USI));
	for (i = 0; i < s; i++)
		FREEMPI(Y[i]);
	ffree(Y, s * sizeof(MPI *));
	return (O);
}

void FERMAT_QUOTIENT(MPI *N)
{
	USI i, n;
	MPI *P, *Q, *O, *T1, *T2, *T3, *T4, *T33, *T44;
	MPI *PP, *A, *B, *R;

	n = CONVERTI(N);
	O = ONEI();
	for (i = 2; i <= n; i++)
	{
		printf("i = %u\n", i);
		R = CHANGE(PRIME[i]);
		P = CHANGE(PRIME[i]+ 1);
		Q = CHANGE(PRIME[i]- 1);
		T1 = POWERI(P, (USI)PRIME[i]);
		FREEMPI(P);
		T2 = POWERI(Q, (USI)PRIME[i]);
		FREEMPI(Q);
		T3 = SUB0I(T1, O);
		T4 = ADD0I(T2, O);
		FREEMPI(T1);
		FREEMPI(T2);
		PP = MULTI(R, R);
		FREEMPI(R);
		T33 = INT0(T3, PP);
		T44 = INT0(T4, PP);
		FREEMPI(T3);
		FREEMPI(T4);
		FREEMPI(PP);
		A = LUCAS(T33);
		FREEMPI(T33);
		if (A->S)
			printf("p = %lu gives a - prime quotient\n", PRIME[i]);
		FREEMPI(A);
		B = LUCAS(T44);
		FREEMPI(T44);
		if (B->S)
			printf("p = %lu gives a + prime quotient\n", PRIME[i]);
		FREEMPI(B);
	}
	FREEMPI(O);
	return;
}

MPI *SQROOT1(MPI *A, MPI *P, USL n)
{
/* Solving x^2=A (mod P^n), P an odd prime not dividing A. */
/* returns -1 if not soluble, otherwise returns x. */

	MPI *T, *K, *U, *V, *Tmp1, *Tmp2, *N;
	USI i, t;

	Tmp1 = INT_(P, 2);
	T= MPOWER(A, Tmp1, P); 
	t = EQONEI(T);
	FREEMPI(T);
	FREEMPI(Tmp1);
	if (t == 0)
		return(MINUS_ONEI());
	if (n == 1)
		return(SQRTM(A, P));
	else{
		U=SQROOT1(A, P, n - 1);
		Tmp1 = MULTI(U, U);
		V = SUBI(Tmp1, A);
		FREEMPI(Tmp1);
		for(i = 0; i < n - 1; i++){
			Tmp1 = V;
			V = INT(V, P);
			FREEMPI(Tmp1);
		}
		Tmp1 = MULT_I(U, 2);
		Tmp2 = MINUSI(V);
		FREEMPI(V);
		K = CONGR(Tmp1, Tmp2, P, &N);
		FREEMPI(Tmp1);
		FREEMPI(Tmp2);
		FREEMPI(N);
		Tmp1 = POWERI(P, (USI)(n - 1));
		Tmp2 = MULTI(K, Tmp1);
		FREEMPI(K);
		FREEMPI(Tmp1);
		Tmp1 = ADDI(U, Tmp2);
		FREEMPI(U);
		FREEMPI(Tmp2);
		return(Tmp1);
	}
}

void SQROOT1_W(Stack S)