zmod.c

来自「Calc Software Package for Number Calc」· C语言 代码 · 共 2,164 行 · 第 1/4 页

		carry.ivalue += (FULL) topdigit;
		hd[-1] = carry.silow;
		topdigit = carry.sihigh;
	}
	if (topdigit == 0) {
		len = modlen;
		while (*--hd == 0 && len > 1) {
			len--;
		}
		res->len = len;
		if  (zrel(*res, rp->mod) < 0) {
			if (ztmp.len)
				zfree(ztmp);
			return;
		}
	}

	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;
	}

	len = modlen;
	hd = &res->v[len - 1];
	while ((*hd == 0) && (len > 1)) {
		hd--;
		len--;
	}
	res->len = len;
	if (ztmp.len)
		zfree(ztmp);
}


/*
 * Compute the result of raising a REDC format number to a power.
 * The result is within the range 0 to the modulus - 1.
 * 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.
 *
 * given:
 *	rp		REDC information
 *	z1		REDC number to be raised
 *	z2		normal number to raise number to
 *	res		result
 */
void
zredcpower(REDC *rp, ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
	HALF *hp;		/* pointer to current word of the power */
	ZVALUE *pp;		/* pointer to low power table */
	ZVALUE ans, temp;	/* calculation values */
	ZVALUE ztmp;
	ZVALUE modpow;		/* current small power */
	ZVALUE lowpowers[POWNUMS];	/* low powers */
	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 sign;
	int i;

	if (zisneg(z2)) {
		math_error("Negative power in zredcpower");
		/*NOTREACHED*/
	}

	if (zisunit(rp->mod)) {
		*res = _zero_;
		return;
	}

	sign = zisodd(z2) ? z1.sign : 0;
	z1.sign = 0;
	ztmp.len = 0;
	if (zrel(z1, rp->mod) >= 0) {
		zmod(z1, rp->mod, &ztmp, 0);
		z1 = ztmp;
	}
	/*
	 * Check for zero or the REDC format for one.
	 */
	if (ziszero(z1)) {
		if (ziszero(z2))
			*res = _one_;
		else
			*res = _zero_;
		if (ztmp.len)
			zfree(ztmp);
		return;
	}
	if (zcmp(z1, rp->one) == 0) {
		if (sign)
			zsub(rp->mod, rp->one, res);
		else
			zcopy(rp->one, res);
		if (ztmp.len)
			zfree(ztmp);
		return;
	}

	/*
	 * See if the number being raised is the REDC format for -1.
	 * If so, then the answer is the REDC format for one or minus one.
	 * To do this check, calculate the REDC format for -1.
	 */
	if (((HALF)(z1.v[0] + rp->one.v[0])) == rp->mod.v[0]) {
		zsub(rp->mod, rp->one, &temp);
		if (zcmp(z1, temp) == 0) {
			if (zisodd(z2) ^ sign) {
				*res = temp;
				if (ztmp.len)
					zfree(ztmp);
				return;
			}
			zfree(temp);
			zcopy(rp->one, res);
			if (ztmp.len)
				zfree(ztmp);
			return;
		}
		zfree(temp);
	}

	for (pp = &lowpowers[2]; pp < &lowpowers[POWNUMS]; pp++)
		pp->len = 0;
	zcopy(rp->one, &lowpowers[0]);
	zcopy(z1, &lowpowers[1]);
	zcopy(rp->one, &ans);

	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
				zcopy(rp->one, &modpow);

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

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

	for (pp = lowpowers; pp < &lowpowers[POWNUMS]; pp++) {
		if (pp->len)
			freeh(pp->v);
	}
	if (sign && !ziszero(ans)) {
		zsub(rp->mod, ans, res);
		zfree(ans);
	} else {
		*res = ans;
	}
	if (ztmp.len)
		zfree(ztmp);
}


