primes1.c

calc大数库

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


extern unsigned long PRIME[];
extern unsigned int reciprocal[][Z0];
extern unsigned int VERBOSE;
extern unsigned int FREE_PRIMES;/* FREE_PRIMES = 1 destroys the Q_[i] below */
extern unsigned int ECMAX;
#define JMAX 500


MPI *PERFECT_POWER(MPI *N)
/*
 * Here N > 1.
 * Returns X if N = X^k, for some X, k > 1,  NULL otherwise.
 * See E. Bach and J. Sorenson, "Sieve algorithms for perfect power testing"
 * Algorithmica 9 (1993) 313-328.
 */
{
	unsigned int i, t;
	int s;
	MPI *X, *Y;

	t = BINARYB(N)-1;
	i = 0;
	while (PRIME[i] <= t)
	{
		X = BIG_MTHROOT(N, (USI) PRIME[i]);
		Y = POWERI(X, (USI) PRIME[i]);
		s = RSV(Y, N);
		FREEMPI(Y);
		if (s == 0)
			return (X);
		else
		{
			i++;
			FREEMPI(X);
		}
	}
	return ((MPI *) NULL);
}

MPI *POLLARD(MPI *Nptr)
/* 
 * Pollard's p-1 method of factoring *Nptr.
 */
{
	unsigned long i, b = 1;
	MPI *T, *P, *TT, *TMP, *PP;
	
	PP = ONEI();
	T = ADD0I(PP, PP);
	while (b <= 2){
		b++;
		printf("b = %lu\n", b);
		i = 2;
	while (i  <= 10000)
	{
		if (i % 100 == 0)	
			printf("i = %lu\n", i);
		TMP = T;
		T = MPOWER_M(T, i, Nptr);
		FREEMPI(TMP);
		TT = SUB0I(T, PP);
		P = GCD(Nptr, TT);
		FREEMPI(TT);
		if (!EQONEI(P) && RSV(P, Nptr) == -1)
		{
			PRINTI(P);
			printf(" is a proper factor of ");
			PRINTI(Nptr);
			printf("\n");
			FREEMPI(PP);
			FREEMPI(T);
			return (P);
		}
		if (EQUALI(P, Nptr))
		{
			FREEMPI(P);
			b++;
			TMP = T;
			printf("GCD = *Nptr; b increased to %lu\n", b);
			T = CHANGE(b);
			FREEMPI(TMP);
			continue;
		}
		FREEMPI(P);
		i++;
	}
}
	printf("no factors returned\n");
	FREEMPI(PP);
	FREEMPI(T);
	return ((MPI *)NULL);
}

MPI *BRENT_POLLARD(MPI *Nptr)
/*
 * the Brent-Pollard method for returning a proper factor of
 * a composite MPI *Nptr. (see R. Brent, BIT 20, 176 - 184).
 */
{
	MPI  *X, *Y, *Z, *Temp, *Gptr, *PRODUCT, *Temp1;
	unsigned int g, k;
	unsigned long a, r, i, s, t;
	
	if (VERBOSE)
	{
		printf("BRENT-POLLARD IS SEARCHING FOR A PROPER FACTOR OF ");
		PRINTI(Nptr);
		printf("\n");
	}
	a = 1;
	g = 1;
	r = 1;
	PRODUCT = ONEI();
	Gptr = ONEI();
	Y = ZEROI();
	while (g == 1)
	{
		X = COPYI(Y);
		for (i = 1; i <= r; i++)
		{
			Temp = Y;
			Y = MPOWER_M(Y, (USL)2, Nptr);
			FREEMPI(Temp);
			Temp = Y;
			Y = ADD0_I(Y, a);
			FREEMPI(Temp);
			Temp = Y;
			Y = MOD0(Y, Nptr);
			FREEMPI(Temp);
		}
		k = 0;
		while (1)
		{
			Temp = Y;
			Y = MULTI(Y, Y);
			FREEMPI(Temp);
			Temp = Y;
			Y = ADD0_I(Y, a);
			FREEMPI(Temp);
			Temp = Y;
			Y = MOD0(Y, Nptr);
			FREEMPI(Temp);
			k++;
			Z = SUBI(X, Y);
			Temp1 = PRODUCT;
			PRODUCT = MULTI(PRODUCT, Z);
			FREEMPI(Z);
			FREEMPI(Temp1);
			if (k % 10 == 0)
			{
				FREEMPI(Gptr);
				Gptr = GCD(PRODUCT, Nptr);
				FREEMPI(PRODUCT);
				PRODUCT = ONEI();
				if (Gptr->D || Gptr->V[0] > 1)
				{
					g = 0;
					break;
				}
			}
			if (k >= r)
				break;
		}
		FREEMPI(X);
		r = 2 * r;
		if (VERBOSE)
			printf("r = %lu\n", r);
		s = r & 8192;
		t = EQUALI(Gptr, Nptr);
		if (s || t)
		{
			g = 1;
			if (t)
			{
				a++;
				if (VERBOSE)
					printf("a increased to %lu\n", a);
			}
 			FREEMPI(Gptr);
                        Gptr = ONEI();
			r = 1;
			FREEMPI(Y);
			Y = ZEROI();
			if (s || (a == 3))
			{
				if (VERBOSE)
				{
					printf(" BRENT_POLLARD failed to find a factor of ");
					PRINTI(Nptr);
					printf("\n");
				}
				FREEMPI(Gptr);
				FREEMPI(PRODUCT);
				FREEMPI(Y);
				return ((MPI *) NULL);
			}
		}
/*
		if (g)
			FREEMPI(Gptr);
*/
	}
	FREEMPI(Y);
	if (VERBOSE)
	{
		printf("--\n");
		printf("BRENT_POLLARD FINISHED: \n");
		PRINTI(Gptr);
		printf(" IS A PROPER FACTOR OF ");
		PRINTI(Nptr);
		printf("--\n");
	}
	FREEMPI(PRODUCT);
	return (Gptr);
}

