zmath.c

早期freebsd实现

		zsub(z2, ztmp1, rem);
		zfree(ztmp1);
	} else
		*rem = ztmp1;
}


/*
 * Calculate the mod of an integer by a small number.
 * This is only defined for positive moduli.
 */
long
zmodi(z, n)
	ZVALUE z;
	long n;
{
	register HALF *h1;
	FULL val;
	HALF divval[2];
	ZVALUE div;
	ZVALUE temp;
	long len;

	if (n == 0)
		math_error("Division by zero");
	if (n < 0)
		math_error("Non-positive modulus");
	if (ziszero(z) || (n == 1))
		return 0;
	if (zisone(z))
		return 1;
	/*
	 * If the modulus is by a large number, then call the normal
	 * modulo routine.
	 */
	if (n & ~BASE1) {
		div.sign = 0;
		div.len = 2;
		div.v = divval;
		divval[0] = (HALF) n;
		divval[1] = ((FULL) n) >> BASEB;
		zmod(z, div, &temp);
		n = (zistiny(temp) ? z1tol(temp) : z2tol(temp));
		zfree(temp);
		return n;
	}
	/*
	 * The modulus is by a small number, so we can do this quickly.
	 */
	len = z.len;
	h1 = z.v + len - 1;
	val = 0;
	while (len--)
		val = ((val << BASEB) + ((FULL) *h1--)) % n;
	if (z.sign)
		val = n - val;
	return val;
}


/*
 * Return whether or not one number exactly divides another one.
 * Returns TRUE if division occurs with no remainder.
 * z1 is the number to be divided by z2.
 */
BOOL
zdivides(z1, z2)
	ZVALUE z1, z2;		/* numbers to test division into and by */
{
	ZVALUE temp;
	long cv;

	z1.sign = 0;
	z2.sign = 0;
	/*
	 * Take care of obvious cases first
	 */
	if (zisleone(z2)) {	/* division by zero or one */
		if (*z2.v == 0)
			math_error("Division by zero");
		return TRUE;
	}
	if (ziszero(z1))	/* everything divides zero */
		return TRUE;
	if (z1.len < z2.len)	/* quick size comparison */
		return FALSE;
	if ((z1.len == z2.len) && (z1.v[z1.len-1] < z2.v[z2.len-1]))	/* more */
		return FALSE;
	if (zisodd(z1) && ziseven(z2))	/* can't divide odd by even */
		return FALSE;
	if (zlowbit(z1) < zlowbit(z2))	/* can't have smaller power of two */
		return FALSE;
	cv = zrel(z1, z2);	/* can't divide smaller number */
	if (cv <= 0)
		return (cv == 0);
	/*
	 * Now do the real work.  Divisor divides dividend if the gcd of the
	 * two numbers equals the divisor.
	 */
	zgcd(z1, z2, &temp);
	cv = zcmp(z2, temp);
	zfree(temp);
	return (cv == 0);
}


/*
 * Compute the logical OR of two numbers
 */
void
zor(z1, z2, res)
	ZVALUE z1, z2, *res;
{
	register HALF *sp, *dp;
	long len;
	ZVALUE bz, lz, dest;

	if (z1.len >= z2.len) {
		bz = z1;
		lz = z2;
	} else {
		bz = z2;
		lz = z1;
	}
	dest.len = bz.len;
	dest.v = alloc(dest.len);
	dest.sign = 0;
	zcopyval(bz, dest);
	len = lz.len;
	sp = lz.v;
	dp = dest.v;
	while (len--)
		*dp++ |= *sp++;
	*res = dest;
}


/*
 * Compute the logical AND of two numbers.
 */
void
zand(z1, z2, res)
	ZVALUE z1, z2, *res;
{
	register HALF *h1, *h2, *hd;
	register long len;
	ZVALUE dest;

	len = ((z1.len <= z2.len) ? z1.len : z2.len);
	h1 = &z1.v[len-1];
	h2 = &z2.v[len-1];
	while ((len > 1) && ((*h1 & *h2) == 0)) {
		h1--;
		h2--;
		len--;
	}
	dest.len = len;
	dest.v = alloc(len);
	dest.sign = 0;
	h1 = z1.v;
	h2 = z2.v;
	hd = dest.v;
	while (len--)
		*hd++ = (*h1++ & *h2++);
	*res = dest;
}


