qmath.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 rational arithmetic primitive routines
 */

#include "qmath.h"


NUMBER _qzero_ =	{ { _zeroval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
NUMBER _qone_ =		{ { _oneval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
static NUMBER _qtwo_ =	{ { _twoval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
static NUMBER _qten_ =	{ { _tenval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
NUMBER _qnegone_ =	{ { _oneval_, 1, 1 }, { _oneval_, 1, 0 }, 1 };
NUMBER _qonehalf_ =	{ { _oneval_, 1, 0 }, { _twoval_, 1, 0 }, 1 };


/*
 * Create another copy of a number.
 *	q2 = qcopy(q1);
 */
NUMBER *
qcopy(q)
	register NUMBER *q;
{
	register NUMBER *r;

	r = qalloc();
	r->num.sign = q->num.sign;
	if (!zisunit(q->num)) {
		r->num.len = q->num.len;
		r->num.v = alloc(r->num.len);
		zcopyval(q->num, r->num);
	}
	if (!zisunit(q->den)) {
		r->den.len = q->den.len;
		r->den.v = alloc(r->den.len);
		zcopyval(q->den, r->den);
	}
	return r;
}


/*
 * Convert a number to a normal integer.
 *	i = qtoi(q);
 */
long
qtoi(q)
	register NUMBER *q;
{
	long i;
	ZVALUE res;

	if (qisint(q))
		return ztoi(q->num);
	zquo(q->num, q->den, &res);
	i = ztoi(res);
	zfree(res);
	return i;
}


/*
 * Convert a normal integer into a number.
 *	q = itoq(i);
 */
NUMBER *
itoq(i)
	long i;
{
	register NUMBER *q;

	if ((i >= -1) && (i <= 10)) {
		switch ((int) i) {
			case 0: q = &_qzero_; break;
			case 1: q = &_qone_; break;
			case 2: q = &_qtwo_; break;
			case 10: q = &_qten_; break;
			case -1: q = &_qnegone_; break;
			default: q = NULL;
		}
		if (q)
			return qlink(q);
	}
	q = qalloc();
	itoz(i, &q->num);
	return q;
}


/*
 * Create a number from the given integral numerator and denominator.
 *	q = iitoq(inum, iden);
 */
NUMBER *
iitoq(inum, iden)
	long inum, iden;
{
	register NUMBER *q;
	long d;
	BOOL sign;

	if (iden == 0)
		math_error("Division by zero");
	if (inum == 0)
		return qlink(&_qzero_);
	sign = 0;
	if (inum < 0) {
		sign = 1;
		inum = -inum;
	}
	if (iden < 0) {
		sign = 1 - sign;
		iden = -iden;
	}
	d = iigcd(inum, iden);
	inum /= d;
	iden /= d;
	if (iden == 1)
		return itoq(sign ? -inum : inum);
	q = qalloc();
	if (inum != 1)
		itoz(inum, &q->num);
	itoz(iden, &q->den);
	q->num.sign = sign;
	return q;
}


/*
 * Add two numbers to each other.
 *	q3 = qadd(q1, q2);
 */
NUMBER *
qadd(q1, q2)
	register NUMBER *q1, *q2;
{
	NUMBER *r;
	ZVALUE t1, t2, temp, d1, d2, vpd1, upd1;

	if (qiszero(q1))
		return qlink(q2);
	if (qiszero(q2))
		return qlink(q1);
	r = qalloc();
	/*
	 * If either number is an integer, then the result is easy.
	 */
	if (qisint(q1) && qisint(q2)) {
		zadd(q1->num, q2->num, &r->num);
		return r;
	}
	if (qisint(q2)) {
		zmul(q1->den, q2->num, &temp);
		zadd(q1->num, temp, &r->num);
		zfree(temp);
		zcopy(q1->den, &r->den);
		return r;
	}
	if (qisint(q1)) {
		zmul(q2->den, q1->num, &temp);
		zadd(q2->num, temp, &r->num);
		zfree(temp);
		zcopy(q2->den, &r->den);
		return r;
	}
	/*
	 * Both arguments are true fractions, so we need more work.
	 * If the denominators are relatively prime, then the answer is the
	 * straightforward cross product result with no need for reduction.
	 */
	zgcd(q1->den, q2->den, &d1);
	if (zisunit(d1)) {
		zfree(d1);
		zmul(q1->num, q2->den, &t1);
		zmul(q1->den, q2->num, &t2);
		zadd(t1, t2, &r->num);
		zfree(t1);
		zfree(t2);
		zmul(q1->den, q2->den, &r->den);
		return r;
	}
	/*
	 * The calculation is now more complicated.
	 * See Knuth Vol 2 for details.
	 */
	zquo(q2->den, d1, &vpd1);
	zquo(q1->den, d1, &upd1);
	zmul(q1->num, vpd1, &t1);
	zmul(q2->num, upd1, &t2);
	zadd(t1, t2, &temp);
	zfree(t1);
	zfree(t2);
	zfree(vpd1);
	zgcd(temp, d1, &d2);
	zfree(d1);
	if (zisunit(d2)) {
		zfree(d2);
		r->num = temp;
		zmul(upd1, q2->den, &r->den);
		zfree(upd1);
		return r;
	}
	zquo(temp, d2, &r->num);
	zfree(temp);
	zquo(q2->den, d2, &temp);
	zfree(d2);
	zmul(temp, upd1, &r->den);
	zfree(temp);
	zfree(upd1);
	return r;
}


/*
 * Subtract one number from another.
 *	q3 = qsub(q1, q2);
 */
NUMBER *
qsub(q1, q2)
	register NUMBER *q1, *q2;
{
	NUMBER *r;

	if (q1 == q2)
		return qlink(&_qzero_);
	if (qiszero(q2))
		return qlink(q1);
	if (qisint(q1) && qisint(q2)) {
		r = qalloc();
		zsub(q1->num, q2->num, &r->num);
		return r;
	}
	q2 = qneg(q2);
	if (qiszero(q1))
		return q2;
	r = qadd(q1, q2);
	qfree(q2);
	return r;
}


/*
 * Increment a number by one.
 */
NUMBER *
qinc(q)
	NUMBER *q;
{
	NUMBER *r;

	r = qalloc();
	if (qisint(q)) {
		zadd(q->num, _one_, &r->num);
		return r;
	}
	zadd(q->num, q->den, &r->num);
	zcopy(q->den, &r->den);
	return r;
}


/*
 * Decrement a number by one.
 */
NUMBER *
qdec(q)
	NUMBER *q;
{
	NUMBER *r;

	r = qalloc();
	if (qisint(q)) {
		zsub(q->num, _one_, &r->num);
		return r;
	}
	zsub(q->num, q->den, &r->num);
	zcopy(q->den, &r->den);
	return r;
}


/*
 * Add a normal small integer value to an arbitrary number.
 */
NUMBER *
qaddi(q1, n)
	NUMBER *q1;
	long n;
{
	NUMBER addnum;		/* temporary number */
	HALF addval[2];		/* value of small number */
	BOOL neg;		/* TRUE if number is neg */

	if (n == 0)
		return qlink(q1);
	if (n == 1)
		return qinc(q1);
	if (n == -1)
		return qdec(q1);
	if (qiszero(q1))
		return itoq(n);
	addnum.num.sign = 0;
	addnum.num.len = 1;
	addnum.num.v = addval;
	addnum.den = _one_;
	neg = (n < 0);
	if (neg)
		n = -n;
	addval[0] = (HALF) n;
	n = (((FULL) n) >> BASEB);
	if (n) {
		addval[1] = (HALF) n;
		addnum.num.len = 2;
	}
	if (neg)
		return qsub(q1, &addnum);
	else
		return qadd(q1, &addnum);
}


/*
 * Multiply two numbers.
 *	q3 = qmul(q1, q2);
 */
NUMBER *
qmul(q1, q2)
	register NUMBER *q1, *q2;
{
	NUMBER *r;			/* returned value */
	ZVALUE n1, n2, d1, d2;		/* numerators and denominators */
	ZVALUE tmp;

	if (qiszero(q1) || qiszero(q2))
		return qlink(&_qzero_);
	if (qisone(q1))
		return qlink(q2);
	if (qisone(q2))
		return qlink(q1);
	if (qisint(q1) && qisint(q2)) {	/* easy results if integers */
		r = qalloc();
		zmul(q1->num, q2->num, &r->num);
		return r;
	}
	n1 = q1->num;
	n2 = q2->num;
	d1 = q1->den;
	d2 = q2->den;
	if (ziszero(d1) || ziszero(d2))
		math_error("Division by zero");
	if (ziszero(n1) || ziszero(n2))
		return qlink(&_qzero_);
	if (!zisunit(n1) && !zisunit(d2)) {	/* possibly reduce */
		zgcd(n1, d2, &tmp);
		if (!zisunit(tmp)) {
			zquo(q1->num, tmp, &n1);
			zquo(q2->den, tmp, &d2);
		}
		zfree(tmp);
	}
	if (!zisunit(n2) && !zisunit(d1)) {	/* again possibly reduce */
		zgcd(n2, d1, &tmp);
		if (!zisunit(tmp)) {
			zquo(q2->num, tmp, &n2);
			zquo(q1->den, tmp, &d1);
		}
		zfree(tmp);
	}
	r = qalloc();
	zmul(n1, n2, &r->num);
	zmul(d1, d2, &r->den);
	if (q1->num.v != n1.v)
		zfree(n1);
	if (q1->den.v != d1.v)
		zfree(d1);
	if (q2->num.v != n2.v)
		zfree(n2);
	if (q2->den.v != d2.v)
		zfree(d2);
	return r;
}


/*
 * Multiply a number by a small integer.
 *	q2 = qmuli(q1, n);
 */
NUMBER *
qmuli(q, n)
	NUMBER *q;
	long n;
{
	NUMBER *r;
	long d;			/* gcd of multiplier and denominator */
	int sign;

	if ((n == 0) || qiszero(q))
		return qlink(&_qzero_);
	if (n == 1)
		return qlink(q);
	r = qalloc();
	if (qisint(q)) {
		zmuli(q->num, n, &r->num);
		return r;
	}
	sign = 1;
	if (n < 0) {
		n = -n;
		sign = -1;
	}
	d = zmodi(q->den, n);
	d = iigcd(d, n);
	zmuli(q->num, (n * sign) / d, &r->num);
	(void) zdivi(q->den, d, &r->den);
	return r;
}


/*
 * Divide two numbers (as fractions).
 *	q3 = qdiv(q1, q2);
 */
NUMBER *
qdiv(q1, q2)
	register NUMBER *q1, *q2;
{
	NUMBER temp;

	if (qiszero(q2))
		math_error("Division by zero");
	if ((q1 == q2) || !qcmp(q1, q2))
		return qlink(&_qone_);
	if (qisone(q1))
		return qinv(q2);
	temp.num = q2->den;
	temp.den = q2->num;
	temp.num.sign = temp.den.sign;
	temp.den.sign = 0;
	temp.links = 1;
	return qmul(q1, &temp);
}


/*
 * Divide a number by a small integer.
 *	q2 = qdivi(q1, n);
 */
NUMBER *
qdivi(q, n)
	NUMBER *q;
	long n;
{
	NUMBER *r;
	long d;			/* gcd of divisor and numerator */
	int sign;

	if (n == 0)
		math_error("Division by zero");
	if ((n == 1) || qiszero(q))
		return qlink(q);
	sign = 1;
	if (n < 0) {
		n = -n;
		sign = -1;
	}
	r = qalloc();
	d = zmodi(q->num, n);
	d = iigcd(d, n);
	(void) zdivi(q->num, d * sign, &r->num);
	zmuli(q->den, n / d, &r->den);
	return r;
}


/*
 * Return the quotient when one number is divided by another.
 * This works for fractional values also, and in all cases:
 *	qquo(q1, q2) = int(q1 / q2).
 */
NUMBER *
qquo(q1, q2)
	register NUMBER *q1, *q2;
{
	ZVALUE num, den, res;
	NUMBER *q;

	if (zisunit(q1->num))
		num = q2->den;
	else if (zisunit(q2->den))
		num = q1->num;
	else
		zmul(q1->num, q2->den, &num);
	if (zisunit(q1->den))
		den = q2->num;
	else if (zisunit(q2->num))
		den = q1->den;
	else
		zmul(q1->den, q2->num, &den);
	zquo(num, den, &res);
	if ((num.v != q2->den.v) && (num.v != q1->num.v))
		zfree(num);
	if ((den.v != q2->num.v) && (den.v != q1->den.v))
		zfree(den);
	if (ziszero(res)) {
		zfree(res);
		return qlink(&_qzero_);
	}
	res.sign = (q1->num.sign != q2->num.sign);
	if (zisunit(res)) {
		q = (res.sign ? &_qnegone_ : &_qone_);
		zfree(res);
		return qlink(q);
	}
	q = qalloc();
	q->num = res;
	return q;
}


/*
 * Return the absolute value of a number.
 *	q2 = qabs(q1);
 */
NUMBER *
qabs(q)
	register NUMBER *q;
{
	register NUMBER *r;

	if (q->num.sign == 0)
		return qlink(q);
	r = qalloc();
	if (!zisunit(q->num))
		zcopy(q->num, &r->num);
	if (!zisunit(q->den))
		zcopy(q->den, &r->den);
	r->num.sign = 0;
	return r;
}


/*
 * Negate a number.
 *	q2 = qneg(q1);
 */
NUMBER *
qneg(q)
	register NUMBER *q;
{
	register NUMBER *r;

	if (qiszero(q))
		return qlink(&_qzero_);
	r = qalloc();
	if (!zisunit(q->num))
		zcopy(q->num, &r->num);
	if (!zisunit(q->den))
		zcopy(q->den, &r->den);
	r->num.sign = !q->num.sign;
	return r;
}


/*
 * Return the sign of a number (-1, 0 or 1)
 */
NUMBER *
qsign(q)
	NUMBER *q;
{
	if (qiszero(q))
		return qlink(&_qzero_);
	if (qisneg(q))
		return qlink(&_qnegone_);
	return qlink(&_qone_);
}


/*
 * Invert a number.
 *	q2 = qinv(q1);
 */
NUMBER *
qinv(q)
	register NUMBER *q;
{
	register NUMBER *r;

	if (qisunit(q)) {
		r = (qisneg(q) ? &_qnegone_ : &_qone_);
		return qlink(r);
	}
	if (qiszero(q))
		math_error("Division by zero");
	r = qalloc();
	if (!zisunit(q->num))
		zcopy(q->num, &r->den);
	if (!zisunit(q->den))
		zcopy(q->den, &r->num);
	r->num.sign = q->num.sign;
	r->den.sign = 0;
	return r;