zmod.c

Calc Software Package for Number Calc

		if (ztmp.len)
			zfree(ztmp);
		*res = _zero_;
		return;
	}
	if (zisone(z1)) {
		if (ztmp.len)
			zfree(ztmp);
		*res = _one_;
		return;
	}

	/*
	 * If modulus is large enough use zmod5
	 */
	if (z3.len >= conf->pow2) {
		if (havelastmod && zcmp(z3, *lastmod)) {
			zfree(*lastmod);
			zfree(*lastmodinv);
			havelastmod = FALSE;
		}
		if (!havelastmod) {
			zcopy(z3, lastmod);
			zbitvalue(2 * z3.len * BASEB, &temp);
			zquo(temp, z3, lastmodinv, 0);
			zfree(temp);
			havelastmod = TRUE;
		}

		/* zzz */
		for (pp = &lowpowers[2]; pp <= &lowpowers[POWNUMS-1]; pp++) {
			pp->len = 0;
			pp->v = NULL;
		}
		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->v == NULL) {
				if (curpow & 0x1)
					zcopy(z1, &modpow);
				else
					modpow = _one_;

				for (curbit = 0x2;
				     curbit <= curpow;
				     curbit *= 2) {
					pp = &lowpowers[curbit];
					if (pp->v == NULL) {
						zsquare(lowpowers[curbit/2],
							&temp);
						zmod5(&temp);
						zcopy(temp, pp);
						zfree(temp);
					}
					if (curbit & curpow) {
						zmul(*pp, modpow, &temp);
						zfree(modpow);
						zmod5(&temp);
						zcopy(temp, &modpow);
						zfree(temp);
					}
				}
				pp = &lowpowers[curpow];
				if (pp->v != NULL) {
					zfree(*pp);
				}
				*pp = modpow;
			}

			/*
			 * If the power is nonzero, then accumulate the small
			 * power into the result.
			 */
			if (curpow) {
				zmul(ans, *pp, &temp);
				zfree(ans);
				zmod5(&temp);
				zcopy(temp, &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);
				zmod5(&temp);
				zcopy(temp, &ans);
				zfree(temp);
			}
		}

		for (pp = &lowpowers[2]; pp <= &lowpowers[POWNUMS-1]; pp++) {
			if (pp->v != NULL)
				freeh(pp->v);
		}
		*res = ans;
		if (ztmp.len)
			zfree(ztmp);
		return;
	}

	/*
	 * If the modulus is odd and small enough then use
	 * the REDC algorithm.	The size where this is done is configurable.
	 */
	if (z3.len < conf->redc2 && zisodd(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-1]; pp++) {
		pp->len = 0;
		pp->v = NULL;
	}
	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->v == NULL) {
			if (curpow & 0x1)
				zcopy(z1, &modpow);
			else
				modpow = _one_;

			for (curbit = 0x2; curbit <= curpow; curbit *= 2) {
				pp = &lowpowers[curbit];
				if (pp->v == NULL) {
					zsquare(lowpowers[curbit/2], &temp);
					zmod(temp, z3, pp, 0);
					zfree(temp);
				}
				if (curbit & curpow) {
					zmul(*pp, modpow, &temp);
					zfree(modpow);
					zmod(temp, z3, &modpow, 0);
					zfree(temp);
				}
			}
			pp = &lowpowers[curpow];
			if (pp->v != NULL) {
				zfree(*pp);
			}
			*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, 0);
			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, 0);
			zfree(temp);
		}
	}

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

/*
 * Given a positive odd N-word integer z, evaluate minv(-z, BASEB^N)
 */
static void
zredcmodinv(ZVALUE z, ZVALUE *res)
{
	ZVALUE tmp;
	HALF *a0, *a, *b;
	HALF bit, h, inv, v;
	FULL f;
	LEN N, i, j, len;

	N = z.len;
	tmp.sign = 0;
	tmp.len = N;
	tmp.v = alloc(N);
	zclearval(tmp);
	*tmp.v = 1;
	h = 1 + *z.v;
	bit = 1;
	inv = 1;
	while (h) {
		bit <<= 1;
		if (bit & h) {
			inv |= bit;
			h += bit * *z.v;
		}
	}

	j = N;
	a0 = tmp.v;
	while (j-- > 0) {
		v = inv * *a0;
		i = j;
		a = a0;
		b = z.v;
		f = (FULL) v * (FULL) *b++ + (FULL) *a++;
		*a0 = v;
		while (i-- > 0) {
			f = (FULL) v * (FULL) *b++  + (FULL) *a + (f >> BASEB);
			*a++ = (HALF) f;
		}
		while (j > 0 && *++a0 == 0)
			j--;
	}
	a = tmp.v + N;
	len = N;
	while (*--a == 0)
		len--;
	tmp.len = len;
	zcopy(tmp, res);
	zfree(tmp);
}