/*
 * Compute the logical XOR of two numbers.
 */
void
zxor(z1, z2, res)
	ZVALUE z1, z2, *res;
{
	register HALF *sp, *dp;
	long len;
	ZVALUE bz, lz, dest;

	if (z1.len == z2.len) {
		for (len = z1.len; ((len > 1) && (z1.v[len-1] == z2.v[len-1])); len--) ;
		z1.len = len;
		z2.len = len;
	}
	if (z1.len >= z2.len) {
		bz = z1;
		lz = z2;
	} else {
		bz = z2;
		lz = z1;
	}
	dest.len = bz.len;
	dest.v = alloc(dest.len);
	dest.sign = 0;
	zcopyval(bz, dest);
	len = lz.len;
	sp = lz.v;
	dp = dest.v;
	while (len--)
		*dp++ ^= *sp++;
	*res = dest;
}


/*
 * Shift a number left (or right) by the specified number of bits.
 * Positive shift means to the left.  When shifting right, rightmost
 * bits are lost.  The sign of the number is preserved.
 */
void
zshift(z, n, res)
	ZVALUE z, *res;
	long n;
{
	ZVALUE ans;
	long hc;		/* number of halfwords shift is by */

	if (ziszero(z)) {
		*res = _zero_;
		return;
	}
	if (n == 0) {
		zcopy(z, res);
		return;
	}
	/*
	 * If shift value is negative, then shift right.
	 * Check for large shifts, and handle word-sized shifts quickly.
	 */
	if (n < 0) {
		n = -n;
		if ((n < 0) || (n >= (z.len * BASEB))) {
			*res = _zero_;
			return;
		}
		hc = n / BASEB;
		n %= BASEB;
		z.v += hc;
		z.len -= hc;
		ans.len = z.len;
		ans.v = alloc(ans.len);
		ans.sign = z.sign;
		zcopyval(z, ans);
		if (n > 0) {
			zshiftr(ans, n);
			ztrim(&ans);
		}
		if (ziszero(ans)) {
			zfree(ans);
			ans = _zero_;
		}
		*res = ans;
		return;
	}
	/*
	 * Shift value is positive, so shift leftwards.
	 * Check specially for a shift of the value 1, since this is common.
	 * Also handle word-sized shifts quickly.
	 */
	if (zisunit(z)) {
		zbitvalue(n, res);
		res->sign = z.sign;
		return;
	}
	hc = n / BASEB;
	n %= BASEB;
	ans.len = z.len + hc + 1;
	ans.v = alloc(ans.len);
	ans.sign = z.sign;
	if (hc > 0)
		memset((char *) ans.v, 0, hc * sizeof(HALF));
	memcpy((char *) (ans.v + hc), 
	    (char *) z.v, z.len * sizeof(HALF));
	ans.v[ans.len - 1] = 0;
	if (n > 0) {
		ans.v += hc;
		ans.len -= hc;
		zshiftl(ans, n);
		ans.v -= hc;
		ans.len += hc;
	}
	ztrim(&ans);
	*res = ans;
}


/*
 * Return the position of the lowest bit which is set in the binary
 * representation of a number (counting from zero).  This is the highest
 * power of two which evenly divides the number.
 */
long
zlowbit(z)
	ZVALUE z;
{
	register HALF *zp;
	long n;
	HALF dataval;
	HALF *bitval;

	n = 0;
	for (zp = z.v; *zp == 0; zp++)
		if (++n >= z.len)
			return 0;
	dataval = *zp;
	bitval = bitmask;
	while ((*(bitval++) & dataval) == 0) {
	}
	return (n*BASEB)+(bitval-bitmask-1);
}


/*
 * Return the position of the highest bit which is set in the binary
 * representation of a number (counting from zero).  This is the highest power
 * of two which is less than or equal to the number (which is assumed nonzero).
 */