/*
 * zhnrmod - compute z mod h*2^n+r
 *
 * We compute v mod h*2^n+r, where h>0, n>0, abs(r) <= 1, as follows:
 *
 *	Let v = b*2^n + a, where 0 <= a < 2^n
 *
 *	Now v mod h*2^n+r == b*2^n + a mod h*2^n+r,
 *	and thus v mod h*2^n+r == b*2^n mod h*2^n+r + a mod h*2^n+r.
 *
 *	Because 0 <= a < 2^n < h*2^n+r, a mod h*2^n+r == a.
 *	Thus v mod h*2^n+r == b*2^n mod h*2^n+r + a.
 *
 *	It can be shown that b*2^n mod h*2^n == 2^n * (b mod h).
 *
 *	Thus for r == 0, v mod h*2^n+r == (2^n)*(b mod h) + a.
 *
 *	It can be shown that v mod 2^n-1 == a+b mod 2^n-1.
 *
 *	Thus for r == -1, v mod h*2^n+r == (2^n)*(b mod h) + a + int(b/h).
 *
 *	It can be shown that v mod 2^n+1 == a-b mod 2^n+1.
 *
 *	Thus for r == +1, v mod h*2^n+r == (2^n)*(b mod h) + a - int(b/h).
 *
 *	Therefore, v mod h*2^n+r == (2^n)*(b mod h) + a - r*int(b/h).
 *
 * The above proof leads to the following calc resource file which computes
 * the value z mod h*2^n+r:
 *
 *    define hnrmod(v,h,n,r)
 *    {
 *	local a,b,modulus,tquo,tmod,lbit,ret;
 *
 *	if (!isint(h) || h < 1) {
 *	    quit "h must be an integer be > 0";
 *	}
 *	if (!isint(n) || n < 1) {
 *	    quit "n must be an integer be > 0";
 *	}
 *	if (r != 1 && r != 0 && r != -1) {
 *	    quit "r must be -1, 0 or 1";
 *	}
 *
 *	lbit = lowbit(h);
 *	if (lbit > 0) {
 *	    n += lbit;
 *	    h >>= lbit;
 *	}
 *
 *	modulus = h<<n+r;
 *	if (modulus <= 2^31-1) {
 *	    return v % modulus;
 *	}
 *	ret = v;
 *
 *	do {
 *	    if (highbit(ret) < n) {
 *		break;
 *	    }
 *	    b = ret>>n;
 *	    a = ret - (b<<n);
 *
 *	    switch (r) {
 *	    case -1:
 *		if (h == 1) {
 *		    ret = a + b;
 *		} else {
 *		    quomod(b, h, tquo, tmod);
 *		    ret = tmod<<n + a + tquo;
 *		}
 *		break;
 *	    case 0:
 *		if (h == 1) {
 *		    ret = a;
 *		} else {
 *		    ret = (b%h)<<n + a;
 *		}
 *		break;
 *	    case 1:
 *		if (h == 1) {
 *		    ret = ((a > b) ? a-b : modulus+a-b);
 *		} else {
 *		    quomod(b, h, tquo, tmod);
 *		    tmod = tmod<<n + a;
 *		    ret = ((tmod >= tquo) ? tmod-tquo : modulus+tmod-tquo);
 *		}
 *		break;
 *	    }
 *	} while (ret > modulus);
 *	ret = ((ret < 0) ? ret+modlus : ((ret == modulus) ? 0 : ret));
 *
 *	return ret;
 *    }
 *
 * This function implements the above calc resource file.
 *
 * given:
 *	v		take mod of this value, v >= 0
 *	zh		h from modulus h*2^n+r, h > 0
 *	zn		n from modulus h*2^n+r, n > 0
 *	zr		r from modulus h*2^n+r, abs(r) <= 1
 *	res		v mod h*2^n+r
 */
