zmath.c

早期freebsd实现

/*
 * Copyright (c) 1994 David I. Bell
 * Permission is granted to use, distribute, or modify this source,
 * provided that this copyright notice remains intact.
 *
 * Extended precision integral arithmetic primitives
 */

#include "zmath.h"


HALF _twoval_[] = { 2 };
HALF _zeroval_[] = { 0 };
HALF _oneval_[] = { 1 };
HALF _tenval_[] = { 10 };

ZVALUE _zero_ = { _zeroval_, 1, 0};
ZVALUE _one_ = { _oneval_, 1, 0 };
ZVALUE _ten_ = { _tenval_, 1, 0 };


/*
 * mask of given bits, rotated thru all bit positions twice
 *
 * bitmask[i] 	 (1 << (i-1)),  for  -BASEB*4<=i<=BASEB*4
 */
static HALF *bmask;		/* actual rotation thru 8 cycles */
static HALF **rmask;		/* actual rotation pointers thru 2 cycles */
HALF *bitmask;			/* bit rotation, norm 0 */

BOOL _math_abort_;		/* nonzero to abort calculations */


static void dadd MATH_PROTO((ZVALUE z1, ZVALUE z2, long y, long n));
static BOOL dsub MATH_PROTO((ZVALUE z1, ZVALUE z2, long y, long n));
static void dmul MATH_PROTO((ZVALUE z, FULL x, ZVALUE *dest));


#ifdef ALLOCTEST
static long nalloc = 0;
static long nfree = 0;
#endif


HALF *
alloc(len)
	long len;
{
	HALF *hp;

	if (_math_abort_)
		math_error("Calculation aborted");
	hp = (HALF *) malloc((len+1) * sizeof(HALF));
	if (hp == 0)
		math_error("Not enough memory");
#ifdef ALLOCTEST
	++nalloc;
#endif
	return hp;
}


#ifdef ALLOCTEST
void
freeh(h)
	HALF *h;
{
	if ((h != _zeroval_) && (h != _oneval_)) {
		free(h);
		++nfree;
	}
}


void
allocStat()
{
	fprintf(stderr, "nalloc: %ld nfree: %ld kept: %ld\n",
		nalloc, nfree, nalloc - nfree);
}
#endif


/*
 * Convert a normal integer to a number.
 */
void
itoz(i, res)
	long i;
	ZVALUE *res;
{
	long diddle, len;

	res->len = 1;
	res->sign = 0;
	diddle = 0;
	if (i == 0) {
		res->v = _zeroval_;
		return;
	}
	if (i < 0) {
		res->sign = 1;
		i = -i;
		if (i < 0) {	/* fix most negative number */
			diddle = 1;
			i--;
		}
	}
	if (i == 1) {
		res->v = _oneval_;
		return;
	}
	len = 1 + (((FULL) i) >= BASE);
	res->len = len;
	res->v = alloc(len);
	res->v[0] = (HALF) (i + diddle);
	if (len == 2)
		res->v[1] = (HALF) (i / BASE);
}


/*
 * Convert a number to a normal integer, as far as possible.
 * If the number is out of range, the largest number is returned.
 */
long
ztoi(z)
	ZVALUE z;
{
	long i;

	if (zisbig(z)) {
		i = MAXFULL;
		return (z.sign ? -i : i);
	}
	i = (zistiny(z) ? z1tol(z) : z2tol(z));
	return (z.sign ? -i : i);
}


/*
 * Make a copy of an integer value
 */
void
zcopy(z, res)
	ZVALUE z, *res;
{
	res->sign = z.sign;
	res->len = z.len;
	if (zisleone(z)) {	/* zero or plus or minus one are easy */
		res->v = (z.v[0] ? _oneval_ : _zeroval_);
		return;
	}
	res->v = alloc(z.len);
	zcopyval(z, *res);
}


/*
 * Add together two integers.
 */
