zprime.c

Calc Software Package for Number Calc

	int type;		/* random, prime or consecutive integers */
	CONST unsigned short *pr;	/* pointer to small prime */


	rp = zredcalloc(z);
	zsub(z, rp->one, &zredcm1);
	if (ziszero(skip)) {
		zbase = _zero_;
		type = 0;
	} else if (zisone(skip)) {
		itoz(2, &zbase);
		type = 1;
		limit = 1 << 16;
		if (!zge16b(z))
			limit = ztolong(z);
	} else {
		zredcencode(rp, skip, &zbase);
		type = 2;
	}
	/*
	 * Loop over various "random" numbers, 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:
			do {
				zfree(zbase);
				zrandrange(_one_, z, &zbase);
			}
			while (!zcmp(zbase, rp->one) ||
					!zcmp(zbase, zredcm1));
			break;
		case 1:
			if (i == 0) {
				break;
			}
			zfree(zbase);
			if (*pr == 1 || (long)*pr >= limit) {
				zfree(z1);
				zfree(zm1);
				if (z.len < conf->redc2) {
					zredcfree(rp);
					zfree(zredcm1);
				}
				return TRUE;
			}
			itoz((long) *pr++, &z3);
			zredcencode(rp, z3, &zbase);
			zfree(z3);
			break;
		default:
			if (i == 0)
				break;
			zadd(zbase, rp->one, &z3);
			zfree(zbase);
			zbase = z3;
			if (zrel(zbase, z) >= 0) {
				zsub(zbase, z, &z3);
				zfree(zbase);
				zbase = z3;
			}
		}

		ij = 0;
		zredcpower(rp, zbase, z1, &z3);
		for (;;) {
			if (!zcmp(z3, rp->one)) {
				if (ij) {
				/* number is definitely not prime */
					zfree(z3);
					zfree(zm1);
					zfree(z1);
					zfree(zbase);
					zredcfree(rp);
					zfree(zredcm1);
					return FALSE;
				}
				break;
			}
			if (!zcmp(z3, zredcm1))
				break;
			if (++ij >= ik) {
				/* number is definitely not prime */
				zfree(z3);
				zfree(zm1);
				zfree(z1);
				zfree(zbase);
				zredcfree(rp);
				zfree(zredcm1);
				return FALSE;
			}
			zredcsquare(rp, z3, &z2);
			zfree(z3);
			z3 = z2;
		}
		zfree(z3);
	}
	zfree(zbase);
	zredcfree(rp);
	zfree(zredcm1);
	zfree(zm1);
	zfree(z1);

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


/*
 * znextcand - find the next integer that passes ptest().
 * The signs of z and mod are ignored.	Result is the least integer
 * greater than abs(z) congruent to res modulo abs(mod), or if there
 * is no such integer, zero.
 *
 * given:
 *	z		search point > 2
 *	count		ptests to perform per candidate
 *	skip		ptests to skip
 *	res		return congruent to res modulo abs(mod)
 *	mod		congruent to res modulo abs(mod)
 *	cand		candidate found
 */
BOOL
znextcand(ZVALUE z, long count, ZVALUE skip, ZVALUE res, ZVALUE mod, ZVALUE *cand)
{
	ZVALUE tmp1;
	ZVALUE tmp2;

	z.sign	= 0;
	mod.sign = 0;
	if (ziszero(mod)) {
		if (zrel(res, z) > 0 && zprimetest(res, count, skip)) {
			zcopy(res, cand);
			return TRUE;
		}
		return FALSE;
	}
	if (ziszero(z) && zisone(mod)) {
		zcopy(_two_, cand);
		return TRUE;
	}
	zsub(res, z, &tmp1);
	if (zmod(tmp1, mod, &tmp2, 0))
		zadd(z, tmp2, cand);
	else
		zadd(z, mod, cand);

	/*
	 * Now *cand is least integer greater than abs(z) and congruent
	 * to res modulo mod.
	 */
	zfree(tmp1);
	zfree(tmp2);
	if (zprimetest(*cand, count, skip))
		return TRUE;
	zgcd(*cand, mod, &tmp1);
	if (!zisone(tmp1)) {
		zfree(tmp1);
		zfree(*cand);
		return FALSE;
	}
	zfree(tmp1);
	if (ziseven(*cand)) {
		zadd(*cand, mod, &tmp1);
		zfree(*cand);
		*cand = tmp1;
		if (zprimetest(*cand, count, skip))
			return TRUE;
	}
	/*
	 * *cand is now least odd integer > abs(z) and congruent to
	 * res modulo mod.
	 */
	if (zisodd(mod))
		zshift(mod, 1, &tmp1);
	else
		zcopy(mod, &tmp1);
	do {
		zadd(*cand, tmp1, &tmp2);
		zfree(*cand);
		*cand = tmp2;
	} while (!zprimetest(*cand, count, skip));
	zfree(tmp1);
	return TRUE;
}