long
zhighbit(z)
	ZVALUE z;
{
	HALF dataval;
	HALF *bitval;

	dataval = z.v[z.len-1];
	bitval = bitmask+BASEB;
	if (dataval) {
		while ((*(--bitval) & dataval) == 0) {
		}
	}
	return (z.len*BASEB)+(bitval-bitmask-BASEB);
}


#if 0
/*
 * Reverse the bits of a particular range of bits of a number.
 *
 * This function returns an integer with bits a thru b swapped.
 * That is, bit a is swapped with bit b, bit a+1 is swapped with b-1,
 * and so on.
 *
 * As a special case, if the ending bit position is < 0, is it taken to 
 * mean the highest bit set.  Thus zbitrev(0, -1, z, &res) will 
 * perform a complete bit reverse of the number 'z'.
 *
 * As a special case, if the starting bit position is < 0, is it taken to 
 * mean the lowest bit set.  Thus zbitrev(-1, -1, z, &res) is the
 * same as zbitrev(lowbit(z), highbit(z), z, &res).
 *
 * Note that the low order bit number is taken to be 0.  Also, bitrev
 * ignores the sign of the number.
 *
 * Bits beyond the highest bit are taken to be zero.  Thus the calling
 * bitrev(0, 100, _one_, &res) will result in a value of 2^100.
 */
void
zbitrev(low, high, z, res)
	long low;	/* lowest bit to reverse, <0 => lowbit(z) */
	long high;	/* highest bit to reverse, <0 => highbit(z) */
	ZVALUE z;	/* value to bit reverse */
	ZVALUE *res;	/* resulting bit reverse number */
{
}
#endif


/*
 * Return whether or not the specifed bit number is set in a number.
 * Rightmost bit of a number is bit 0.
 */
BOOL
zisset(z, n)
	ZVALUE z;
	long n;
{
	if ((n < 0) || ((n / BASEB) >= z.len))
		return FALSE;
	return ((z.v[n / BASEB] & (((HALF) 1) << (n % BASEB))) != 0);
}


/*
 * Check whether or not a number has exactly one bit set, and
 * thus is an exact power of two.  Returns TRUE if so.
 */
BOOL
zisonebit(z)
	ZVALUE z;
{
	register HALF *hp;
	register LEN len;

	if (ziszero(z) || zisneg(z))
		return FALSE;
	hp = z.v;
	len = z.len;
	while (len > 4) {
		len -= 4;
		if (*hp++ || *hp++ || *hp++ || *hp++)
			return FALSE;
	}
	while (--len > 0) {
		if (*hp++)
			return FALSE;
	}
	return ((*hp & -*hp) == *hp);		/* NEEDS 2'S COMPLEMENT */
}


/*
 * Check whether or not a number has all of its bits set below some
 * bit position, and thus is one less than an exact power of two.
 * Returns TRUE if so.
 */
BOOL
zisallbits(z)
	ZVALUE z;
{
	register HALF *hp;
	register LEN len;
	HALF digit;

	if (ziszero(z) || zisneg(z))
		return FALSE;
	hp = z.v;
	len = z.len;
	while (len > 4) {
		len -= 4;
		if ((*hp++ != BASE1) || (*hp++ != BASE1) ||
			(*hp++ != BASE1) || (*hp++ != BASE1))
				return FALSE;
	}
	while (--len > 0) {
		if (*hp++ != BASE1)
			return FALSE;
	}
	digit = (HALF)(*hp + 1);
	return ((digit & -digit) == digit);	/* NEEDS 2'S COMPLEMENT */
}


/*
 * Return the number whose binary representation contains only one bit which
 * is in the specified position (counting from zero).  This is equivilant
 * to raising two to the given power.
 */
void
zbitvalue(n, res)
	long n;
	ZVALUE *res;
{
	ZVALUE z;

	if (n < 0) n = 0;
	z.sign = 0;
	z.len = (n / BASEB) + 1;
	z.v = alloc(z.len);
	zclearval(z);
	z.v[z.len-1] = (((HALF) 1) << (n % BASEB));
	*res = z;
}


/*
 * Compare a number against zero.
 * Returns the sgn function of the number (-1, 0, or 1).
 */
