zmod.c

早期freebsd实现

		if (tmp1.v != z1.v)
			zfree(tmp1);
		if (tmp2.v != z2.v)
			zfree(tmp2);
		return TRUE;
	}

	/*
	 * Either one of the numbers is negative or is large.
	 * So do the standard thing and subtract the two numbers.
	 * Then they are equal if the result is 0 (mod z3).
	 */
	zsub(tmp1, tmp2, &tmp3);
	if (tmp1.v != z1.v)
		zfree(tmp1);
	if (tmp2.v != z2.v)
		zfree(tmp2);

	/*
	 * Compare the result with the modulus to see if it is equal to
	 * or less than the modulus.  If so, we know the mod result.
	 */
	tmp3.sign = 0;
	cv = zrel(tmp3, z3);
	if (cv == 0) {
		zfree(tmp3);
		return FALSE;
	}
	if (cv < 0) {
		zfree(tmp3);
		return TRUE;
	}

	/*
	 * We are forced to actually do the division.
	 * The numbers are congruent if the result is zero.
	 */
	zmod(tmp3, z3, &tmp1);
	zfree(tmp3);
	if (ziszero(tmp1)) {
		zfree(tmp1);
		return FALSE;
	} else {
		zfree(tmp1);
		return TRUE;
	}
}


/*
 * Compute the result of raising one number to a power modulo another number.
 * That is, this computes:  a^b (modulo c).
 * This calculates the result by examining the power POWBITS bits at a time,
 * using a small table of POWNUMS low powers to calculate powers for those bits,
 * and repeated squaring and multiplying by the partial powers to generate
 * the complete power.  If the power being raised to is high enough, then
 * this uses the REDC algorithm to avoid doing many divisions.  When using
 * REDC, multiple calls to this routine using the same modulus will be
 * slightly faster.
 */
