zprime.c

Calc Software Package for Number Calc

 *
 * NOTE: This routine will not return a factor == the test value
 *	 except when the test value is 1 or -1.
 */
FLAG
zfactor(ZVALUE n, ZVALUE zlimit, ZVALUE *res)
{
	FULL f;			/* factor found, or 0 */

	/*
	 * determine the limit
	 */
	if (zge32b(zlimit)) {
		/* limit is too large to be reasonable */
		return -1;
	}
	n.sign = 0;		/* ignore sign of n */

	/*
	 * find the smallest factor <= limit, if possible
	 */
	f = small_factor(n, ztofull(zlimit));

	/*
	 * report the results
	 */
	if (f > 0) {
		/* return factor if requested */
		if (res) {
			utoz(f, res);
		}
		/* report a factor was found */
		return 1;
	}
	/* no factor was found */
	return 0;
}


/*
 * Find a smallest prime factor <= some small (32 bit) limit of a value
 *
 * given:
 *	z		number to factor
 *	limit		largest factor we will test
 *
 * Returns:
 *	0	no prime <= the limit was found
 *	!= 0	the smallest prime factor
 */
static FULL
small_factor(ZVALUE z, FULL limit)
{
	FULL top;			/* current max factor level */
	CONST unsigned short *tp;	/* pointer to a tiny prime */
	FULL factlim;			/* highest factor to test */
	CONST unsigned short *p;	/* test factor */
	FULL factor;			/* test factor */
	HALF tlim;			/* limit on prime table use */
	HALF divval[2];			/* divisor value */
	ZVALUE div;			/* test factor/divisor */
	ZVALUE tmp;
	CONST unsigned char *j;

	/*
	 * catch impossible ranges
	 */
	if (limit < 2) {
		/* range is too small */
		return 0;
	}

	/*
	 * perform the even test
	 */
	if (ziseven(z)) {
		if (zistwo(z)) {
			/* z is 2, so don't return 2 as a factor */
			return 0;
		}
		return 2;

	/*
	 * value is odd
	 */
	} else if (limit == 2) {
		/* limit is 2, value is odd, no factors will ever be found */
		return 0;
	}

	/*
	 * force the factor limit to be odd
	 */
	if ((limit & 0x1) == 0) {
		--limit;
	}

	/*
	 * case: number to factor fits into a FULL
	 */
	if (!zgtmaxufull(z)) {
		FULL val = ztofull(z);	/* find the smallest factor of val */
		FULL isqr;		/* sqrt of val */

		/*
		 * special case: val is a prime <= MAX_MAP_PRIME
		 */
		if (val <= MAX_MAP_PRIME && pr_map_bit(val)) {
			/* z is prime, so no factors will be found */
			return 0;
		}

		/*
		 * we need not search above the sqrt of val
		 */
		isqr = fsqrt(val);
		if (limit > isqr) {
			/* limit is largest odd value <= sqrt of val */
			limit = ((isqr & 0x1) ? isqr : isqr-1);
		}

		/*
		 * search for a small prime factor
		 */
		top = ((limit < MAX_MAP_VAL) ? limit : MAX_MAP_VAL);
		for (tp = prime; *tp <= top && *tp != 1; ++tp) {
			if (val%(*tp) == 0) {
				return ((FULL)*tp);
			}
		}

#if FULL_BITS == 64
		/*
		 * Our search will carry us beyond the prime table.  We will
		 * continue to values until we reach our limit or until a
		 * factor is found.
		 *
		 * It is faster to simply test odd values and ignore non-prime
		 * factors because the work needed to find the next prime is
		 * more than the work one saves in not factor with non-prime
		 * values.
		 *
		 * We can improve on this method by skipping odd values that
		 * are a multiple of 3, 5, 7 and 11.  We use a table of
		 * bytes that indicate the offsets between odd values that
		 * are not a multiple of 3,4,5,7 & 11.
		 */
		/* XXX - speed up test for large z by using gcds */
		j = jmp + jmpptr(NXT_MAP_PRIME);
		for (top=NXT_MAP_PRIME; top <= limit; top += nxtjmp(j)) {
			if ((val % top) == 0) {
				return top;
			}
		}
#endif /* FULL_BITS == 64 */

		/* no prime factors found */
		return 0;
	}

	/*
	 * Find a factor of a value that is too large to fit into a FULL.
	 *
	 * determine if/what our sqrt factor limit will be
	 */
	if (zge64b(z)) {
		/* we have no factor limit, avoid highest factor */
		factlim = MAX_SM_PRIME-1;
	} else if (zge32b(z)) {
		/* determine if limit is too small to matter */
		if (limit < BASE) {
			factlim = limit;
		} else {
			/* find the isqrt(z) */
			if (!zsqrt(z, &tmp, 0)) {
				/* sqrt is exact */
				factlim = ztofull(tmp);
			} else {
				/* sqrt is inexact */
				factlim = ztofull(tmp)+1;
			}
			zfree(tmp);

			/* avoid highest factor */
			if (factlim >= MAX_SM_PRIME) {
				factlim = MAX_SM_PRIME-1;
			}
		}
	} else {
		/* determine our factor limit */
		factlim = fsqrt(ztofull(z));
		if (factlim >= MAX_SM_PRIME) {
			factlim = MAX_SM_PRIME-1;
		}
	}
	if (factlim > limit) {
		factlim = limit;
	}

	/*
	 * walk the prime table looking for factors
	 *
	 * XXX - consider using gcd of products of primes to speed this
	 *	 section up
	 */
	tlim = (HALF)((factlim >= MAX_MAP_PRIME) ? MAX_MAP_PRIME-1 : factlim);
	div.sign = 0;
	div.v = divval;
	div.len = 1;
	for (p=prime; (HALF)*p <= tlim; ++p) {

		/* setup factor */
		div.v[0] = (HALF)(*p);

		if (zdivides(z, div))
			return (FULL)(*p);
	}
	if ((FULL)*p > factlim) {
		/* no factor found */
		return (FULL)0;
	}

	/*
	 * test the highest factor possible
	 */
	div.v[0] = MAX_MAP_PRIME;

	if (zdivides(z, div))
		return (FULL)MAX_MAP_PRIME;

	/*
	 * generate higher test factors as needed
	 *
	 * XXX - consider using gcd of products of primes to speed this
	 *	 section up
	 */
#if BASEB == 16
	div.len = 2;
#endif
	factor = NXT_MAP_PRIME;
	j = jmp + jmpptr(factor);
	for(; factor <= factlim; factor += nxtjmp(j)) {

		/* setup factor */
#if BASEB == 32
		div.v[0] = (HALF)factor;
#else
		div.v[0] = (HALF)(factor & BASE1);
		div.v[1] = (HALF)(factor >> BASEB);
#endif

		if (zdivides(z, div))
			return (FULL)(factor);
	}
	if (factor >= factlim) {
		/* no factor found */
		return (FULL)0;
	}

	/*
	 * test the highest factor possible
	 */
#if BASEB == 32
	div.v[0] = MAX_SM_PRIME;
#else
	div.v[0] = (MAX_SM_PRIME & BASE1);
	div.v[1] = (MAX_SM_PRIME >> BASEB);
#endif
	if (zdivides(z, div))
		return (FULL)MAX_SM_PRIME;

	/*
	 * no factor found
	 */
	return (FULL)0;
}