void
zadd(z1, z2, res)
	ZVALUE z1, z2, *res;
{
	ZVALUE dest;
	HALF *p1, *p2, *pd;
	long len;
	FULL carry;
	SIUNION sival;

	if (z1.sign && !z2.sign) {
		z1.sign = 0;
		zsub(z2, z1, res);
		return;
	}
	if (z2.sign && !z1.sign) {
		z2.sign = 0;
		zsub(z1, z2, res);
		return;
	}
	if (z2.len > z1.len) {
		pd = z1.v; z1.v = z2.v; z2.v = pd;
		len = z1.len; z1.len = z2.len; z2.len = len;
	}
	dest.len = z1.len + 1;
	dest.v = alloc(dest.len);
	dest.sign = z1.sign;
	carry = 0;
	pd = dest.v;
	p1 = z1.v;
	p2 = z2.v;
	len = z2.len;
	while (len--) {
		sival.ivalue = ((FULL) *p1++) + ((FULL) *p2++) + carry;
		*pd++ = sival.silow;
		carry = sival.sihigh;
	}
	len = z1.len - z2.len;
	while (len--) {
		sival.ivalue = ((FULL) *p1++) + carry;
		*pd++ = sival.silow;
		carry = sival.sihigh;
	}
	*pd = (HALF)carry;
	zquicktrim(dest);
	*res = dest;
}


/*
 * Subtract two integers.
 */
void
zsub(z1, z2, res)
	ZVALUE z1, z2, *res;
{
	register HALF *h1, *h2, *hd;
	long len1, len2;
	FULL carry;
	SIUNION sival;
	ZVALUE dest;

	if (z1.sign != z2.sign) {
		z2.sign = z1.sign;
		zadd(z1, z2, res);
		return;
	}
	len1 = z1.len;
	len2 = z2.len;
	if (len1 == len2) {
		h1 = z1.v + len1 - 1;
		h2 = z2.v + len2 - 1;
		while ((len1 > 0) && ((FULL)*h1 == (FULL)*h2)) {
			len1--;
			h1--;
			h2--;
		}
		if (len1 == 0) {
			*res = _zero_;
			return;
		}
		len2 = len1;
		carry = ((FULL)*h1 < (FULL)*h2);
	} else {
		carry = (len1 < len2);
	}
	dest.sign = z1.sign;
	h1 = z1.v;
	h2 = z2.v;
	if (carry) {
		carry = len1;
		len1 = len2;
		len2 = carry;
		h1 = z2.v;
		h2 = z1.v;
		dest.sign = !dest.sign;
	}
	hd = alloc(len1);
	dest.v = hd;
	dest.len = len1;
	len1 -= len2;
	carry = 0;
	while (--len2 >= 0) {
		sival.ivalue = (BASE1 - ((FULL) *h1++)) + *h2++ + carry;
		*hd++ = BASE1 - sival.silow;
		carry = sival.sihigh;
	}
	while (--len1 >= 0) {
		sival.ivalue = (BASE1 - ((FULL) *h1++)) + carry;
		*hd++ = BASE1 - sival.silow;
		carry = sival.sihigh;
	}
	if (hd[-1] == 0)
		ztrim(&dest);
	*res = dest;
}


/*
 * Multiply an integer by a small number.
 */