void
zhnrmod(ZVALUE v, ZVALUE zh, ZVALUE zn, ZVALUE zr, ZVALUE *res)
{
	ZVALUE a;		/* lower n bits of v */
	ZVALUE b;		/* bits above the lower n bits of v */
	ZVALUE h;		/* working zh value */
	ZVALUE modulus;		/* h^2^n + r */
	ZVALUE tquo;		/* b // h */
	ZVALUE tmod;		/* b % h or (b%h)<<n + a */
	ZVALUE t;		/* temp ZVALUE */
	ZVALUE t2;		/* temp ZVALUE */
	ZVALUE ret;		/* return value, what *res is set to */
	long n;			/* integer value of zn */
	long r;			/* integer value of zr */
	long hbit;		/* highbit(res) */
	long lbit;		/* lowbit(h) */
	int zrelval;		/* return value of zrel() */
	int hisone;		/* 1 => h == 1, 0 => h != 1 */

	/*
	 * firewall
	 */
	if (zisneg(zh) || ziszero(zh)) {
		math_error("h must be > 0");
		/*NOTREACHED*/
	}
	if (zisneg(zn) || ziszero(zn)) {
		math_error("n must be > 0");
		/*NOTREACHED*/
	}
	if (zge31b(zn)) {
		math_error("n must be < 2^31");
		/*NOTREACHED*/
	}
	if (!zisabsleone(zr)) {
		math_error("r must be -1, 0 or 1");
		/*NOTREACHED*/
	}


	/*
	 * setup for loop
	 */
	n = ztolong(zn);
	r = ztolong(zr);
	if (zisneg(zr)) {
		r = -r;
	}
	/* lbit = lowbit(h); */
	lbit = zlowbit(zh);
	/* if (lbit > 0) { n += lbit; h >>= lbit; } */
	if (lbit > 0) {
		n += lbit;
		zshift(zh, -lbit, &h);
	} else {
		h = zh;
	}
	/* modulus = h<<n+r; */
	zshift(h, n, &t);
	switch (r) {
	case 1:
		zadd(t, _one_, &modulus);
		zfree(t);
		break;
	case 0:
		modulus = t;
		break;
	case -1:
		zsub(t, _one_, &modulus);
		zfree(t);
		break;
	}
	/* if (modulus <= MAXLONG) { return v % modulus; } */
	if (!zgtmaxlong(modulus)) {
		itoz(zmodi(v, ztolong(modulus)), res);
		zfree(modulus);
		if (lbit > 0) {
			zfree(h);
		}
		return;
	}
	/* ret = v; */
	zcopy(v, &ret);

	/*
	 * shift-add modulus loop
	 */
	hisone = zisone(h);
	do {

		/*
		 * split ret into to chunks, the lower n bits
		 * and everything above the lower n bits
		 */
		/* if (highbit(ret) < n) { break; } */
		hbit = (long)zhighbit(ret);
		if (hbit < n) {
			zrelval = (zcmp(ret, modulus) ? -1 : 0);
			break;
		}
		/* b = ret>>n; */
		zshift(ret, -n, &b);
		b.sign = ret.sign;
		/* a = ret - (b<<n); */
		a.sign = ret.sign;
		a.len = (n+BASEB-1)/BASEB;
		a.v = alloc(a.len);
		memcpy(a.v, ret.v, a.len*sizeof(HALF));
		if (n % BASEB) {
			a.v[a.len - 1] &= lowhalf[n % BASEB];
		}
		ztrim(&a);

		/*
		 * switch depending on r == -1, 0 or 1
		 */
		switch (r) {
		case -1:	/* v mod h*2^h-1 */
			/* if (h == 1) ... */
			if (hisone) {
				/* ret = a + b; */
				zfree(ret);
				zadd(a, b, &ret);

			/* ... else ... */
			} else {
				/* quomod(b, h, tquo, tmod); */
				(void) zdiv(b, h, &tquo, &tmod, 0);
				/* ret = tmod<<n + a + tquo; */
				zshift(tmod, n, &t);
				zfree(tmod);
				zadd(a, tquo, &t2);
				zfree(tquo);
				zfree(ret);
				zadd(t, t2, &ret);
				zfree(t);
				zfree(t2);
			}
			break;

		case 0:		/* v mod h*2^h-1 */
			/* if (h == 1) ... */
			if (hisone) {
				/* ret = a; */
				zfree(ret);
				zcopy(a, &ret);

			/* ... else ... */
			} else {
				/* ret = (b%h)<<n + a; */
				(void) zmod(b, h, &tmod, 0);
				zshift(tmod, n, &t);
				zfree(tmod);
				zfree(ret);
				zadd(t, a, &ret);
				zfree(t);
			}
			break;

		case 1:		/* v mod h*2^h-1 */
			/* if (h == 1) ... */
			if (hisone) {
				/* ret = a-b; */
				zfree(ret);
				zsub(a, b, &ret);

			/* ... else ... */
			} else {
				/* quomod(b, h, tquo, tmod); */
				(void) zdiv(b, h, &tquo, &tmod, 0);
				/* tmod = tmod<<n + a; */
				zshift(tmod, n, &t);
				zfree(tmod);
				zadd(t, a, &tmod);
				zfree(t);
				/* ret = tmod-tquo; */
				zfree(ret);
				zsub(tmod, tquo, &ret);
				zfree(tquo);
				zfree(tmod);
			}
			break;
		}
		zfree(a);
		zfree(b);

	/* ... while (abs(ret) > modulus); */
	} while ((zrelval = zabsrel(ret, modulus)) > 0);
	/* ret = ((ret < 0) ? ret+modlus : ((ret == modulus) ? 0 : ret)); */
	if (ret.sign) {
		zadd(ret, modulus, &t);
		zfree(ret);
		ret = t;
	} else if (zrelval == 0) {
		zfree(ret);
		ret = _zero_;
	}
	zfree(modulus);
	if (lbit > 0) {
		zfree(h);
	}

	/*
	 * return ret
	 */
	*res = ret;
	return;
}