qmath.c

早期freebsd实现

}


/*
 * Return just the numerator of a number.
 *	q2 = qnum(q1);
 */
NUMBER *
qnum(q)
	register NUMBER *q;
{
	register NUMBER *r;

	if (qisint(q))
		return qlink(q);
	if (zisunit(q->num)) {
		r = (qisneg(q) ? &_qnegone_ : &_qone_);
		return qlink(r);
	}
	r = qalloc();
	zcopy(q->num, &r->num);
	return r;
}


/*
 * Return just the denominator of a number.
 *	q2 = qden(q1);
 */
NUMBER *
qden(q)
	register NUMBER *q;
{
	register NUMBER *r;

	if (qisint(q))
		return qlink(&_qone_);
	r = qalloc();
	zcopy(q->den, &r->num);
	return r;
}


/*
 * Return the fractional part of a number.
 *	q2 = qfrac(q1);
 */
NUMBER *
qfrac(q)
	register NUMBER *q;
{
	register NUMBER *r;
	ZVALUE z;

	if (qisint(q))
		return qlink(&_qzero_);
	if ((q->num.len < q->den.len) || ((q->num.len == q->den.len) &&
		(q->num.v[q->num.len - 1] < q->den.v[q->num.len - 1])))
			return qlink(q);
	r = qalloc();
	if (qisneg(q)) {
		zmod(q->num, q->den, &z);
		zsub(q->den, z, &r->num);
		zfree(z);
	} else {
		zmod(q->num, q->den, &r->num);
	}
	zcopy(q->den, &r->den);
	r->num.sign = q->num.sign;
	return r;
}


/*
 * Return the integral part of a number.
 *	q2 = qint(q1);
 */
NUMBER *
qint(q)
	register NUMBER *q;
{
	register NUMBER *r;

	if (qisint(q))
		return qlink(q);
	if ((q->num.len < q->den.len) || ((q->num.len == q->den.len) &&
		(q->num.v[q->num.len - 1] < q->den.v[q->num.len - 1])))
			return qlink(&_qzero_);
	r = qalloc();
	zquo(q->num, q->den, &r->num);
	return r;
}


/*
 * Compute the square of a number.
 */
NUMBER *
qsquare(q)
	register NUMBER *q;
{
	ZVALUE num, den;

	if (qiszero(q))
		return qlink(&_qzero_);
	if (qisunit(q))
		return qlink(&_qone_);
	num = q->num;
	den = q->den;
	q = qalloc();
	if (!zisunit(num))
		zsquare(num, &q->num);
	if (!zisunit(den))
		zsquare(den, &q->den);
	return q;
}


/*
 * Shift an integer by a given number of bits. This multiplies the number
 * by the appropriate power of two.  Positive numbers shift left, negative
 * ones shift right.  Low bits are truncated when shifting right.
 */
NUMBER *
qshift(q, n)
	NUMBER *q;
	long n;
{
	register NUMBER *r;

	if (qisfrac(q))
		math_error("Shift of non-integer");
	if (qiszero(q) || (n == 0))
		return qlink(q);
	if (n <= -(q->num.len * BASEB))
		return qlink(&_qzero_);
	r = qalloc();
	zshift(q->num, n, &r->num);
	return r;
}


/*
 * Scale a number by a power of two, as in:
 *	ans = q * 2^n.
 * This is similar to shifting, except that fractions work.
 */
NUMBER *
qscale(q, pow)
	NUMBER *q;
	long pow;
{
	long numshift, denshift, tmp;
	NUMBER *r;

	if (qiszero(q) || (pow == 0))
		return qlink(q);
	if ((pow > 1000000L) || (pow < -1000000L))
		math_error("Very large scale value");
	numshift = zisodd(q->num) ? 0 : zlowbit(q->num);
	denshift = zisodd(q->den) ? 0 : zlowbit(q->den);
	if (pow > 0) {
		tmp = pow;
		if (tmp > denshift)
		tmp = denshift;
		denshift = -tmp;
		numshift = (pow - tmp);
	} else {
		pow = -pow;
		tmp = pow;
		if (tmp > numshift)
			tmp = numshift;
		numshift = -tmp;
		denshift = (pow - tmp);
	}
	r = qalloc();
	if (numshift)
		zshift(q->num, numshift, &r->num);
	else
		zcopy(q->num, &r->num);
	if (denshift)
		zshift(q->den, denshift, &r->den);
	else
		zcopy(q->den, &r->den);
	return r;
}


/*
 * Return the minimum of two numbers.
 */
NUMBER *
qmin(q1, q2)
	NUMBER *q1, *q2;
{
	if (q1 == q2)
		return qlink(q1);
	if (qrel(q1, q2) > 0)
		q1 = q2;
	return qlink(q1);
}


/*
 * Return the maximum of two numbers.
 */
