zmod.c

早期freebsd实现

/*
 * Copyright (c) 1994 David I. Bell
 * Permission is granted to use, distribute, or modify this source,
 * provided that this copyright notice remains intact.
 *
 * Routines to do modulo arithmetic both normally and also using the REDC
 * algorithm given by Peter L. Montgomery in Mathematics of Computation,
 * volume 44, number 170 (April, 1985).  For multiple multiplies using
 * the same large modulus, the REDC algorithm avoids the usual division
 * by the modulus, instead replacing it with two multiplies or else a
 * special algorithm.  When these two multiplies or the special algorithm
 * are faster then the division, then the REDC algorithm is better.
 */

#include "zmath.h"


#define	POWBITS	4		/* bits for power chunks (must divide BASEB) */
#define	POWNUMS	(1<<POWBITS)	/* number of powers needed in table */


LEN _pow2_ = POW_ALG2;		/* modulo size to use REDC for powers */
LEN _redc2_ = REDC_ALG2;	/* modulo size to use second REDC algorithm */

static REDC *powermodredc = NULL;	/* REDC info for raising to power */

#if 0
extern void zaddmod MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res));
extern void znegmod MATH_PROTO((ZVALUE z1, ZVALUE z2, ZVALUE *res));

/*
 * Multiply two numbers together and then mod the result with a third number.
 * The two numbers to be multiplied can be negative or out of modulo range.
 * The result will be in the range 0 to the modulus - 1.
 */
void
zmulmod(z1, z2, z3, res)
	ZVALUE z1;		/* first number to be multiplied */
	ZVALUE z2;		/* second number to be multiplied */
	ZVALUE z3;		/* number to take mod with */
	ZVALUE *res;		/* result */
{
	ZVALUE tmp;
	FULL prod;
	FULL digit;
	BOOL neg;

	if (ziszero(z3) || zisneg(z3))
		math_error("Mod of non-positive integer");
	if (ziszero(z1) || ziszero(z2) || zisunit(z3)) {
		*res = _zero_;
		return;
	}

	/*
	 * If the modulus is a single digit number, then do the result
	 * cheaply.  Check especially for a small power of two.
	 */
	if (zistiny(z3)) {
		neg = (z1.sign != z2.sign);
		digit = z3.v[0];
		if ((digit & -digit) == digit) {	/* NEEDS 2'S COMP */
			prod = ((FULL) z1.v[0]) * ((FULL) z2.v[0]);
			prod &= (digit - 1);
		} else {
			z1.sign = 0;
			z2.sign = 0;
			prod = (FULL) zmodi(z1, (long) digit);
			prod *= (FULL) zmodi(z2, (long) digit);
			prod %= digit;
		}
		if (neg && prod)
			prod = digit - prod;
		itoz((long) prod, res);
		return;
	}

	/*
	 * The modulus is more than one digit.
	 * Actually do the multiply and divide if necessary.
	 */
	zmul(z1, z2, &tmp);
	if (zispos(tmp) && ((tmp.len < z3.len) || ((tmp.len == z3.len) &&
		(tmp.v[tmp.len-1] < z2.v[z3.len-1]))))
	{
		*res = tmp;
		return;
	}
	zmod(tmp, z3, res);
	zfree(tmp);
}


/*
 * Square a number and then mod the result with a second number.
 * The number to be squared can be negative or out of modulo range.
 * The result will be in the range 0 to the modulus - 1.
 */
void
zsquaremod(z1, z2, res)
	ZVALUE z1;		/* number to be squared */
	ZVALUE z2;		/* number to take mod with */
	ZVALUE *res;		/* result */
{
	ZVALUE tmp;
	FULL prod;
	FULL digit;

	if (ziszero(z2) || zisneg(z2))
		math_error("Mod of non-positive integer");
	if (ziszero(z1) || zisunit(z2)) {
		*res = _zero_;
		return;
	}

	/*
	 * If the modulus is a single digit number, then do the result
	 * cheaply.  Check especially for a small power of two.
	 */
	if (zistiny(z2)) {
		digit = z2.v[0];
		if ((digit & -digit) == digit) {	/* NEEDS 2'S COMP */
			prod = (FULL) z1.v[0];
			prod = (prod * prod) & (digit - 1);
		} else {
			z1.sign = 0;
			prod = (FULL) zmodi(z1, (long) digit);
			prod = (prod * prod) % digit;
		}
		itoz((long) prod, res);
		return;
	}

	/*
	 * The modulus is more than one digit.
	 * Actually do the square and divide if necessary.
	 */
	zsquare(z1, &tmp);
	if ((tmp.len < z2.len) ||
		((tmp.len == z2.len) && (tmp.v[tmp.len-1] < z2.v[z2.len-1]))) {
			*res = tmp;
			return;
	}
	zmod(tmp, z2, res);
	zfree(tmp);
}