void
zmuli(z, n, res)
	ZVALUE z;
	long n;
	ZVALUE *res;
{
	register HALF *h1, *sd;
	FULL low;
	FULL high;
	FULL carry;
	long len;
	SIUNION sival;
	ZVALUE dest;

	if ((n == 0) || ziszero(z)) {
		*res = _zero_;
		return;
	}
	if (n < 0) {
		n = -n;
		z.sign = !z.sign;
	}
	if (n == 1) {
		zcopy(z, res);
		return;
	}
	low = ((FULL) n) & BASE1;
	high = ((FULL) n) >> BASEB;
	dest.len = z.len + 2;
	dest.v = alloc(dest.len);
	dest.sign = z.sign;
	/*
	 * Multiply by the low digit.
	 */
	h1 = z.v;
	sd = dest.v;
	len = z.len;
	carry = 0;
	while (len--) {
		sival.ivalue = ((FULL) *h1++) * low + carry;
		*sd++ = sival.silow;
		carry = sival.sihigh;
	}
	*sd = (HALF)carry;
	/*
	 * If there was only one digit, then we are all done except
	 * for trimming the number if there was no last carry.
	 */
	if (high == 0) {
		dest.len--;
		if (carry == 0)
			dest.len--;
		*res = dest;
		return;
	}
	/*
	 * Need to multiply by the high digit and add it into the
	 * previous value.  Clear the final word of rubbish first.
	 */
	*(++sd) = 0;
	h1 = z.v;
	sd = dest.v + 1;
	len = z.len;
	carry = 0;
	while (len--) {
		sival.ivalue = ((FULL) *h1++) * high + ((FULL) *sd) + carry;
		*sd++ = sival.silow;
		carry = sival.sihigh;
	}
	*sd = (HALF)carry;
	zquicktrim(dest);
	*res = dest;
}


/*
 * Divide two numbers by their greatest common divisor.
 * This is useful for reducing the numerator and denominator of
 * a fraction to its lowest terms.
 */
void
zreduce(z1, z2, z1res, z2res)
	ZVALUE z1, z2, *z1res, *z2res;
{
	ZVALUE tmp;

	if (zisleone(z1) || zisleone(z2))
		tmp = _one_;
	else
		zgcd(z1, z2, &tmp);
	if (zisunit(tmp)) {
		zcopy(z1, z1res);
		zcopy(z2, z2res);
	} else {
		zquo(z1, tmp, z1res);
		zquo(z2, tmp, z2res);
	}
	zfree(tmp);
}


/*
 * Divide two numbers to obtain a quotient and remainder.
 * This algorithm is taken from
 * Knuth, The Art of Computer Programming, vol 2: Seminumerical Algorithms.
 * Slight modifications were made to speed this mess up.
 */
void
zdiv(z1, z2, res, rem)
	ZVALUE z1, z2, *res, *rem;
{
	long i, j, k;
	register HALF *q, *pp;
	SIUNION pair;		/* pair of halfword values */
	HALF h2, v2;
	long y;
	FULL x;
	ZVALUE ztmp1, ztmp2, ztmp3, quo;

	if (ziszero(z2))
		math_error("Division by zero");
	if (ziszero(z1)) {
		*res = _zero_;
		*rem = _zero_;
		return;
	}
	if (zisone(z2)) {
		zcopy(z1, res);
		*rem = _zero_;
		return;
	}
	i = BASE / 2;
	j = 0;
	k = (long) z2.v[z2.len - 1];
	while (! (k & i)) {
		j ++;
		i >>= 1;
	}
	ztmp1.v = alloc(z1.len + 1);
	ztmp1.len = z1.len + 1;
	zcopyval(z1, ztmp1);
	ztmp1.v[z1.len] = 0;
	ztmp1.sign = 0;
	ztmp2.v = alloc(z2.len);
	ztmp2.len = z2.len;
	ztmp2.sign = 0;
	zcopyval(z2, ztmp2);
	if (zrel(ztmp1, ztmp2) < 0) {
		rem->v = ztmp1.v;
		rem->sign = z1.sign;
		rem->len = z1.len;
		zfree(ztmp2);
		*res = _zero_;
		return;
	}
	quo.len = z1.len - z2.len + 1;
	quo.v = alloc(quo.len);
	quo.sign = z1.sign != z2.sign;
	zclearval(quo);

	ztmp3.v = zalloctemp(z2.len + 1);

	/*
	 * Normalize z1 and z2
	 */
	zshiftl(ztmp1, j);
	zshiftl(ztmp2, j);

	k = ztmp1.len - ztmp2.len;
	q = quo.v + quo.len;
	y = ztmp1.len - 1;
	h2 = ztmp2.v [ztmp2.len - 1];
	v2 = 0;
	if (ztmp2.len >= 2)
		v2 = ztmp2.v [ztmp2.len - 2];
	for (;k--; --y) {
		pp = ztmp1.v + y - 1;
		pair.silow = pp[0];
		pair.sihigh = pp[1];
		if (ztmp1.v[y] == h2)
			x = BASE1;
		else
			x = pair.ivalue / h2;
		if (x) {
			while (pair.ivalue - x * h2 < BASE && y > 1 &&
				v2 * x > (pair.ivalue - x * h2) * BASE + ztmp1.v [y-2]) {
					--x;
			}
			dmul(ztmp2, x, &ztmp3);
#ifdef divblab
			printf(" x: %ld\n", x);
			printf("ztmp1: ");
			printz(ztmp1);
			printf("ztmp2: ");
			printz(ztmp2);
			printf("ztmp3: ");
			printz(ztmp3);
#endif
			if (dsub(ztmp1, ztmp3, y, ztmp2.len)) {
				--x;
				/*
				printf("adding back\n");
				*/
				dadd(ztmp1, ztmp2, y, ztmp2.len);
			}
		}
		ztrim(&ztmp1);
		*--q = (HALF)x;
	}
	zshiftr(ztmp1, j);
	*rem = ztmp1;
	ztrim(rem);
	zfree(ztmp2);
	ztrim(&quo);
	*res = quo;
}