/*
 * 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.
 *
 * given:
 *	z1		modulus to initialize for
 */
REDC *
zredcalloc(ZVALUE z1)
{
	REDC *rp;		/* REDC information */
	ZVALUE tmp;
	long bit;

	if (ziseven(z1) || zisneg(z1)) {
		math_error("REDC requires positive odd modulus");
		/*NOTREACHED*/
	}

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

	/*
	 * 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.
	 */
	zcopy(z1, &rp->mod);
	zredcmodinv(z1, &rp->inv);
	bit = zhighbit(z1) + 1;
	if (bit % BASEB)
		bit += (BASEB - (bit % BASEB));
	zbitvalue(bit, &tmp);
	zmod(tmp, rp->mod, &rp->one, 0);
	zfree(tmp);
	rp->len = (LEN)(bit / BASEB);
	return rp;
}


/*
 * Free any numbers associated with the specified REDC structure,
 * and then the REDC structure itself.
 *
 * given:
 *	rp		REDC information to be cleared
 */
void
zredcfree(REDC *rp)
{
	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.
 *
 * given:
 *	rp		REDC information
 *	z1		number to be converted
 *	res		returned converted number
 */
void
zredcencode(REDC *rp, ZVALUE z1, ZVALUE *res)
{
	ZVALUE tmp1;

	/*
	 * Confirm or initialize lastmod information when modulus is a
	 * big number.
	 */

	if (rp->len >= conf->pow2) {
		if (havelastmod && zcmp(rp->mod, *lastmod)) {
			zfree(*lastmod);
			zfree(*lastmodinv);
			havelastmod = FALSE;
		}
		if (!havelastmod) {
			zcopy(rp->mod, lastmod);
			zbitvalue(2 * rp->len * BASEB, &tmp1);
			zquo(tmp1, rp->mod, lastmodinv, 0);
			zfree(tmp1);
			havelastmod = TRUE;
		}
	}
	/*
	 * 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.
	 */
	zshift(z1, rp->len * BASEB, &tmp1);
	if (rp->len < conf->pow2)
		zmod(tmp1, rp->mod, res, 0);
	else
		zmod6(tmp1, res);
	zfree(tmp1);
}


/*
 * 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.
 *
 * given:
 *	rp		REDC information
 *	z1		number to be transformed
 *	res		returned transformed number
 */
void
zredcdecode(REDC *rp, ZVALUE z1, ZVALUE *res)
{
	ZVALUE tmp1, tmp2;
	ZVALUE ztmp;
	ZVALUE ztop;
	ZVALUE zp1;
	FULL muln;
	HALF *h1;
	HALF *h3;
	HALF *hd = NULL;
	HALF Ninv;
	LEN modlen;
	LEN len;
	FULL f;
	int sign;
	int i, j;

	/*
	 * 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;
	}
	ztop.len = 0;
	ztmp.len = 0;
	modlen = rp->len;
	sign = z1.sign;
	z1.sign = 0;
	if (z1.len > modlen) {
		ztop.v = z1.v + modlen;
		ztop.len = z1.len - modlen;
		ztop.sign = 0;
		if (zrel(ztop, rp->mod) >= 0) {
			zmod(ztop, rp->mod, &ztmp, 0);
			ztop = ztmp;
		}
		len = modlen;
		h1 = z1.v + len;
		while (len > 0 && *--h1 == 0)
			len--;
		if (len == 0) {
			if (ztmp.len)
				*res = ztmp;
			else
				zcopy(ztop, res);
			return;
		}
		z1.len = len;

	}
	if (rp->mod.len < conf->pow2) {
		Ninv = rp->inv.v[0];
		res->sign = 0;
		res->len = modlen;
		res->v = alloc(modlen);
		zclearval(*res);
		h1 = z1.v;
		for (i = 0; i < modlen; i++) {
			h3 = rp->mod.v;
			hd = res->v;
			f = (FULL) *hd++;
			if (i < z1.len)
				f += (FULL) *h1++;
			muln = (HALF) ((f & BASE1) * Ninv);
			f = ((muln * (FULL) *h3++) + f) >> BASEB;
			j = modlen;
			while (--j > 0) {
				f += (muln * (FULL) *h3++) + (FULL) *hd;