zmath.c

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

	return ans;
}


/*
 * Compute the bitwise OR of two integers
 */
void
zor(ZVALUE z1, ZVALUE z2, ZVALUE *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 bitwise AND of two integers
 */
void
zand(ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
	HALF *h1, *h2, *hd;
	LEN 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 bitwise XOR of two integers.
 */
void
zxor(ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
	HALF *dp, *h1, *h2;
	LEN len, j, k;
	ZVALUE dest;

	h1 = z1.v;
	h2 = z2.v;
	len = z1.len;
	j = z2.len;
	if (z1.len < z2.len) {
		len = z2.len;
		j = z1.len;
		h1 = z2.v;
		h2 = z1.v;
	} else if (z1.len == z2.len) {
		while (len > 1 && z1.v[len-1] == z2.v[len-1])
			len--;
		j = len;
	}
	k = len - j;
	dest.len = len;
	dest.v = alloc(len);
	dest.sign = 0;
	dp = dest.v;
	while (j-- > 0)
		*dp++ = *h1++ ^ *h2++;
	while (k-- > 0)
		*dp++ = *h1++;
	*res = dest;
}


/*
 * Compute the bitwise ANDNOT of two integers.
 */
void
zandnot(ZVALUE z1, ZVALUE z2, ZVALUE *res)
{
	HALF *dp, *h1, *h2;
	LEN len, j, k;
	ZVALUE dest;

	len = z1.len;
	if (z2.len >= len) {
		while (len > 1 && (z1.v[len-1] & ~z2.v[len-1]) == 0)
			len--;
		j = len;
		k = 0;
	} else {
		j = z2.len;
		k = len - z2.len;
	}
	dest.len = len;
	dest.v = alloc(len);
	dest.sign = 0;
	dp = dest.v;
	h1 = z1.v;
	h2 = z2.v;
	while (j-- > 0)
		*dp++ = *h1++ & ~*h2++;
	while (k-- > 0)
		*dp++ = *h1++;
	*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(ZVALUE z, long n, ZVALUE *res)
{
	ZVALUE ans;
	LEN 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 = (LEN)(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 = (LEN)(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(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;
	/* ignore Saber-C warning #530 about empty while statement */
	/*	  ok to ignore in proc zlowbit */
	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).
 */
LEN
zhighbit(ZVALUE z)
{
	HALF dataval;
	HALF *bitval;

	dataval = z.v[z.len-1];
	if (dataval == 0) {
		return 0;
	}
	bitval = bitmask+BASEB;
	if (dataval) {
		/* ignore Saber-C warning #530 about empty while statement */
		/*	  ok to ignore in proc zhighbit */
		while ((*(--bitval) & dataval) == 0) {
		}
	}
	return (z.len*BASEB)+(LEN)(bitval-bitmask-BASEB);
}


/*
 * Return whether or not the specifed bit number is set in a number.
 * Rightmost bit of a number is bit 0.
 */
BOOL
zisset(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(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(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(long n, ZVALUE *res)
{
	ZVALUE z;

	if (n < 0) n = 0;
	z.sign = 0;
	z.len = (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(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;
}


/*
 * Return the sign of z1 - z2, i.e. 1 if the first integer is greater,
 * 0 if they are equal, -1 otherwise.
 */
FLAG
zrel(ZVALUE z1, ZVALUE z2)
{
	HALF *h1, *h2;
	LEN len;
	int sign;

	sign = 1;
	if (z1.sign < z2.sign)
		return 1;
	if (z2.sign < z1.sign)
		return -1;
	if (z2.sign)
		sign = -1;
	if (z1.len != z2.len)
		return (z1.len > z2.len) ? sign : -sign;
	len = z1.len;
	h1 = z1.v + len;
	h2 = z2.v + len;

	while (len > 0) {
		if (*--h1 != *--h2)
			break;
		len--;
	}
	if (len > 0)
		return (*h1 > *h2) ? sign : -sign;
	return 0;
}


/*
 * Return the sign of abs(z1) - abs(z2), i.e. 1 if the first integer
 * has greater absolute value, 0 is they have equal absolute value,
 * -1 otherwise.
 */
FLAG
zabsrel(ZVALUE z1, ZVALUE z2)
{
	HALF *h1, *h2;
	LEN len;

	if (z1.len != z2.len)
		return (z1.len > z2.len) ? 1 : -1;

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


/*
 * Compare two numbers to see if they are equal or not.
 * Returns TRUE if they differ.
 */
BOOL
zcmp(ZVALUE z1, ZVALUE 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;
}


/*
 * Utility to calculate the gcd of two FULL integers.
 */
FULL
uugcd(FULL f1, FULL f2)
{
	FULL temp;

	while (f1) {
		temp = f2 % f1;
		f2 = f1;
		f1 = temp;
	}
	return (FULL) f2;
}


/*
 * Utility to calculate the gcd of two small integers.
 */
long
iigcd(long i1, long 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(ZVALUE *z)
{
	HALF *h;
	LEN 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.
 *
 * NOTE: The ZVALUE length is not adjusted instead, the value is
 *	 zero padded from the left.  One may need to call ztrim()
 *	 or use zshift() instead.
 */
void
zshiftr(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;
		mask = 0;
		while (len--) {
			maskt = (((FULL) *--h) << (BASEB - n)) & BASE1;
			*h = ((*h >> n) | (HALF)mask);
			mask = maskt;
		}
	}
}


/*
 * Utility routine to shift left.
 *
 * NOTE: The ZVALUE length is not adjusted.  The bits in the upper
 *	 HALF are simply tossed.  You may want to use zshift() instead.
 */
void
zshiftl(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;
		}
	}
}


/*
 * popcnt - count the number of 0 or 1 bits in an integer
 *
 * We ignore all 0 bits above the highest bit.
 */
long
zpopcnt(ZVALUE z, int bitval)
{
	long cnt = 0;	/* number of times found */
	HALF h;		/* HALF to count */
	int i;

	/*
	 * count 1's
	 */
	if (bitval) {

		/*
		 * count each HALF
		 */
		for (i=0; i < z.len; ++i) {
			/* count each octet */
			for (h = z.v[i]; h; h >>= 8) {
				cnt += (long)popcnt[h & 0xff];
			}
		}

	/*
	 * count 0's
	 */
	} else {

		/*
		 * count each HALF up until the last
		 */
		for (i=0; i < z.len-1; ++i) {

			/* count each octet */
			cnt += BASEB;
			for (h = z.v[i]; h; h >>= 8) {
				cnt -= (long)popcnt[h & 0xff];
			}
		}

		/*
		 * count the last octet up until the highest 1 bit
		 */
		for (h = z.v[z.len-1]; h; h>>=1) {
			/* count each 0 bit */
			if ((h & 0x1) == 0) {
				++cnt;
			}
		}
	}

	/*
	 * return count
	 */
	return cnt;
}