/*
 * Return the quotient and remainder of an integer divided by a small
 * number.  A nonzero remainder is only meaningful when both numbers
 * are positive.
 */
long
zdivi(z, n, res)
	ZVALUE z, *res;
	long n;
{
	register HALF *h1, *sd;
	FULL val;
	HALF divval[2];
	ZVALUE div;
	ZVALUE dest;
	long len;

	if (n == 0)
		math_error("Division by zero");
	if (ziszero(z)) {
		*res = _zero_;
		return 0;
	}
	if (n < 0) {
		n = -n;
		z.sign = !z.sign;
	}
	if (n == 1) {
		zcopy(z, res);
		return 0;
	}
	/*
	 * If the division is by a large number, then call the normal
	 * divide routine.
	 */
	if (n & ~BASE1) {
		div.sign = 0;
		div.len = 2;
		div.v = divval;
		divval[0] = (HALF) n;
		divval[1] = ((FULL) n) >> BASEB;
		zdiv(z, div, res, &dest);
		n = (zistiny(dest) ? z1tol(dest) : z2tol(dest));
		zfree(dest);
		return n;
	}
	/*
	 * Division is by a small number, so we can be quick about it.
	 */
	len = z.len;
	dest.sign = z.sign;
	dest.len = len;
	dest.v = alloc(len);
	h1 = z.v + len - 1;
	sd = dest.v + len - 1;
	val = 0;
	while (len--) {
		val = ((val << BASEB) + ((FULL) *h1--));
		*sd-- = val / n;
		val %= n;
	}
	zquicktrim(dest);
	*res = dest;
	return val;
}


/*
 * Return the quotient of two numbers.
 * This works the same as zdiv, except that the remainer is not returned.
 */