FLAG
ztest(z)
	ZVALUE z;
{
	register int sign;
	register HALF *h;
	register long len;

	sign = 1;
	if (z.sign)
		sign = -sign;
	h = z.v;
	len = z.len;
	while (len--)
		if (*h++)
			return sign;
	return 0;
}


/*
 * Compare two numbers to see which is larger.
 * Returns -1 if first number is smaller, 0 if they are equal, and 1 if
 * first number is larger.  This is the same result as ztest(z2-z1).
 */
FLAG
zrel(z1, z2)
	ZVALUE z1, z2;
{
	register HALF *h1, *h2;
	register long len1, len2;
	int sign;

	sign = 1;
	if (z1.sign < z2.sign)
		return 1;
	if (z2.sign < z1.sign)
		return -1;
	if (z2.sign)
		sign = -1;
	len1 = z1.len;
	len2 = z2.len;
	h1 = z1.v + z1.len - 1;
	h2 = z2.v + z2.len - 1;
	while (len1 > len2) {
		if (*h1--)
			return sign;
		len1--;
	}
	while (len2 > len1) {
		if (*h2--)
			return -sign;
		len2--;
	}
	while (len1--) {
		if (*h1-- != *h2--)
			break;
	}
	if ((len1 = *++h1) > (len2 = *++h2))
		return sign;
	if (len1 < len2)
		return -sign;
	return 0;
}


/*
 * Compare two numbers to see if they are equal or not.
 * Returns TRUE if they differ.
 */
BOOL
zcmp(z1, z2)
	ZVALUE z1, z2;
{
	register HALF *h1, *h2;
	register long len;

	if ((z1.sign != z2.sign) || (z1.len != z2.len) || (*z1.v != *z2.v))
		return TRUE;
	len = z1.len;
	h1 = z1.v;
	h2 = z2.v;
	while (len-- > 0) {
		if (*h1++ != *h2++)
			return TRUE;
	}
	return FALSE;
}


/*
 * Internal utility subroutines
 */
static void
dadd(z1, z2, y, n)
	ZVALUE z1, z2;
	long y, n;
{
	HALF *s1, *s2;
	short carry;
	long sum;

	s1 = z1.v + y - n;
	s2 = z2.v;
	carry = 0;
	while (n--) {
		sum = (long)*s1 + (long)*s2 + carry;
		carry = 0;
		if (sum >= BASE) {
			sum -= BASE;
			carry = 1;
		}
		*s1 = (HALF)sum;
		++s1;
		++s2;
	}
	sum = (long)*s1 + carry;
	*s1 = (HALF)sum;
}


/*
 * Do subtract for divide, returning TRUE if subtraction went negative.
 */
static BOOL
dsub(z1, z2, y, n)
	ZVALUE z1, z2;
	long y, n;
{
	HALF *s1, *s2, *s3;
	FULL i1;
	BOOL neg;

	neg = FALSE;
	s1 = z1.v + y - n;
	s2 = z2.v;
	if (++n > z2.len)
		n = z2.len;
	while (n--) {
		i1 = (FULL) *s1;
		if (i1 < (FULL) *s2) {
			s3 = s1 + 1;
			while (s3 < z1.v + z1.len && !(*s3)) {
				*s3 = BASE1;
				++s3;
			}
			if (s3 >= z1.v + z1.len)
				neg = TRUE;
			else
				--(*s3);
			i1 += BASE;
		}
		*s1 = i1 - (FULL) *s2;
		++s1;
		++s2;
	}
	return neg;
}


/*
 * Multiply a number by a single 'digit'.
 * This is meant to be used only by the divide routine, and so the
 * destination area must already be allocated and be large enough.
 */
