zmod.c

Calc Software Package for Number Calc

				hd[-1] = (HALF) f;
				f >>= BASEB;
				hd++;
			}
			hd[-1] = (HALF) f;
		}
		len = modlen;
		while (*--hd == 0 && len > 1)
			len--;
		if (len == 0)
			len = 1;
		res->len = len;
	} else {
		/* Here 0 < z1 < 2^bitnum */

		/*
		 * First calculate the following:
		 *	tmp2 = ((z1 * 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.
		 */
		zmul(z1, rp->inv, &tmp2);
		if (tmp2.len > modlen) {
			h1 = tmp2.v + modlen;
			len = modlen;
			while (len > 0 && *--h1 == 0)
				len--;
			tmp2.len = len;
		}

		/*
		 * Next calculate the following:
		 *	res = (z1 + tmp2 * modulus) / 2^bitnum
		 * Since 0 < z1 < 2^bitnum and the division is always exact,
		 * the quotient can be evaluated by rounding up
		 * (tmp2 * modulus)/2^bitnum.  This can be achieved by defining
		 * zp1 by an appropriate shift and then adding one.
		 */
		zmul(tmp2, rp->mod, &tmp1);
		zfree(tmp2);
		if (tmp1.len > modlen) {
			zp1.v = tmp1.v + modlen;
			zp1.len = tmp1.len - modlen;
			zp1.sign = 0;
			zadd(zp1, _one_, res);
		} else {
			*res = _one_;
		}
		zfree(tmp1);
	}
	if (ztop.len) {
		zadd(*res, ztop, &tmp1);
		zfree(*res);
		if (ztmp.len)
			zfree(ztmp);
		*res = tmp1;
	}

	/*
	 * 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(*res, rp->mod) >= 0) {
		zsub(*res, rp->mod, &tmp1);
		zfree(*res);
		*res = tmp1;
	}
	if (sign && !ziszero(*res)) {
		zsub(rp->mod, *res, &tmp1);
		zfree(*res);
		*res = tmp1;
	}
	return;
}


/*
 * Multiply two numbers in REDC format together producing a result also
 * in REDC format.  If the result is converted back to a normal number,
 * then the result is the same as the modulo'd multiplication of the
 * original numbers before they were converted to REDC format.	This
 * calculation is done in one of two ways, depending on the size of the
 * modulus.  For large numbers, the REDC definition is used directly
 * which involves three multiplies overall.  For small numbers, a
 * complicated routine is used which does the indicated multiplication
 * and the REDC algorithm at the same time to produce the result.
 *
 * given:
 *	rp		REDC information
 *	z1		first REDC number to be multiplied
 *	z2		second REDC number to be multiplied
 *	res		resulting REDC number
 */