void
zpowermod(z1, z2, z3, res)
	ZVALUE z1, z2, z3, *res;
{
	HALF *hp;		/* pointer to current word of the power */
	REDC *rp;		/* REDC information to be used */
	ZVALUE *pp;		/* pointer to low power table */
	ZVALUE ans, temp;	/* calculation values */
	ZVALUE modpow;		/* current small power */
	ZVALUE lowpowers[POWNUMS];	/* low powers */
	int sign;		/* original sign of number */
	int curshift;		/* shift value for word of power */
	HALF curhalf;		/* current word of power */
	unsigned int curpow;	/* current low power */
	unsigned int curbit;	/* current bit of low power */
	int i;

	if (zisneg(z3) || ziszero(z3))
		math_error("Non-positive modulus in zpowermod");
	if (zisneg(z2))
		math_error("Negative power in zpowermod");

	sign = z1.sign;
	z1.sign = 0;

	/*
	 * Check easy cases first.
	 */
	if ((ziszero(z1) && !ziszero(z2)) || zisunit(z3)) {
		/* 0^(non_zero) or x^y mod 1 always produces zero */
		*res = _zero_;
		return;
	}
	if (ziszero(z2)) {			/* x^0 == 1 */
		*res = _one_;
		return;
	}
	if (zistwo(z3)) {			/* mod 2 */
		if (zisodd(z1))
			*res = _one_;
		else
			*res = _zero_;
		return;
	}
	if (zisunit(z1) && (!sign || ziseven(z2))) {
		/* 1^x or (-1)^(2x) */
		*res = _one_;
		return;
	}

	/*
	 * Normalize the number being raised to be non-negative and to lie
	 * within the modulo range.  Then check for zero or one specially.
	 */
	zmod(z1, z3, &temp);
	if (ziszero(temp)) {
		zfree(temp);
		*res = _zero_;
		return;
	}
	z1 = temp;
	if (sign) {
		zsub(z3, z1, &temp);
		zfree(z1);
		z1 = temp;
	}
	if (zisunit(z1)) {
		zfree(z1);
		*res = _one_;
		return;
	}

	/*
	 * If the modulus is odd, large enough, is not one less than an
	 * exact power of two, and if the power is large enough, then use
	 * the REDC algorithm.  The size where this is done is configurable.
	 */
	if ((z2.len > 1) && (z3.len >= _pow2_) && zisodd(z3)
		&& !zisallbits(z3))
	{
		if (powermodredc && zcmp(powermodredc->mod, z3)) {
			zredcfree(powermodredc);
			powermodredc = NULL;
		}
		if (powermodredc == NULL)
			powermodredc = zredcalloc(z3);
		rp = powermodredc;
		zredcencode(rp, z1, &temp);
		zredcpower(rp, temp, z2, &z1);
		zfree(temp);
		zredcdecode(rp, z1, res);
		zfree(z1);
		return;
	}

	/*
	 * Modulus or power is small enough to perform the power raising
	 * directly.  Initialize the table of powers.
	 */
	for (pp = &lowpowers[2]; pp < &lowpowers[POWNUMS]; pp++)
		pp->len = 0;
	lowpowers[0] = _one_;
	lowpowers[1] = z1;
	ans = _one_;

	hp = &z2.v[z2.len - 1];
	curhalf = *hp;
	curshift = BASEB - POWBITS;
	while (curshift && ((curhalf >> curshift) == 0))
		curshift -= POWBITS;

	/*
	 * Calculate the result by examining the power POWBITS bits at a time,
	 * and use the table of low powers at each iteration.
	 */
	for (;;) {
		curpow = (curhalf >> curshift) & (POWNUMS - 1);
		pp = &lowpowers[curpow];

		/*
		 * If the small power is not yet saved in the table, then
		 * calculate it and remember it in the table for future use.
		 */
		if (pp->len == 0) {
			if (curpow & 0x1)
				zcopy(z1, &modpow);
			else
				modpow = _one_;

			for (curbit = 0x2; curbit <= curpow; curbit *= 2) {
				pp = &lowpowers[curbit];
				if (pp->len == 0) {
					zsquare(lowpowers[curbit/2], &temp);
					zmod(temp, z3, pp);
					zfree(temp);
				}
				if (curbit & curpow) {
					zmul(*pp, modpow, &temp);
					zfree(modpow);
					zmod(temp, z3, &modpow);
					zfree(temp);
				}
			}
			pp = &lowpowers[curpow];
			*pp = modpow;
		}

		/*
		 * If the power is nonzero, then accumulate the small power
		 * into the result.
		 */
		if (curpow) {
			zmul(ans, *pp, &temp);
			zfree(ans);
			zmod(temp, z3, &ans);
			zfree(temp);
		}

		/*
		 * Select the next POWBITS bits of the power, if there is
		 * any more to generate.
		 */
		curshift -= POWBITS;
		if (curshift < 0) {
			if (hp-- == z2.v)
				break;
			curhalf = *hp;
			curshift = BASEB - POWBITS;
		}

		/*
		 * Square the result POWBITS times to make room for the next
		 * chunk of bits.
		 */
		for (i = 0; i < POWBITS; i++) {
			zsquare(ans, &temp);
			zfree(ans);
			zmod(temp, z3, &ans);
			zfree(temp);
		}
	}

	for (pp = &lowpowers[2]; pp < &lowpowers[POWNUMS]; pp++) {
		if (pp->len)
			freeh(pp->v);
	}
	*res = ans;
}


/*
 * Initialize the REDC algorithm for a particular modulus,
 * returning a pointer to a structure that is used for other
 * REDC calls.  An error is generated if the structure cannot
 * be allocated.  The modulus must be odd and positive.
 */
REDC *
zredcalloc(z1)
	ZVALUE z1;		/* modulus to initialize for */
{
	REDC *rp;		/* REDC information */
	ZVALUE tmp;
	long bit;

	if (ziseven(z1) || zisneg(z1))
		math_error("REDC requires positive odd modulus");

	rp = (REDC *) malloc(sizeof(REDC));
	if (rp == NULL)
		math_error("Cannot allocate REDC structure");

	/*
	 * Round up the binary modulus to the next power of two
	 * which is at a word boundary.  Then the shift and modulo
	 * operations mod the binary modulus can be done very cheaply.
	 * Calculate the REDC format for the number 1 for future use.
	 */
	bit = zhighbit(z1) + 1;
	if (bit % BASEB)
		bit += (BASEB - (bit % BASEB));
	zcopy(z1, &rp->mod);
	zbitvalue(bit, &tmp);
	z1.sign = 1;
	(void) zmodinv(z1, tmp, &rp->inv);
	zmod(tmp, rp->mod, &rp->one);
	zfree(tmp);
	rp->len = bit / BASEB;
	return rp;
}


/*
 * Free any numbers associated with the specified REDC structure,
 * and then the REDC structure itself.
 */
void
zredcfree(rp)
	REDC *rp;		/* REDC information to be cleared */
{
	zfree(rp->mod);
	zfree(rp->inv);
	zfree(rp->one);
	free(rp);
}