unsigned long Q[JMAX];
unsigned int j, j_, K[JMAX], K_[JMAX], ltotal;
MPI *Q_[JMAX], *Q_PRIME[JMAX];

MPI *BABY_DIVIDE(MPI *Nptr)
/*
 * *Eptr = 1 if *Nptr is composed of primes < 1000, otherwise *Eptr
 * is the part of *Nptr composed of primes > 1000,
 * The prime factors of *Nptr < 1000 are stored in the global array Q[]
 * and the corresponding exponents in the global array K[].
 */
{
	MPI *X, *Temp;
	unsigned int a = Y0, i, k;
	unsigned long y;

	X = COPYI(Nptr);
	for (i = 0; i < a; i++)
	{
		k = 0;
		y = MOD0_(X, PRIME[i]);
		if (y == 0)
		{
 			while (y == 0)
			{
				k++;
				Temp = X;
				X = INT0_(X, PRIME[i]);
				FREEMPI(Temp);
				y = MOD0_(X, PRIME[i]);
			}
			Q[j] = PRIME[i];
			K[j] = k;
			j++;
			if (j == JMAX)
				execerror("j = JMAX, increase JMAX", "");
			if(VERBOSE){
				printf("%lu is a prime factor of ", PRIME[i]);
				PRINTI(Nptr);
				printf(" exponent: %u\n", k);
				printf("--\n");
			}
		}
	}
	return (X);
}

unsigned int MILLER(MPI *Nptr, USL b)
/*
 * *Nptr is odd, > 1, and does not divide b, 0 < b < R0.
 *  if 1 is returned, then *Nptr passes Miller's test for base b.
 * if 0 is returned, then *Nptr is composite.
 */
{
	MPI *A, *C, *D, *Temp;
	unsigned int i = 0, j = 0;

	D = SUB0_I(Nptr, (USL)1);
	A = INT0_(D, (USL)2);
	while (MOD0_(A, (USL)2) == 0)
	{
		i++;
		Temp = A;
		A = INT0_(A, (USL)2);
		FREEMPI(Temp);
	}
	C = MPOWER_((long)b, A, Nptr);
	if (C->D == 0 && C->V[0] == 1) 
	{
		FREEMPI(C);
		FREEMPI(D);
		FREEMPI(A);
		return (1);
	}
	while (1)
	{
		if (EQUALI(C, D))
		{
			FREEMPI(C);
			FREEMPI(D);
			FREEMPI(A);
			return (1);
		}
		j++;
		if (i < j)
		{
			FREEMPI(C);
			FREEMPI(D);
			FREEMPI(A);
			return (0);
		}
		Temp = C;
		C = MULTI(C, C);
		FREEMPI(Temp);
		Temp = C;
		C = MOD0(C, Nptr);
		FREEMPI(Temp);
		Temp = A;
		A = MULT_I(A, 2L);
		FREEMPI(Temp);
	}
}