/*
 * Compute the product of the primes up to the specified number.
 */
void
zpfact(ZVALUE z, ZVALUE *dest)
{
	long n;				/* limiting number to multiply by */
	long p;				/* current prime */
	CONST unsigned short *tp;	/* pointer to a tiny prime */
	CONST unsigned char *j;		/* current jump increment */
	ZVALUE res, temp;

	/* firewall */
	if (zisneg(z)) {
		math_error("Negative argument for factorial");
		/*NOTREACHED*/
	}
	if (zge24b(z)) {
		math_error("Very large factorial");
		/*NOTREACHED*/
	}
	n = ztolong(z);

	/*
	 * Deal with table lookup pfact values
	 */
	if (n <= MAX_PFACT_VAL) {
		utoz(pfact_tbl[n], dest);
		return;
	}

	/*
	 * Multiply by the primes in the static table
	 */
	utoz(pfact_tbl[MAX_PFACT_VAL], &res);
	for (tp=(&prime[NXT_PFACT_VAL]); *tp != 1 && (long)(*tp) <= n; ++tp) {
		zmuli(res, *tp, &temp);
		zfree(res);
		res = temp;
	}

	/*
	 * if needed, multiply by primes beyond the static table
	 */
	j = jmp + jmpptr(NXT_MAP_PRIME);
	for (p = NXT_MAP_PRIME; p <= n; p += nxtjmp(j)) {
		FULL isqr;	/* isqrt(p) */

		/* our factor limit - see next_prime for why this works */
		isqr = fsqrt(p)+1;
		if ((isqr & 0x1) == 0) {
			--isqr;
		}

		/* ignore Saber-C warning #530 about empty for statement */
		/*	  ok to ignore in proc zpfact */
		/* find the next prime */
		for (tp=prime; (*tp <= isqr) && (p % (long)(*tp)); ++tp) {
		}
		if (*tp <= isqr && *tp != 1) {
			continue;
		}

		/* multiply by the next prime */
		zmuli(res, p, &temp);
		zfree(res);
		res = temp;
	}
	*dest = res;
}