/*
 * zprevcand - find the nearest previous integer that passes ptest().
 * The signs of z and mod are ignored.	Result is greatest positive integer
 * less than abs(z) congruent to res modulo abs(mod), or if there
 * is no such integer, zero.
 *
 * given:
 *	z		search point > 2
 *	count		ptests to perform per candidate
 *	skip		ptests to skip
 *	res		return congruent to res modulo abs(mod)
 *	mod		congruent to res modulo abs(mod)
 *	cand		candidate found
 */
BOOL
zprevcand(ZVALUE z, long count, ZVALUE skip, ZVALUE res, ZVALUE mod, ZVALUE *cand)
{
	ZVALUE tmp1;
	ZVALUE tmp2;

	z.sign	= 0;
	mod.sign = 0;
	if (ziszero(mod)) {
		if (zispos(res)&&zrel(res, z)<0 && zprimetest(res,count,skip)) {
			zcopy(res, cand);
			return TRUE;
		}
		return FALSE;
	}
	zsub(z, res, &tmp1);
	if (zmod(tmp1, mod, &tmp2, 0))
		zsub(z, tmp2, cand);
	else
		zsub(z, mod, cand);
	/*
	 * *cand is now the greatest integer < z that is congruent to res
	 * modulo mod.
	 */
	zfree(tmp1);
	zfree(tmp2);
	if (zisneg(*cand)) {
		zfree(*cand);
		return FALSE;
	}
	if (zprimetest(*cand, count, skip))
		return TRUE;
	zgcd(*cand, mod, &tmp1);
	if (!zisone(tmp1)) {
		zfree(tmp1);
		zmod(*cand, mod, &tmp1, 0);
		zfree(*cand);
		if (zprimetest(tmp1, count, skip)) {
			*cand = tmp1;
			return TRUE;
		}
		if (ziszero(tmp1)) {
			zfree(tmp1);
			if (zprimetest(mod, count, skip)) {
				zcopy(mod, cand);
				return TRUE;
			}
			return FALSE;
		}
		zfree(tmp1);
		return FALSE;
	}
	zfree(tmp1);
	if (ziseven(*cand)) {
		zsub(*cand, mod, &tmp1);
		zfree(*cand);
		if (zisneg(tmp1)) {
			zfree(tmp1);
			return FALSE;
		}
		*cand = tmp1;
		if (zprimetest(*cand, count, skip))
			return TRUE;
	}
	/*
	 * *cand is now the greatest odd integer < z that is congruent to
	 * res modulo mod.
	 */
	if (zisodd(mod))
		zshift(mod, 1, &tmp1);
	else
		zcopy(mod, &tmp1);

	do {
		zsub(*cand, tmp1, &tmp2);
		zfree(*cand);
		*cand = tmp2;
	} while (!zprimetest(*cand, count, skip) && !zisneg(*cand));
	zfree(tmp1);
	if (zisneg(*cand)) {
			zadd(*cand, mod, &tmp1);
			zfree(*cand);
			*cand = tmp1;
			if (zistwo(*cand))
				return TRUE;
			zfree(*cand);
			return FALSE;
	}
	return TRUE;
}


/*
 * Find the lowest prime factor of a number if one can be found.
 * Search is conducted for the first count primes.
 *
 * Returns:
 *	1	no factor found or z < 3
 *	>1	factor found
 */