NUMBER *
qmax(q1, q2)
	NUMBER *q1, *q2;
{
	if (q1 == q2)
		return qlink(q1);
	if (qrel(q1, q2) < 0)
		q1 = q2;
	return qlink(q1);
}


/*
 * Perform the logical OR of two integers.
 */
NUMBER *
qor(q1, q2)
	NUMBER *q1, *q2;
{
	register NUMBER *r;

	if (qisfrac(q1) || qisfrac(q2))
		math_error("Non-integers for logical or");
	if ((q1 == q2) || qiszero(q2))
		return qlink(q1);
	if (qiszero(q1))
		return qlink(q2);
	r = qalloc();
	zor(q1->num, q2->num, &r->num);
	return r;
}


/*
 * Perform the logical AND of two integers.
 */
NUMBER *
qand(q1, q2)
	NUMBER *q1, *q2;
{
	register NUMBER *r;
	ZVALUE res;

	if (qisfrac(q1) || qisfrac(q2))
		math_error("Non-integers for logical and");
	if (q1 == q2)
		return qlink(q1);
	if (qiszero(q1) || qiszero(q2))
		return qlink(&_qzero_);
	zand(q1->num, q2->num, &res);
	if (ziszero(res)) {
		zfree(res);
		return qlink(&_qzero_);
	}
	r = qalloc();
	r->num = res;
	return r;
}


/*
 * Perform the logical XOR of two integers.
 */
NUMBER *
qxor(q1, q2)
	NUMBER *q1, *q2;
{
	register NUMBER *r;
	ZVALUE res;

	if (qisfrac(q1) || qisfrac(q2))
		math_error("Non-integers for logical xor");
	if (q1 == q2)
		return qlink(&_qzero_);
	if (qiszero(q1))
		return qlink(q2);
	if (qiszero(q2))
		return qlink(q1);
	zxor(q1->num, q2->num, &res);
	if (ziszero(res)) {
		zfree(res);
		return qlink(&_qzero_);
	}
	r = qalloc();
	r->num = res;
	return r;
}


#if 0
/*
 * Return the number whose binary representation only has the specified
 * bit set (counting from zero).  This thus produces a given power of two.
 */
NUMBER *
qbitvalue(n)
	long n;
{
	register NUMBER *r;

	if (n <= 0)
		return qlink(&_qone_);
	r = qalloc();
	zbitvalue(n, &r->num);
	return r;
}


/*
 * Test to see if the specified bit of a number is on (counted from zero).
 * Returns TRUE if the bit is set, or FALSE if it is not.
 *	i = qbittest(q, n);
 */
BOOL
qbittest(q, n)
	register NUMBER *q;
	long n;
{
	int x, y;

	if ((n < 0) || (n >= (q->num.len * BASEB)))
		return FALSE;
	x = q->num.v[n / BASEB];
	y = (1 << (n % BASEB));
	return ((x & y) != 0);
}
#endif


/*
 * Return the precision of a number (usually for examining an epsilon value).
 * This is the largest power of two whose reciprocal is not smaller in absolute
 * value than the specified number.  For example, qbitprec(1/100) = 6.
 * Numbers larger than one have a precision of zero.
 */
long
qprecision(q)
	NUMBER *q;
{
	long r;

	if (qisint(q))
		return 0;
	if (zisunit(q->num))
		return zhighbit(q->den);
	r = zhighbit(q->den) - zhighbit(q->num) - 1;
	if (r < 0)
		r = 0;
	return r;
}


#if 0
/*
 * Return an integer indicating the sign of a number (-1, 0, or 1).
 *	i = qtst(q);
 */
FLAG
qtest(q)
	register NUMBER *q;
{
	if (!ztest(q->num))
		return 0;
	if (q->num.sign)
		return -1;
	return 1;
}
#endif


/*
 * Determine whether or not one number exactly divides another one.
 * Returns TRUE if the first number is an integer multiple of the second one.
 */
BOOL
qdivides(q1, q2)
	NUMBER *q1, *q2;
{
	if (qiszero(q1))
		return TRUE;
	if (qisint(q1) && qisint(q2)) {
		if (qisunit(q2))
			return TRUE;
		return zdivides(q1->num, q2->num);
	}
	return zdivides(q1->num, q2->num) && zdivides(q2->den, q1->den);
}


/*
 * Compare two numbers and return an integer indicating their relative size.
 *	i = qrel(q1, q2);
 */