/*
 * Perform a probabilistic primality test (algorithm P in Knuth vol2, 4.5.4).
 * Returns FALSE if definitely not prime, or TRUE if probably prime.
 * Count determines how many times to check for primality.
 * The chance of a non-prime passing this test is less than (1/4)^count.
 * For example, a count of 100 fails for only 1 in 10^60 numbers.
 *
 * It is interesting to note that ptest(a,1,x) (for any x >= 0) of this
 * test will always return TRUE for a prime, and rarely return TRUE for
 * a non-prime.	 The 1/4 is appears in practice to be a poor upper
 * bound.  Even so the only result that is EXACT and TRUE is when
 * this test returns FALSE for a non-prime.  When ptest returns TRUE,
 * one cannot determine if the value in question is prime, or the value
 * is one of those rare non-primes that produces a false positive.
 *
 * The absolute value of count determines how many times to check
 * for primality.  If count < 0, then the trivial factor check is
 * omitted.
 * skip = 0 uses random bases
 * skip = 1 uses prime bases 2, 3, 5, ...
 * skip > 1 or < 0 uses bases skip, skip + 1, ...
 */
BOOL
zprimetest(ZVALUE z, long count, ZVALUE skip)
{
	long limit = 0;		/* test odd values from skip up to limit */
	ZVALUE zbase;		/* base as a ZVALUE */
	long i, ij, ik;
	ZVALUE zm1, z1, z2, z3;
	int type;		/* random, prime or consecutive integers */
	CONST unsigned short *pr;	/* pointer to small prime */

	/*
	 * firewall - ignore sign of z, values 0 and 1 are not prime
	 */
	z.sign = 0;
	if (zisleone(z)) {
		return 0;
	}

	/*
	 * firewall - All even values, except 2, are not prime
	 */
	if (ziseven(z))
		return zistwo(z);

	if (z.len == 1 && *z.v == 3)
		return 1;		/* 3 is prime */

	/*
	 * we know that z is an odd value > 1
	 */

	/*
	 * Perform trivial checks if count is not negative
	 */
	if (count >= 0) {

		/*
		 * If the number is a small (32 bit) value, do a direct test
		 */
		if (!zge32b(z)) {
			return zisprime(z);
		}

		/*
		 * See if the number has a tiny factor.
		 */
		if (small_factor(z, PTEST_PRECHECK) != 0) {
			/* a tiny factor was found */
			return FALSE;
		}

		/*
		 * If our count is zero, do nothing more
		 */
		if (count == 0) {
			/* no test was done, so no test failed! */
			return TRUE;
		}

	} else {
		/* use the absolute value of count */
		count = -count;
	}
	if (z.len < conf->redc2) {
		return zredcprimetest(z, count, skip);
	}

	if (ziszero(skip)) {
		type = 0;
		zbase = _zero_;
	} else if (zisone(skip)) {
		type = 1;
		itoz(2, &zbase);
		limit = 1 << 16;
		if (!zge16b(z))
			limit = ztolong(z);
	} else {
		type = 2;
		if (zrel(skip, z) >= 0 || zisneg(skip))
			zmod(skip, z, &zbase, 0);
		else
			zcopy(skip, &zbase);
	}
	/*
	 * Loop over various bases, testing each one.
	 */
	zsub(z, _one_, &zm1);
	ik = zlowbit(zm1);
	zshift(zm1, -ik, &z1);
	pr = prime;
	for (i = 0; i < count; i++) {
		switch (type) {
		case 0:
			zfree(zbase);
			zrandrange(_two_, zm1, &zbase);
			break;
		case 1:
			if (i == 0)
				break;
			zfree(zbase);
			if (*pr == 1 || (long)*pr >= limit) {
				zfree(z1);
				zfree(zm1);
				return TRUE;
			}
			itoz((long) *pr++, &zbase);
			break;
		default:
			if (i == 0)
				break;
			zadd(zbase, _one_, &z3);
			zfree(zbase);
			zbase = z3;
		}

		ij = 0;
		zpowermod(zbase, z1, z, &z3);
		for (;;) {
			if (zisone(z3)) {
				if (ij) {
				/* number is definitely not prime */
					zfree(z3);
					zfree(zm1);
					zfree(z1);
					zfree(zbase);
					return FALSE;
				}
				break;
			}
			if (!zcmp(z3, zm1))
				break;
			if (++ij >= ik) {
				/* number is definitely not prime */
				zfree(z3);
				zfree(zm1);
				zfree(z1);
				zfree(zbase);
				return FALSE;
			}
			zsquare(z3, &z2);
			zfree(z3);
			zmod(z2, z, &z3, 0);
			zfree(z2);
		}
		zfree(z3);
	}
	zfree(zm1);
	zfree(z1);
	zfree(zbase);

	/* number might be prime */
	return TRUE;
}


/*
 * Called by zprimetest when simple cases have been eliminated
 * and z.len < conf->redc2.  Here count > 0, z is odd and > 3.
 */
BOOL
zredcprimetest(ZVALUE z, long count, ZVALUE skip)
{
	long limit = 0;		/* test odd values from skip up to limit */
	ZVALUE zbase;		/* base as a ZVALUE */
	REDC *rp;
	long i, ij, ik;
	ZVALUE zm1, z1, z2, z3;
	ZVALUE zredcm1;