void
zredcmul(REDC *rp, ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
	FULL mulb;
	FULL muln;
	HALF *h1;
	HALF *h2;
	HALF *h3;
	HALF *hd;
	HALF Ninv;
	HALF topdigit = 0;
	LEN modlen;
	LEN len;
	LEN len2;
	SIUNION sival1;
	SIUNION sival2;
	SIUNION carry;
	ZVALUE tmp;
	ZVALUE z1tmp, z2tmp;
	int sign;

	sign = z1.sign ^ z2.sign;
	z1.sign = 0;
	z2.sign = 0;
	z1tmp.len = 0;
	if (zrel(z1, rp->mod) >= 0) {
		zmod(z1, rp->mod, &z1tmp, 0);
		z1 = z1tmp;
	}
	z2tmp.len = 0;
	if (zrel(z2, rp->mod) >= 0) {
		zmod(z2, rp->mod, &z2tmp, 0);
		z2 = z2tmp;
	}


	/*
	 * Check for special values which we easily know the answer.
	 */
	if (ziszero(z1) || ziszero(z2)) {
		*res = _zero_;
		if (z1tmp.len)
			zfree(z1tmp);
		if (z2tmp.len)
			zfree(z2tmp);
		return;
	}

	if ((z1.len == rp->one.len) && (z1.v[0] == rp->one.v[0]) &&
		(zcmp(z1, rp->one) == 0)) {
			if (sign)
				zsub(rp->mod, z2, res);
			else
				zcopy(z2, res);
			if (z1tmp.len)
				zfree(z1tmp);
			if (z2tmp.len)
				zfree(z2tmp);
			return;
	}

	if ((z2.len == rp->one.len) && (z2.v[0] == rp->one.v[0]) &&
		(zcmp(z2, rp->one) == 0)) {
			if (sign)
				zsub(rp->mod, z1, res);
			else
				zcopy(z1, res);
			if (z1tmp.len)
				zfree(z1tmp);
			if (z2tmp.len)
				zfree(z2tmp);
			return;
	}

	/*
	 * If the size of the modulus is large, then just do the multiply,
	 * followed by the two multiplies contained in the REDC routine.
	 * This will be quicker than directly doing the REDC calculation
	 * because of the O(N^1.585) speed of the multiplies.  The size
	 * of the number which this is done is configurable.
	 */
	if (rp->mod.len >= conf->redc2) {
		zmul(z1, z2, &tmp);
		zredcdecode(rp, tmp, res);
		zfree(tmp);
		if (sign && !ziszero(*res)) {
			zsub(rp->mod, *res, &tmp);
			zfree(*res);
			*res = tmp;
		}
		if (z1tmp.len)
			zfree(z1tmp);
		if (z2tmp.len)
			zfree(z2tmp);
		return;
	}

	/*
	 * The number is small enough to calculate by doing the O(N^2) REDC
	 * algorithm directly.	This algorithm performs the multiplication and
	 * the reduction at the same time.  Notice the obscure facts that
	 * only the lowest word of the inverse value is used, and that
	 * there is no shifting of the partial products as there is in a
	 * normal multiply.
	 */
	modlen = rp->mod.len;
	Ninv = rp->inv.v[0];

	/*
	 * Allocate the result and clear it.
	 * The size of the result will be equal to or smaller than
	 * the modulus size.
	 */
	res->sign = 0;
	res->len = modlen;
	res->v = alloc(modlen);

	hd = res->v;
	len = modlen;
	zclearval(*res);

	/*
	 * Do this outermost loop over all the digits of z1.
	 */
	h1 = z1.v;
	len = z1.len;
	while (len--) {
		/*
		 * Start off with the next digit of z1, the first
		 * digit of z2, and the first digit of the modulus.
		 */
		mulb = (FULL) *h1++;
		h2 = z2.v;
		h3 = rp->mod.v;
		hd = res->v;
		sival1.ivalue = mulb * ((FULL) *h2++) + ((FULL) *hd++);
		muln = ((HALF) (sival1.silow * Ninv));
		sival2.ivalue = muln * ((FULL) *h3++) + ((FULL) sival1.silow);
		carry.ivalue = ((FULL) sival1.sihigh) + ((FULL) sival2.sihigh);

		/*
		 * Do this innermost loop for each digit of z2, except
		 * for the first digit which was just done above.
		 */
		len2 = z2.len;
		while (--len2 > 0) {
			sival1.ivalue = mulb * ((FULL) *h2++)
				+ ((FULL) *hd) + ((FULL) carry.silow);
			sival2.ivalue = muln * ((FULL) *h3++)
				+ ((FULL) sival1.silow);
			carry.ivalue = ((FULL) sival1.sihigh)
				+ ((FULL) sival2.sihigh)
				+ ((FULL) carry.sihigh);

			hd[-1] = sival2.silow;
			hd++;
		}

		/*
		 * Now continue the loop as necessary so the total number
		 * of iterations is equal to the size of the modulus.
		 * This acts as if the innermost loop was repeated for
		 * high digits of z2 that are zero.
		 */
		len2 = modlen - z2.len;
		while (len2--) {
			sival2.ivalue = muln * ((FULL) *h3++)
				+ ((FULL) *hd)
				+ ((FULL) carry.silow);
			carry.ivalue = ((FULL) sival2.sihigh)
				+ ((FULL) carry.sihigh);

			hd[-1] = sival2.silow;
			hd++;
		}

		carry.ivalue += topdigit;
		hd[-1] = carry.silow;
		topdigit = carry.sihigh;
	}

	/*
	 * Now continue the loop as necessary so the total number
	 * of iterations is equal to the size of the modulus.
	 * This acts as if the outermost loop was repeated for high
	 * digits of z1 that are zero.
	 */
	len = modlen - z1.len;
	while (len--) {
		/*
		 * Start off with the first digit of the modulus.
		 */
		h3 = rp->mod.v;
		hd = res->v;
		muln = ((HALF) (*hd * Ninv));
		sival2.ivalue = muln * ((FULL) *h3++) + (FULL) *hd++;
		carry.ivalue = ((FULL) sival2.sihigh);

		/*
		 * Do this innermost loop for each digit of the modulus,
		 * except for the first digit which was just done above.
		 */
		len2 = modlen;
		while (--len2 > 0) {
			sival2.ivalue = muln * ((FULL) *h3++)
				+ ((FULL) *hd) + ((FULL) carry.silow);
			carry.ivalue = ((FULL) sival2.sihigh)
				+ ((FULL) carry.sihigh);

			hd[-1] = sival2.silow;
			hd++;
		}
		carry.ivalue += topdigit;
		hd[-1] = carry.silow;
		topdigit = carry.sihigh;
	}

	/*
	 * Determine the true size of the result, taking the top digit of
	 * the current result into account.  The top digit is not stored in
	 * the number because it is temporary and would become zero anyway
	 * after the final subtraction is done.
	 */
	if (topdigit == 0) {
		len = modlen;
		while (*--hd == 0 && len > 1) {
			len--;
		}
		res->len = len;

	/*
	 * Compare the result with the modulus.
	 * If it is less than the modulus, then the calculation is complete.
	 */

		if (zrel(*res, rp->mod) < 0) {
			if (z1tmp.len)
				zfree(z1tmp);
			if (z2tmp.len)
				zfree(z2tmp);
			if (sign && !ziszero(*res)) {
				zsub(rp->mod, *res, &tmp);
				zfree(*res);
				*res = tmp;
			}
			return;
		}
	}

	/*
	 * Do a subtraction to reduce the result to a value less than
	 * the modulus.	 The REDC algorithm guarantees that a single subtract
	 * is all that is needed.  Ignore any borrowing from the possible
	 * highest word of the current result because that would affect
	 * only the top digit value that was not stored and would become
	 * zero anyway.
	 */
	carry.ivalue = 0;
	h1 = rp->mod.v;
	hd = res->v;
	len = modlen;
	while (len--) {
		carry.ivalue = BASE1 - ((FULL) *hd) + ((FULL) *h1++)
			+ ((FULL) carry.silow);
		*hd++ = (HALF)(BASE1 - carry.silow);
		carry.silow = carry.sihigh;
	}

	/*
	 * Now finally recompute the size of the result.
	 */
	len = modlen;
	hd = &res->v[len - 1];
	while ((*hd == 0) && (len > 1)) {
		hd--;
		len--;
	}
	res->len = len;
	if (z1tmp.len)
		zfree(z1tmp);
	if (z2tmp.len)
		zfree(z2tmp);
	if (sign && !ziszero(*res)) {
		zsub(rp->mod, *res, &tmp);
		zfree(*res);
		*res = tmp;
	}

}