void
zquo(z1, z2, res)
	ZVALUE z1, z2, *res;
{
	long i, j, k;
	register HALF *q, *pp;
	SIUNION pair;			/* pair of halfword values */
	HALF h2, v2;
	long y;
	FULL x;
	ZVALUE ztmp1, ztmp2, ztmp3, quo;

	if (ziszero(z2))
		math_error("Division by zero");
	if (ziszero(z1)) {
		*res = _zero_;
		return;
	}
	if (zisone(z2)) {
		zcopy(z1, res);
		return;
	}
	i = BASE / 2;
	j = 0;
	k = (long) z2.v[z2.len - 1];
	while (! (k & i)) {
		j ++;
		i >>= 1;
	}
	ztmp1.v = alloc(z1.len + 1);
	ztmp1.len = z1.len + 1;
	zcopyval(z1, ztmp1);
	ztmp1.v[z1.len] = 0;
	ztmp1.sign = 0;
	ztmp2.v = alloc(z2.len);
	ztmp2.len = z2.len;
	ztmp2.sign = 0;
	zcopyval(z2, ztmp2);
	if (zrel(ztmp1, ztmp2) < 0) {
		zfree(ztmp1);
		zfree(ztmp2);
		*res = _zero_;
		return;
	}
	quo.len = z1.len - z2.len + 1;
	quo.v = alloc(quo.len);
	quo.sign = z1.sign != z2.sign;
	zclearval(quo);

	ztmp3.v = zalloctemp(z2.len + 1);

	/*
	 * Normalize z1 and z2
	 */
	zshiftl(ztmp1, j);
	zshiftl(ztmp2, j);

	k = ztmp1.len - ztmp2.len;
	q = quo.v + quo.len;
	y = ztmp1.len - 1;
	h2 = ztmp2.v [ztmp2.len - 1];
	v2 = 0;
	if (ztmp2.len >= 2)
		v2 = ztmp2.v [ztmp2.len - 2];
	for (;k--; --y) {
		pp = ztmp1.v + y - 1;
		pair.silow = pp[0];
		pair.sihigh = pp[1];
		if (ztmp1.v[y] == h2)
			x = BASE1;
		else
			x = pair.ivalue / h2;
		if (x) {
			while (pair.ivalue - x * h2 < BASE && y > 1 &&
				v2 * x > (pair.ivalue - x * h2) * BASE + ztmp1.v [y-2]) {
					--x;
			}
			dmul(ztmp2, x, &ztmp3);
			if (dsub(ztmp1, ztmp3, y, ztmp2.len)) {
				--x;
				dadd(ztmp1, ztmp2, y, ztmp2.len);
			}
		}
		ztrim(&ztmp1);
		*--q = (HALF)x;
	}
	zfree(ztmp1);
	zfree(ztmp2);
	ztrim(&quo);
	*res = quo;
}


/*
 * Compute the remainder after dividing one number by another.
 * This is only defined for positive z2 values.
 * The result is normalized to lie in the range 0 to z2-1.
 */