unsigned int Q_PRIME_TEST(MPI *Nptr)
/*
 * *Nptr > 1 and not equal to PRIME[0],...,PRIME[4].
 * if 1 is returned, then *Nptr passes Miller's test for bases PRIME[0],
 * ...,PRIME[4] and is likely to be prime.
 * if 0 is returned, then *Nptr is composite.
 */
{
	unsigned int i;

	for (i = 0; i <= 4; i++)
	{
		if (MILLER(Nptr, PRIME[i]) == 0)
		{
			if (VERBOSE)
			{
				printf("MILLER'S TEST FINISHED: \n");
				PRINTI(Nptr);
				printf(" is composite \n");
				printf("--\n");
			}
			return (0);
		}
	}
	if (VERBOSE)
	{
		PRINTI(Nptr);
		printf(" passed MILLER'S TEST\n");
		printf("--\n");
	}
	return (1);
	
}

MPI *BIG_FACTOR(MPI *Nptr)
/*
 * *Nptr > 1 is not divisible by PRIMES[0], ..., PRIMES[167].
 * BIG_FACTOR returns a factor *Eptr > 1000 of *Nptr which is either 
 * < 1,000,000 (and hence prime) or which passes Miller's test for bases 
 * PRIMES[0],...,PRIMES[4] and is therefore likely to be prime.
 */
{
	MPI *B, *X, *Y, *Temp;
	unsigned int f;

	B = BUILDMPI(2);	
	B->V[0] = 16960;
	B->V[1] = 15;
	B->S = 1; /* *B = 1,000,000. */

	if (RSV(Nptr, B) == -1 || Q_PRIME_TEST(Nptr))
	{
		FREEMPI(B);
		return (COPYI(Nptr));
	}
	f = 1;
	X = COPYI(Nptr);
	while (f)
	{
		Y = BRENT_POLLARD(X);
		if (Y == NULL)
			Y = PERFECT_POWER(X);
		if (Y == NULL)
			Y = MPQS1(X);
		if (Y == NULL)
		{
			printf("Switching to ECF:%u elliptic curves\n", ECMAX);
			Y = EFACTOR(X,1279,1);
		}
		if (Y == NULL)
		{
			printf("Switching to POLLARD p-1\n");
			Y = POLLARD(X);
		}
		if (Y == NULL)
			execerror("no factor found", "");
		Temp = X;
		X = INT0(X, Y);
		FREEMPI(Temp);
		if (RSV(X, Y) == 1)
		{
			Temp = X;
			X = COPYI(Y); 
			FREEMPI(Temp);
		}
		if (RSV(X, B) == -1)
		{
			FREEMPI(B);
			FREEMPI(Y);
			return (X);
		}
		if (Q_PRIME_TEST(X))
			f = 0;
		FREEMPI(Y);
	}
	FREEMPI(B);
	return (X);
}