FULL
zlowfactor(ZVALUE z, long count)
{
	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;

	z.sign = 0;

	/*
	 * firewall
	 */
	if (count <= 0 || zisleone(z) || zistwo(z)) {
		/* number is < 3 or count is <= 0 */
		return (FULL)1;
	}

	/*
	 * test for the first factor
	 */
	if (ziseven(z)) {
		return (FULL)2;
	}
	if (count <= 1) {
		/* count was 1, tested the one and only factor */
		return (FULL)1;
	}

	/*
	 * 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)) {
		/* 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;
	}

	/*
	 * walk the prime table looking for factors
	 */
	tlim = (HALF)((factlim >= MAX_MAP_PRIME) ? MAX_MAP_PRIME-1 : factlim);
	div.sign = 0;
	div.v = divval;
	div.len = 1;
	for (p=prime, --count; count > 0 && (HALF)*p <= tlim; ++p, --count) {

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

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

	/*
	 * 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
	 */
#if BASEB == 16
	div.len = 2;
#endif
	for(factor = NXT_MAP_PRIME;
	    count > 0 && factor <= factlim;
	    factor = next_prime(factor), --count) {

		/* 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 (count <= 0 || factor >= factlim) {
		/* no factor found */
		return (FULL)1;
	}

	/*
	 * 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)1;
}


/*
 * Compute the least common multiple of all the numbers up to the
 * specified number.
 */
void
zlcmfact(ZVALUE z, ZVALUE *dest)
{
	long n;				/* limiting number to multiply by */
	long p;				/* current prime */
	long pp = 0;			/* power of prime */
	long i;				/* test value */
	CONST unsigned short *pr;	/* pointer to a small prime */
	ZVALUE res, temp;

	if (zisneg(z) || ziszero(z)) {
		math_error("Non-positive argument for lcmfact");
		/*NOTREACHED*/
	}
	if (zge24b(z)) {
		math_error("Very large lcmfact");
		/*NOTREACHED*/
	}
	n = ztolong(z);
	/*
	 * Multiply by powers of the necessary odd primes in order.
	 * The power for each prime is the highest one which is not
	 * more than the specified number.
	 */
	res = _one_;
	for (pr=prime; (long)(*pr) <= n && *pr > 1; ++pr) {
		i = p = *pr;
		while (i <= n) {
			pp = i;
			i *= p;
		}
		zmuli(res, pp, &temp);
		zfree(res);
		res = temp;
	}
	for (p = NXT_MAP_PRIME; p <= n; p = (long)next_prime(p)) {
		i = p;
		while (i <= n) {
			pp = i;
			i *= p;
		}
		zmuli(res, pp, &temp);
		zfree(res);
		res = temp;
	}
	/*
	 * Finish by scaling by the necessary power of two.
	 */
	zshift(res, zhighbit(z), dest);
	zfree(res);
}


/*
 * fsqrt - fast square root of a FULL value
 *
 * We will determine the square root of a given value.
 * Starting with the integer square root of the largest power of
 * two <= the value, we will perform 3 Newton interations to
 * arive at our guess.
 *
 * We have verified that fsqrt(x) == (FULL)sqrt((double)x), or
 * fsqrt(x)-1 == (FULL)sqrt((double)x) for all 0 <= x < 2^32.
 *
 * given:
 *	x		compute the integer square root of x
 */
static FULL
fsqrt(FULL x)
{
	FULL y;		/* (FULL)temporary value */
	int i;

	/* firewall - deal with 0 */
	if (x == 0) {
	    return 0;
	}

	/* ignore Saber-C warning #530 about empty for statement */
	/*	  ok to ignore in proc fsqrt */
	/* determine our initial guess */
	for (i=0, y=x; y >= (FULL)256; i+=8, y>>=8) {
	}
	y = isqrt_pow2[i + topbit[y]];

	/* perform 3 Newton interations */
	y = (y+x/y)>>1;
	y = (y+x/y)>>1;
	y = (y+x/y)>>1;
#if FULL_BITS == 64
	y = (y+x/y)>>1;
#endif

	/* return the result */
	return y;
}