static void
dmul(z, mul, dest)
	ZVALUE z;
	FULL mul;
	ZVALUE *dest;
{
	register HALF *zp, *dp;
	SIUNION sival;
	FULL carry;
	long len;

	dp = dest->v;
	dest->sign = 0;
	if (mul == 0) {
		dest->len = 1;
		*dp = 0;
		return;
	}
	len = z.len;
	zp = z.v + len - 1;
	while ((*zp == 0) && (len > 1)) {
		len--;
		zp--;
	}
	dest->len = len;
	zp = z.v;
	carry = 0;
	while (len >= 4) {
		len -= 4;
		sival.ivalue = (mul * ((FULL) *zp++)) + carry;
		*dp++ = sival.silow;
		sival.ivalue = (mul * ((FULL) *zp++)) + ((FULL) sival.sihigh);
		*dp++ = sival.silow;
		sival.ivalue = (mul * ((FULL) *zp++)) + ((FULL) sival.sihigh);
		*dp++ = sival.silow;
		sival.ivalue = (mul * ((FULL) *zp++)) + ((FULL) sival.sihigh);
		*dp++ = sival.silow;
		carry = sival.sihigh;
	}
	while (--len >= 0) {
		sival.ivalue = (mul * ((FULL) *zp++)) + carry;
		*dp++ = sival.silow;
		carry = sival.sihigh;
	}
	if (carry) {
		*dp = (HALF)carry;
		dest->len++;
	}
}


/*
 * Utility to calculate the gcd of two small integers.
 */
long
iigcd(i1, i2)
	long i1, i2;
{
	FULL f1, f2, temp;

	f1 = (FULL) ((i1 >= 0) ? i1 : -i1);
	f2 = (FULL) ((i2 >= 0) ? i2 : -i2);
	while (f1) {
		temp = f2 % f1;
		f2 = f1;
		f1 = temp;
	}
	return (long) f2;
}


void
ztrim(z)
	ZVALUE *z;
{
	register HALF *h;
	register long len;

	h = z->v + z->len - 1;
	len = z->len;
	while (*h == 0 && len > 1) {
		--h;
		--len;
	}
	z->len = len;
}


/*
 * Utility routine to shift right.
 */
void
zshiftr(z, n)
	ZVALUE z;
	long n;
{
	register HALF *h, *lim;
	FULL mask, maskt;
	long len;

	if (n >= BASEB) {
		len = n / BASEB;
		h = z.v;
		lim = z.v + z.len - len;
		while (h < lim) {
			h[0] = h[len];
			++h;
		}
		n -= BASEB * len;
		lim = z.v + z.len;
		while (h < lim)
			*h++ = 0;
	}
	if (n) {
		len = z.len;
		h = z.v + len - 1;
		mask = 0;
		while (len--) {
			maskt = (((FULL) *h) << (BASEB - n)) & BASE1;
			*h = (*h >> n) | mask;
			mask = maskt;
			--h;
		}
	}
}


/*
 * Utility routine to shift left.
 */
void
zshiftl(z, n)
	ZVALUE z;
	long n;
{
	register HALF *h;
	FULL mask, i;
	long len;

	if (n >= BASEB) {
		len = n / BASEB;
		h = z.v + z.len - 1;
		while (!*h)
			--h;
		while (h >= z.v) {
			h[len] = h[0];
			--h;
		}
		n -= BASEB * len;
		while (len)
			h[len--] = 0;
	}
	if (n > 0) {
		len = z.len;
		h = z.v;
		mask = 0;
		while (len--) {
			i = (((FULL) *h) << n) | mask;
			if (i > BASE1) {
				mask = i >> BASEB;
				i &= BASE1;
			} else
				mask = 0;
			*h = (HALF) i;
			++h;
		}
	}
}

/*
 * initmasks - init the bitmask rotation arrays
 *
 * bitmask[i] 	 (1 << (i-1)),  for  -BASEB*4<=i<=BASEB*4
 *
 * The bmask array contains 8 cycles of rotations of a bit mask.
 */
void
initmasks()
{
	int i;

	/*
	 * setup the bmask array
	 */
	bmask = alloc((long)((8*BASEB)+1));
	for (i=0; i < (8*BASEB)+1; ++i) {
		bmask[i] = 1 << (i%BASEB);
	}

	/*
	 * setup the rmask pointers
	 */
	rmask = (HALF **)malloc(sizeof(HALF *)*((BASEB*4)+2));
	for (i = 0; i <= (4*BASEB)+1; ++i) {
		rmask[i] = &bmask[(2*BASEB)+i];
	}

	/*
	 * setup the bitmask array to allow -4*BASEB thru 4*BASEB indexing
	 */
	bitmask = &bmask[4*BASEB];
	return;
}

/* END CODE */