/*
 * Convert a normal number into the specified REDC format.
 * The number to be converted can be negative or out of modulo range.
 * The resulting number can be used for multiplying, adding, subtracting,
 * or comparing with any other such converted numbers, as if the numbers
 * were being calculated modulo the number which initialized the REDC
 * information.  When the final value is unconverted, the result is the
 * same as if the usual operations were done with the original numbers.
 */
void
zredcencode(rp, z1, res)
	REDC *rp;		/* REDC information */
	ZVALUE z1;		/* number to be converted */
	ZVALUE *res;		/* returned converted number */
{
	ZVALUE tmp1, tmp2;

	/*
	 * Handle the cases 0, 1, -1, and 2 specially since these are
	 * easy to calculate.  Zero transforms to zero, and the others
	 * can be obtained from the precomputed REDC format for 1 since
	 * addition and subtraction act normally for REDC format numbers.
	 */
	if (ziszero(z1)) {
		*res = _zero_;
		return;
	}
	if (zisone(z1)) {
		zcopy(rp->one, res);
		return;
	}
	if (zisunit(z1)) {
		zsub(rp->mod, rp->one, res);
		return;
	}
	if (zistwo(z1)) {
		zadd(rp->one, rp->one, &tmp1);
		if (zrel(tmp1, rp->mod) < 0) {
			*res = tmp1;
			return;
		}
		zsub(tmp1, rp->mod, res);
		zfree(tmp1);
		return;
	}

	/*
	 * Not a trivial number to convert, so do the full transformation.
	 * Convert negative numbers to positive numbers before converting.
	 */
	tmp1.len = 0;
	if (zisneg(z1)) {
		zmod(z1, rp->mod, &tmp1);
		z1 = tmp1;
	}
	zshift(z1, rp->len * BASEB, &tmp2);
	if (tmp1.len)
		zfree(tmp1);
	zmod(tmp2, rp->mod, res);
	zfree(tmp2);
}


/*
 * The REDC algorithm used to convert numbers out of REDC format and also
 * used after multiplication of two REDC numbers.  Using this routine
 * avoids any divides, replacing the divide by two multiplications.
 * If the numbers are very large, then these two multiplies will be
 * quicker than the divide, since dividing is harder than multiplying.
 */
void
zredcdecode(rp, z1, res)
	REDC *rp;		/* REDC information */
	ZVALUE z1;		/* number to be transformed */
	ZVALUE *res;		/* returned transformed number */
{
	ZVALUE tmp1, tmp2;	/* temporaries */
	HALF *hp;		/* saved pointer to tmp2 value */

	if (zisneg(z1))
		math_error("Negative number for zredc");

	/*
	 * Check first for the special values for 0 and 1 that are easy.
	 */
	if (ziszero(z1)) {
		*res = _zero_;
		return;
	}
	if ((z1.len == rp->one.len) && (z1.v[0] == rp->one.v[0]) &&
		(zcmp(z1, rp->one) == 0)) {
			*res = _one_;
			return;
	}

	/*
	 * First calculate the following:
	 * 	tmp2 = ((z1 % 2^bitnum) * inv) % 2^bitnum.
	 * The mod operations can be done with no work since the bit
	 * number was selected as a multiple of the word size.  Just
	 * reduce the sizes of the numbers as required.
	 */
	tmp1 = z1;
	if (tmp1.len > rp->len)
		tmp1.len = rp->len;
	zmul(tmp1, rp->inv, &tmp2);
	if (tmp2.len > rp->len)
		tmp2.len = rp->len;

	/*
	 * Next calculate the following:
	 *	res = (z1 + tmp2 * modulus) / 2^bitnum
	 * The division by a power of 2 is always exact, and requires no
	 * work.  Just adjust the address and length of the number to do
	 * the divide, but save the original pointer for freeing later.
	 */
	zmul(tmp2, rp->mod, &tmp1);
	zfree(tmp2);
	zadd(z1, tmp1, &tmp2);
	zfree(tmp1);
	hp = tmp2.v;
	if (tmp2.len <= rp->len) {
		freeh(hp);
		*res = _zero_;
		return;
	}
	tmp2.v += rp->len;
	tmp2.len -= rp->len;

	/*
	 * Finally do a final modulo by a simple subtraction if necessary.
	 * This is all that is needed because the previous calculation is
	 * guaranteed to always be less than twice the modulus.
	 */
	if (zrel(tmp2, rp->mod) < 0)
		zcopy(tmp2, res);
	else
		zsub(tmp2, rp->mod, res);
	freeh(hp);
}