/*
 * Square a number in REDC format producing a result also in REDC format.
 *
 * given:
 *	rp		REDC information
 *	z1		REDC number to be squared
 *	res		resulting REDC number
 */
void
zredcsquare(REDC *rp, ZVALUE z1, ZVALUE *res)
{
	FULL mulb;
	FULL muln;
	HALF *h1;
	HALF *h2;
	HALF *h3;
	HALF *hd = NULL;
	HALF Ninv;
	HALF topdigit = 0;
	LEN modlen;
	LEN len;
	SIUNION sival1;
	SIUNION sival2;
	SIUNION sival3;
	SIUNION carry;
	ZVALUE tmp, ztmp;
	FULL f;
	int i, j;

	ztmp.len = 0;
	z1.sign = 0;
	if (zrel(z1, rp->mod) >= 0) {
		zmod(z1, rp->mod, &ztmp, 0);
		z1 = ztmp;
	}
	if (ziszero(z1)) {
		*res = _zero_;
		if (ztmp.len)
			zfree(ztmp);
		return;
	}
	if ((z1.len == rp->one.len) && (z1.v[0] == rp->one.v[0]) &&
		(zcmp(z1, rp->one) == 0)) {
			zcopy(z1, res);
			if (ztmp.len)
				zfree(ztmp);
			return;
	}


	/*
	 * If the modulus is small enough, then call the multiply
	 * routine to produce the result.  Otherwise call the O(N^1.585)
	 * routines to get the answer.
	 */
	if (rp->mod.len >= conf->redc2
			|| 3 * z1.len < 2 * rp->mod.len) {
		zsquare(z1, &tmp);
		zredcdecode(rp, tmp, res);
		zfree(tmp);
		if (ztmp.len)
			zfree(ztmp);
		return;
	}
	modlen = rp->mod.len;
	Ninv = rp->inv.v[0];

	res->sign = 0;
	res->len = modlen;
	res->v = alloc(modlen);

	zclearval(*res);

	h1 = z1.v;

	for (i = 0; i < z1.len; i++) {
		mulb = (FULL) *h1++;
		h2 = h1;
		h3 = rp->mod.v;
		hd = res->v;
		if (i == 0) {
			sival1.ivalue = mulb * mulb;
			muln = (HALF) (sival1.silow * Ninv);
			sival2.ivalue = muln * ((FULL) *h3++)
				+ (FULL) sival1.silow;
			carry.ivalue = (FULL) sival1.sihigh
				+ (FULL) sival2.sihigh;
			hd++;
		} else {
			muln = (HALF) (*hd * Ninv);
			f = (muln * ((FULL) *h3++) + (FULL) *hd++) >> BASEB;
			j = i;
			while (--j > 0) {
				f += muln * ((FULL) *h3++) + *hd;
				hd[-1] = (HALF) f;
				f >>= BASEB;
				hd++;
			}
			carry.ivalue = f;
			sival1.ivalue = mulb * mulb + (FULL) carry.silow;
			sival2.ivalue = muln * ((FULL) *h3++)
				+ (FULL) *hd
				+ (FULL) sival1.silow;
			carry.ivalue = (FULL) sival1.sihigh
				+ (FULL) sival2.sihigh
				+ (FULL) carry.sihigh;
			hd[-1] = sival2.silow;
			hd++;
		}
		j = z1.len - i;
		while (--j > 0) {
			sival1.ivalue = mulb * ((FULL) *h2++);
			sival2.ivalue = ((FULL) sival1.silow << 1)
				+ muln * ((FULL) *h3++);
			sival3.ivalue = (FULL) sival2.silow
				+ (FULL) *hd
				+ (FULL) carry.silow;
			carry.ivalue = ((FULL) sival1.sihigh << 1)
				+ (FULL) sival2.sihigh
				+ (FULL) sival3.sihigh
				+ (FULL) carry.sihigh;
			hd[-1] = sival3.silow;
			hd++;
		}
		j = modlen - z1.len;
		while (j-- > 0) {
			sival1.ivalue = muln * ((FULL) *h3++)
				+ (FULL) *hd
				+ (FULL) carry.silow;
			carry.ivalue = (FULL) sival1.sihigh
				+ (FULL) carry.sihigh;
			hd[-1] = sival1.silow;
			hd++;
		}
		carry.ivalue += (FULL) topdigit;
		hd[-1] = carry.silow;
		topdigit = carry.sihigh;
	}
	i = modlen - z1.len;
	while (i-- > 0) {
		h3 = rp->mod.v;
		hd = res->v;
		muln = (HALF) (*hd * Ninv);
		sival1.ivalue = muln * ((FULL) *h3++) + (FULL) *hd++;
		carry.ivalue = (FULL) sival1.sihigh;
		j = modlen;
		while (--j > 0) {
			sival1.ivalue = muln * ((FULL) *h3++)
				+ (FULL) *hd
				+ (FULL) carry.silow;
			carry.ivalue = (FULL) sival1.sihigh
				+ (FULL) carry.sihigh;
			hd[-1] = sival1.silow;
			hd++;
		}