unsigned int PRIME_FACTORS(MPI *Nptr)
/*
 * A quasi-prime (q-prime) factor of *Nptr is > 1,000,000,
 * is not divisible by PRIMES[0],...,PRIMES[167], passes Millers' test
 * for bases PRIMES[0],...,PRIMES[10] and is hence likely to be prime.
 * PRIME_FACTORS returns the number of q-prime factors of *Nptr, stores
 * them in the global array Q_PRIME[].
 * Any prime factors < 1000 of *Nptr and corresponding exponents 
 * are printed and placed in the arrays Q[] and K[], while any prime factors
 * > 1000 and all q-prime factors and corresponding exponents of *Nptr are 
 * printed and placed in the arrays Q_[] and K_[].
 */
{
	MPI *B, *P, *X, *Z, *Temp;
	unsigned int k, t;

	B = BUILDMPI(2);	
	B->V[0] = 16960;
	B->V[1] = 15;
	B->S = 1;
	if (VERBOSE)
	{
		printf("FACTORIZING ");
		PRINTI(Nptr);
		printf("--\n");
	}
	X = BABY_DIVIDE(Nptr);
	t = ltotal;
	while (X->D || X->V[0] != 1)
	{
		k = 0;
		P = BIG_FACTOR(X);
		while (1)
		{
			Z = MOD0(X, P);
			if (Z->S != 0)
			{
				FREEMPI(Z);
				break;
			}
			FREEMPI(Z);
			k++;
			Temp = X;
			X = INT0(X, P);
			FREEMPI(Temp);
		}		
		if (VERBOSE)
		{
			if (RSV(P, B) == -1)
			{
				PRINTI(P);
				printf(" is a prime factor of ");
				PRINTI(Nptr);
				printf("\n");
			}
		}
		if (RSV(P, B) == 1)
		{
			if (VERBOSE)
			{
				PRINTI(P);
				printf(" is a q-prime factor of ");
				PRINTI(Nptr);
				printf("\n");
			}
			Q_PRIME[ltotal] = COPYI(P);
			ltotal++;
		}
		Q_[j_] = COPYI(P);
		FREEMPI(P);
		K_[j_] = k;
		j_++;
		if (j_ == JMAX)
			execerror("j_ = JMAX, increase JMAX", "");
		if (VERBOSE)
		{
			printf(" exponent: %u\n", k);
			printf("--\n");
		}
	}
	if (VERBOSE)
	{
		printf("FACTORIZATION INTO PRIMES AND Q-PRIMES COMPLETED FOR ");
		PRINTI(Nptr);
		printf("\n------------------------------------\n");
	}
	FREEMPI(B);
	FREEMPI(X);
	return (ltotal - t);
}

unsigned int SELFRIDGE(MPI *Nptr, USI *uptr)
/*
 * Selfridges's test for primality - see "Prime Numbers and Computer
 * Methods for Factorization" by H. Riesel, Theorem 4.4, p.106. 
 * *Nptr > 1 is a q-prime.
 * The prime and q-prime factors of n - 1 are first found. If no q-prime
 * factor is present and 1 is returned, then n is prime. However if at
 * least one q-prime factor is present and 1 is returned, then n retains its
 * q-prime status. If 0 is returned, then n is either composite or likely 
 * to be composite.
 */
{
	MPI *S, *T, *M, *N;
	unsigned int i, i_;
	unsigned long x;

	i = j;
	i_ = j_;
	M = SUB0_I(Nptr, (USL)1);
	if (VERBOSE)
	{
		printf("SELFRIDGE'S EXPONENT TEST IN PROGRESS FOR ");
		PRINTI(Nptr);
		printf("\n");
	}
	*uptr = PRIME_FACTORS(M);
	while (i <= j - 1)
	{
		for (x = 2; x < (USL)R0; x++)
		{
			S = MPOWER_((long)x, M, Nptr);
			if (S->D || S->V[0] != 1)
			{
				if (VERBOSE)
				{
				printf("SELFRIDGE'S TEST IS FINISHED:\n");
				PRINTI(Nptr);        
				printf(" is not a pseudo-prime to base %lu\n", x);
				printf(" and is hence composite\n");
				}
				FREEMPI(S);
				FREEMPI(M);
				return (0);
			}
			FREEMPI(S);
			N = INT0_(M, Q[i]);
			T = MPOWER_((long)x, N, Nptr);
			FREEMPI(N);
			if (T->D || T->V[0] != 1)
			{
				i++;
				FREEMPI(T);
				break;
			}
			FREEMPI(T);
		}
		if (x == R0)
		{
			if (VERBOSE)
			{
				printf("SELFRIDGE'S TEST IS FINISHED:\n");
				PRINTI(Nptr);
				printf(" is likely to be composite\n");
			}
			FREEMPI(M);
			return (0); 
		}
	}
	while (i_ <= j_ - 1)
	{
		for (x = 2; x < (USL)R0; x++)
		{
			S = MPOWER_((long)x, M, Nptr);
			if (S->D || S->V[0] != 1)
			{
				if (VERBOSE)
				{
				printf("SELFRIDGE'S TEST IS FINISHED:\n");
				PRINTI(Nptr);        
				printf(" is not a pseudo-prime to base %lu\n", x);
				printf(" and is hence composite\n");
				}
				FREEMPI(M);
				FREEMPI(S);
				return (0);
			}
			FREEMPI(S);
			N = INT0(M, Q_[i_]);
			T = MPOWER_((long)x, N, Nptr);
			FREEMPI(N);
			if (T->D || T->V[0] != 1)
			{
				i_++;
				FREEMPI(T);
				break;
			}
			FREEMPI(T);
		}