FLAG
qrel(q1, q2)
	register NUMBER *q1, *q2;
{
	ZVALUE z1, z2;
	long wc1, wc2;
	int sign;
	int z1f = 0, z2f = 0;

	if (q1 == q2)
		return 0;
	sign = q2->num.sign - q1->num.sign;
	if (sign)
		return sign;
	if (qiszero(q2))
		return !qiszero(q1);
	if (qiszero(q1))
		return -1;
	/*
	 * Make a quick comparison by calculating the number of words resulting as
	 * if we multiplied through by the denominators, and then comparing the
	 * word counts.
	 */
	sign = 1;
	if (qisneg(q1))
		sign = -1;
	wc1 = q1->num.len + q2->den.len;
	wc2 = q2->num.len + q1->den.len;
	if (wc1 < wc2 - 1)
		return -sign;
	if (wc2 < wc1 - 1)
		return sign;
	/*
	 * Quick check failed, must actually do the full comparison.
	 */
	if (zisunit(q2->den))
		z1 = q1->num;
	else if (zisone(q1->num))
		z1 = q2->den;
	else {
		z1f = 1;
		zmul(q1->num, q2->den, &z1);
	}
	if (zisunit(q1->den))
		z2 = q2->num;
	else if (zisone(q2->num))
		z2 = q1->den;
	else {
		z2f = 1;
		zmul(q2->num, q1->den, &z2);
	}
	sign = zrel(z1, z2);
	if (z1f)
		zfree(z1);
	if (z2f)
		zfree(z2);
	return sign;
}


/*
 * Compare two numbers to see if they are equal.
 * This differs from qrel in that the numbers are not ordered.
 * Returns TRUE if they differ.
 */
BOOL
qcmp(q1, q2)
	register NUMBER *q1, *q2;
{
	if (q1 == q2)
		return FALSE;
	if ((q1->num.sign != q2->num.sign) || (q1->num.len != q2->num.len) ||
		(q2->den.len != q2->den.len) || (*q1->num.v != *q2->num.v) ||
		(*q1->den.v != *q2->den.v))
			return TRUE;
	if (zcmp(q1->num, q2->num))
		return TRUE;
	if (qisint(q1))
		return FALSE;
	return zcmp(q1->den, q2->den);
}


/*
 * Compare a number against a normal small integer.
 * Returns 1, 0, or -1, according to whether the first number is greater,
 * equal, or less than the second number.
 *	n = qreli(q, n);
 */
FLAG
qreli(q, n)
	NUMBER *q;
	long n;
{
	int sign;
	ZVALUE num;
	HALF h2[2];
	NUMBER q2;

	sign = ztest(q->num);		/* do trivial sign checks */
	if (sign == 0) {
		if (n > 0)
			return -1;
		return (n < 0);
	}
	if ((sign < 0) && (n >= 0))
		return -1;
	if ((sign > 0) && (n <= 0))
		return 1;
	n *= sign;
	if (n == 1) {			/* quick check against 1 or -1 */
		num = q->num;
		num.sign = 0;
		return (sign * zrel(num, q->den));
	}
	num.sign = (sign < 0);
	num.len = 1 + (n >= BASE);
	num.v = h2;
	h2[0] = (n & BASE1);
	h2[1] = (n >> BASEB);
	if (zisunit(q->den))	/* integer compare if no denominator */
		return zrel(q->num, num);
	q2.num = num;
	q2.den = _one_;
	q2.links = 1;
	return qrel(q, &q2);	/* full fractional compare */
}


/*
 * Compare a number against a small integer to see if they are equal.
 * Returns TRUE if they differ.
 */
BOOL
qcmpi(q, n)
	NUMBER *q;
	long n;
{
	long len;

	len = q->num.len;
	if ((len > 2) || qisfrac(q) || (q->num.sign != (n < 0)))
		return TRUE;
	if (n < 0)
		n = -n;
	if (((HALF)(n)) != q->num.v[0])
		return TRUE;
	n = ((FULL) n) >> BASEB;
	return (((n != 0) != (len == 2)) || (n != q->num.v[1]));
}


/*
 * Number node allocation routines
 */

#define	NNALLOC	1000

union allocNode {
	NUMBER	num;
	union allocNode	*link;
};

static union allocNode	*freeNum;


NUMBER *
qalloc()
{
	register union allocNode *temp;

	if (freeNum == NULL) {
		freeNum = (union allocNode *)
		malloc(sizeof (NUMBER) * NNALLOC);
		if (freeNum == NULL)
			math_error("Not enough memory");
		freeNum[NNALLOC-1].link = NULL;
		for (temp=freeNum+NNALLOC-2; temp >= freeNum; --temp) {
			temp->link = temp+1;
		}
	}
	temp = freeNum;
	freeNum = temp->link;
	temp->num.links = 1;
	temp->num.num = _one_;
	temp->num.den = _one_;
	return &temp->num;
}


void
qfreenum(q)
	register NUMBER *q;
{
	union allocNode *a;

	if (q == NULL)
		return;
	zfree(q->num);
	zfree(q->den);
	a = (union allocNode *) q;
	a->link = freeNum;
	freeNum = a;
}

/* END CODE */