/*
 * Add two numbers together and then mod the result with a third number.
 * The two numbers to be added can be negative or out of modulo range.
 * The result will be in the range 0 to the modulus - 1.
 */
static void
zaddmod(z1, z2, z3, res)
	ZVALUE z1;		/* first number to be added */
	ZVALUE z2;		/* second number to be added */
	ZVALUE z3;		/* number to take mod with */
	ZVALUE *res;		/* result */
{
	ZVALUE tmp;
	FULL sumdigit;
	FULL moddigit;

	if (ziszero(z3) || zisneg(z3))
		math_error("Mod of non-positive integer");
	if ((ziszero(z1) && ziszero(z2)) || zisunit(z3)) {
		*res = _zero_;
		return;
	}
	if (zistwo(z2)) {
		if ((z1.v[0] + z2.v[0]) & 0x1)
			*res = _one_;
		else
			*res = _zero_;
		return;
	}
	zadd(z1, z2, &tmp);
	if (zisneg(tmp) || (tmp.len > z3.len)) {
		zmod(tmp, z3, res);
		zfree(tmp);
		return;
	}
	sumdigit = tmp.v[tmp.len - 1];
	moddigit = z3.v[z3.len - 1];
	if ((tmp.len < z3.len) || (sumdigit < moddigit)) {
		*res = tmp;
		return;
	}
	if (sumdigit < 2 * moddigit) {
		zsub(tmp, z3, res);
		zfree(tmp);
		return;
	}
	zmod(tmp, z2, res);
	zfree(tmp);
}


/*
 * Subtract two numbers together and then mod the result with a third number.
 * The two numbers to be subtract can be negative or out of modulo range.
 * The result will be in the range 0 to the modulus - 1.
 */
void
zsubmod(z1, z2, z3, res)
	ZVALUE z1;		/* number to be subtracted from */
	ZVALUE z2;		/* number to be subtracted */
	ZVALUE z3;		/* number to take mod with */
	ZVALUE *res;		/* result */
{
	if (ziszero(z3) || zisneg(z3))
		math_error("Mod of non-positive integer");
	if (ziszero(z2)) {
		zmod(z1, z3, res);
		return;
	}
	if (ziszero(z1)) {
		znegmod(z2, z3, res);
		return;
	}
	if ((z1.sign == z2.sign) && (z1.len == z2.len) &&
		(z1.v[0] == z2.v[0]) && (zcmp(z1, z2) == 0)) {
			*res = _zero_;
			return;
	}
	z2.sign = !z2.sign;
	zaddmod(z1, z2, z3, res);
}


/*
 * Calculate the negative of a number modulo another number.
 * The number to be negated can be negative or out of modulo range.
 * The result will be in the range 0 to the modulus - 1.
 */
static void
znegmod(z1, z2, res)
	ZVALUE z1;		/* number to take negative of */
	ZVALUE z2;		/* number to take mod with */
	ZVALUE *res;		/* result */
{
	int sign;
	int cv;

	if (ziszero(z2) || zisneg(z2))
		math_error("Mod of non-positive integer");
	if (ziszero(z1) || zisunit(z2)) {
		*res = _zero_;
		return;
	}
	if (zistwo(z2)) {
		if (z1.v[0] & 0x1)
			*res = _one_;
		else
			*res = _zero_;
		return;
	}

	/*
	 * If the absolute value of the number is within the modulo range,
	 * then the result is just a copy or a subtraction.  Otherwise go
	 * ahead and negate and reduce the result.
	 */
	sign = z1.sign;
	z1.sign = 0;
	cv = zrel(z1, z2);
	if (cv == 0) {
		*res = _zero_;
		return;
	}
	if (cv < 0) {
		if (sign)
			zcopy(z1, res);
		else
			zsub(z2, z1, res);
		return;
	}
	z1.sign = !sign;
	zmod(z1, z2, res);
}
#endif


/*
 * Calculate the number congruent to the given number whose absolute
 * value is minimal.  The number to be reduced can be negative or out of
 * modulo range.  The result will be within the range -int((modulus-1)/2)
 * to int(modulus/2) inclusive.  For example, for modulus 7, numbers are
 * reduced to the range [-3, 3], and for modulus 8, numbers are reduced to
 * the range [-3, 4].
 */