void
zmod(z1, z2, rem)
	ZVALUE z1, z2, *rem;
{
	long i, j, k, neg;
	register HALF *pp;
	SIUNION pair;			/* pair of halfword values */
	HALF h2, v2;
	long y;
	FULL x;
	ZVALUE ztmp1, ztmp2, ztmp3;

	if (ziszero(z2))
		math_error("Division by zero");
	if (zisneg(z2))
		math_error("Non-positive modulus");
	if (ziszero(z1) || zisunit(z2)) {
		*rem = _zero_;
		return;
	}
	if (zistwo(z2)) {
		if (zisodd(z1))
			*rem = _one_;
		else
			*rem = _zero_;
		return;
	}
	neg = z1.sign;
	z1.sign = 0;

	/*
	 * Do a quick check to see if the absolute value of the number
	 * is less than the modulus.  If so, then the result is just a
	 * subtract or a copy.
	 */
	h2 = z1.v[z1.len - 1];
	v2 = z2.v[z2.len - 1];
	if ((z1.len < z2.len) || ((z1.len == z2.len) && (h2 < v2))) {
		if (neg)
			zsub(z2, z1, rem);
		else
			zcopy(z1, rem);
		return;
	}

	/*
	 * Do another quick check to see if the number is positive and
	 * between the size of the modulus and twice the modulus.
	 * If so, then the answer is just another subtract.
	 */
	if (!neg && (z1.len == z2.len) && (h2 > v2) &&
		(((FULL) h2) < 2 * ((FULL) v2)))
	{
		zsub(z1, z2, rem);
		return;
	}

	/*
	 * If the modulus is an exact power of two, then the result
	 * can be obtained by ignoring the high bits of the number.
	 * This truncation assumes that the number of words for the
	 * number is at least as large as the number of words in the
	 * modulus, which is true at this point.
	 */
	if (((v2 & -v2) == v2) && zisonebit(z2)) {	/* ASSUMES 2'S COMP */
		i = zhighbit(z2);
		z1.len = (i + BASEB - 1) / BASEB;
		zcopy(z1, &ztmp1);
		i %= BASEB;
		if (i)
			ztmp1.v[ztmp1.len - 1] &= ((((HALF) 1) << i) - 1);
		ztmp2.len = 0;
		goto gotanswer;
	}

	/*
	 * If the modulus is one less than an exact power of two, then
	 * the result can be simplified similarly to "casting out 9's".
	 * Only do this simplification for large enough modulos.
	 */
	if ((z2.len > 1) && (z2.v[0] == BASE1) && zisallbits(z2)) {
		i = -(zhighbit(z2) + 1);
		zcopy(z1, &ztmp1);
		z1 = ztmp1;
		while ((k = zrel(z1, z2)) > 0) {
			ztmp1 = _zero_;
			while (!ziszero(z1)) {
				zand(z1, z2, &ztmp2);
				zadd(ztmp2, ztmp1, &ztmp3);
				zfree(ztmp1);
				zfree(ztmp2);
				ztmp1 = ztmp3;
				zshift(z1, i, &ztmp2);
				zfree(z1);
				z1 = ztmp2;
			}
			zfree(z1);
			z1 = ztmp1;
		}
		if (k == 0) {
			zfree(ztmp1);
			*rem = _zero_;
			return;
		}
		ztmp2.len = 0;
		goto gotanswer;
	}

	/*
	 * Must actually do the divide.
	 */
	i = BASE / 2;
	j = 0;
	k = (long) z2.v[z2.len - 1];
	while (! (k & i)) {
		j ++;
		i >>= 1;
	}
	ztmp1.v = alloc(z1.len + 1);
	ztmp1.len = z1.len + 1;
	zcopyval(z1, ztmp1);
	ztmp1.v[z1.len] = 0;
	ztmp1.sign = 0;
	ztmp2.v = alloc(z2.len);
	ztmp2.len = z2.len;
	ztmp2.sign = 0;
	zcopyval(z2, ztmp2);
	if (zrel(ztmp1, ztmp2) < 0)
		goto gotanswer;

	ztmp3.v = zalloctemp(z2.len + 1);

	/*
	 * Normalize z1 and z2
	 */
	zshiftl(ztmp1, j);
	zshiftl(ztmp2, j);

	k = ztmp1.len - ztmp2.len;
	y = ztmp1.len - 1;
	h2 = ztmp2.v [ztmp2.len - 1];
	v2 = 0;
	if (ztmp2.len >= 2)
		v2 = ztmp2.v [ztmp2.len - 2];
	for (;k--; --y) {
		pp = ztmp1.v + y - 1;
		pair.silow = pp[0];
		pair.sihigh = pp[1];
		if (ztmp1.v[y] == h2)
			x = BASE1;
		else
			x = pair.ivalue / h2;
		if (x) {
			while (pair.ivalue - x * h2 < BASE && y > 1 &&
				v2 * x > (pair.ivalue - x * h2) * BASE + ztmp1.v [y-2]) {
					--x;
			}
			dmul(ztmp2, x, &ztmp3);
			if (dsub(ztmp1, ztmp3, y, ztmp2.len))
				dadd(ztmp1, ztmp2, y, ztmp2.len);
		}
		ztrim(&ztmp1);
	}
	zshiftr(ztmp1, j);

gotanswer:
	ztrim(&ztmp1);
	if (ztmp2.len)
		zfree(ztmp2);
	if (neg && !ziszero(ztmp1)) {