void
zminmod(z1, z2, res)
	ZVALUE z1;		/* number to find minimum congruence of */
	ZVALUE z2;		/* number to take mod with */
	ZVALUE *res;		/* result */
{
	ZVALUE tmp1, tmp2;
	int sign;
	int cv;

	if (ziszero(z2) || zisneg(z2))
		math_error("Mod of non-positive integer");
	if (ziszero(z1) || zisunit(z2)) {
		*res = _zero_;
		return;
	}
	if (zistwo(z2)) {
		if (zisodd(z1))
			*res = _one_;
		else
			*res = _zero_;
		return;
	}

	/*
	 * Do a quick check to see if the number is very small compared
	 * to the modulus.  If so, then the result is obvious.
	 */
	if (z1.len < z2.len - 1) {
		zcopy(z1, res);
		return;
	}

	/*
	 * Now make sure the input number is within the modulo range.
	 * If not, then reduce it to be within range and make the
	 * quick check again.
	 */
	sign = z1.sign;
	z1.sign = 0;
	cv = zrel(z1, z2);
	if (cv == 0) {
		*res = _zero_;
		return;
	}
	tmp1 = z1;
	if (cv > 0) {
		z1.sign = (BOOL)sign;
		zmod(z1, z2, &tmp1);
		if (tmp1.len < z2.len - 1) {
			*res = tmp1;
			return;
		}
		sign = 0;
	}

	/*
	 * Now calculate the difference of the modulus and the absolute
	 * value of the original number.  Compare the original number with
	 * the difference, and return the one with the smallest absolute
	 * value, with the correct sign.  If the two values are equal, then
	 * return the positive result.
	 */
	zsub(z2, tmp1, &tmp2);
	cv = zrel(tmp1, tmp2);
	if (cv < 0) {
		zfree(tmp2);
		tmp1.sign = (BOOL)sign;
		if (tmp1.v == z1.v)
			zcopy(tmp1, res);
		else
			*res = tmp1;
	} else {
		if (cv)
			tmp2.sign = !sign;
		if (tmp1.v != z1.v)
			zfree(tmp1);
		*res = tmp2;
	}
}


/*
 * Compare two numbers for equality modulo a third number.
 * The two numbers to be compared can be negative or out of modulo range.
 * Returns TRUE if the numbers are not congruent, and FALSE if they are
 * congruent.
 */
BOOL
zcmpmod(z1, z2, z3)
	ZVALUE z1;		/* first number to be compared */
	ZVALUE z2;		/* second number to be compared */
	ZVALUE z3;		/* modulus */
{
	ZVALUE tmp1, tmp2, tmp3;
	FULL digit;
	LEN len;
	int cv;

	if (zisneg(z3) || ziszero(z3))
		math_error("Non-positive modulus in zcmpmod");
	if (zistwo(z3))
		return (((z1.v[0] + z2.v[0]) & 0x1) != 0);

	/*
	 * If the two numbers are equal, then their mods are equal.
	 */
	if ((z1.sign == z2.sign) && (z1.len == z2.len) &&
		(z1.v[0] == z2.v[0]) && (zcmp(z1, z2) == 0))
			return FALSE;

	/*
	 * If both numbers are negative, then we can make them positive.
	 */
	if (zisneg(z1) && zisneg(z2)) {
		z1.sign = 0;
		z2.sign = 0;
	}

	/*
	 * For small negative numbers, make them positive before comparing.
	 * In any case, the resulting numbers are in tmp1 and tmp2.
	 */
	tmp1 = z1;
	tmp2 = z2;
	len = z3.len;
	digit = z3.v[len - 1];

	if (zisneg(z1) && ((z1.len < len) ||
		((z1.len == len) && (z1.v[z1.len - 1] < digit))))
			zadd(z1, z3, &tmp1);

	if (zisneg(z2) && ((z2.len < len) ||
		((z2.len == len) && (z2.v[z2.len - 1] < digit))))
			zadd(z2, z3, &tmp2);

	/*
	 * Now compare the two numbers for equality.
	 * If they are equal we are all done.
	 */
	if (zcmp(tmp1, tmp2) == 0) {
		if (tmp1.v != z1.v)
			zfree(tmp1);
		if (tmp2.v != z2.v)
			zfree(tmp2);
		return FALSE;
	}

	/*
	 * They are not identical.  Now if both numbers are positive
	 * and less than the modulus, then they are definitely not equal.
	 */
	if ((tmp1.sign == tmp2.sign) &&
		((tmp1.len < len) || (zrel(tmp1, z3) < 0)) &&
		((tmp2.len < len) || (zrel(tmp2, z3) < 0)))
	{