qfunc.c

早期freebsd实现

	register NUMBER *q1, *q2;
{
	register NUMBER *r;

	if (qisfrac(q1) || qisfrac(q2))
		math_error("Non-integral arguments for combinatorial");
	r = qalloc();
	zcomb(q1->num, q2->num, &r->num);
	return r;
}


/*
 * Compute the Jacobi function (a / b).
 * -1 => a is not quadratic residue mod b
 * 1 => b is composite, or a is quad residue of b
 * 0 => b is even or b < 0
 */
NUMBER *
qjacobi(q1, q2)
	register NUMBER *q1, *q2;
{
	if (qisfrac(q1) || qisfrac(q2))
		math_error("Non-integral arguments for jacobi");
	return itoq((long) zjacobi(q1->num, q2->num));
}


/*
 * Compute the Fibonacci number F(n).
 */
NUMBER *
qfib(q)
	register NUMBER *q;
{
	register NUMBER *r;

	if (qisfrac(q))
		math_error("Non-integral Fibonacci number");
	r = qalloc();
	zfib(q->num, &r->num);
	return r;
}


/*
 * Truncate a number to the specified number of decimal places.
 * Specifying zero places makes the result identical to qint.
 * Example: qtrunc(2/3, 3) = .666
 */
NUMBER *
qtrunc(q1, q2)
	NUMBER *q1, *q2;
{
	long places;
	NUMBER *r;
	ZVALUE tenpow, tmp1, tmp2;

	if (qisfrac(q2) || qisneg(q2) || !zistiny(q2->num))
		math_error("Bad number of places for qtrunc");
	if (qisint(q1))
		return qlink(q1);
	r = qalloc();
	places = z1tol(q2->num);
	/*
	 * Ok, produce the result.
	 * First see if we want no places, in which case just take integer part.
	 */
	if (places == 0) {
		zquo(q1->num, q1->den, &r->num);
		return r;
	}
	ztenpow(places, &tenpow);
	zmul(q1->num, tenpow, &tmp1);
	zquo(tmp1, q1->den, &tmp2);
	zfree(tmp1);
	if (ziszero(tmp2)) {
		zfree(tmp2);
		return qlink(&_qzero_);
	}
	/*
	 * Now reduce the result to the lowest common denominator.
	 */
	zgcd(tmp2, tenpow, &tmp1);
	if (zisunit(tmp1)) {
		zfree(tmp1);
		r->num = tmp2;
		r->den = tenpow;
		return r;
	}
	zquo(tmp2, tmp1, &r->num);
	zquo(tenpow, tmp1, &r->den);
	zfree(tmp1);
	zfree(tmp2);
	zfree(tenpow);
	return r;
}


/*
 * Round a number to the specified number of decimal places.
 * Zero decimal places means round to the nearest integer.
 * Example: qround(2/3, 3) = .667
 */
NUMBER *
qround(q, places)
	NUMBER *q;		/* number to be rounded */
	long places;		/* number of decimal places to round to */
{
	NUMBER *r;
	ZVALUE tenpow, roundval, tmp1, tmp2;

	if (places < 0)
		math_error("Negative places for qround");
	if (qisint(q))
		return qlink(q);
	/*
	 * Calculate one half of the denominator, ignoring fractional results.
	 * This is the value we will add in order to cause rounding.
	 */
	zshift(q->den, -1L, &roundval);
	roundval.sign = q->num.sign;
	/*
	 * Ok, now do the actual work to produce the result.
	 */
	r = qalloc();
	ztenpow(places, &tenpow);
	zmul(q->num, tenpow, &tmp2);
	zadd(tmp2, roundval, &tmp1);
	zfree(tmp2);
	zfree(roundval);
	zquo(tmp1, q->den, &tmp2);
	zfree(tmp1);
	if (ziszero(tmp2)) {
		zfree(tmp2);
		return qlink(&_qzero_);
	}
	/*
	 * Now reduce the result to the lowest common denominator.
	 */
	zgcd(tmp2, tenpow, &tmp1);
	if (zisunit(tmp1)) {
		zfree(tmp1);
		r->num = tmp2;
		r->den = tenpow;
		return r;
	}
	zquo(tmp2, tmp1, &r->num);
	zquo(tenpow, tmp1, &r->den);
	zfree(tmp1);
	zfree(tmp2);
	zfree(tenpow);
	return r;
}


/*
 * Truncate a number to the specified number of binary places.
 * Specifying zero places makes the result identical to qint.
 */
NUMBER *
qbtrunc(q1, q2)
	NUMBER *q1, *q2;
{
	long places, twopow;
	NUMBER *r;
	ZVALUE tmp1, tmp2;

	if (qisfrac(q2) || qisneg(q2) || !zistiny(q2->num))
		math_error("Bad number of places for qtrunc");
	if (qisint(q1))
		return qlink(q1);
	r = qalloc();
	places = z1tol(q2->num);
	/*
	 * Ok, produce the result.
	 * First see if we want no places, in which case just take integer part.
	 */
	if (places == 0) {
		zquo(q1->num, q1->den, &r->num);
		return r;
	}
	zshift(q1->num, places, &tmp1);
	zquo(tmp1, q1->den, &tmp2);
	zfree(tmp1);
	if (ziszero(tmp2)) {
		zfree(tmp2);
		return qlink(&_qzero_);
	}
	/*
	 * Now reduce the result to the lowest common denominator.
	 */
	if (zisodd(tmp2)) {
		r->num = tmp2;
		zbitvalue(places, &r->den);
		return r;
	}
	twopow = zlowbit(tmp2);
	if (twopow > places)
		twopow = places;
	places -= twopow;
	zshift(tmp2, -twopow, &r->num);
	zfree(tmp2);
	zbitvalue(places, &r->den);
	return r;
}


/*
 * Round a number to the specified number of binary places.
 * Zero binary places means round to the nearest integer.
 */
NUMBER *
qbround(q, places)
	NUMBER *q;		/* number to be rounded */
	long places;		/* number of binary places to round to */
{
	long twopow;
	NUMBER *r;
	ZVALUE roundval, tmp1, tmp2;

	if (places < 0)
		math_error("Negative places for qbround");
	if (qisint(q))
		return qlink(q);
	r = qalloc();
	/*
	 * Calculate one half of the denominator, ignoring fractional results.
	 * This is the value we will add in order to cause rounding.
	 */
	zshift(q->den, -1L, &roundval);
	roundval.sign = q->num.sign;
	/*
	 * Ok, produce the result.
	 */
	zshift(q->num, places, &tmp1);
	zadd(tmp1, roundval, &tmp2);
	zfree(roundval);
	zfree(tmp1);
	zquo(tmp2, q->den, &tmp1);
	zfree(tmp2);
	if (ziszero(tmp1)) {
		zfree(tmp1);
		return qlink(&_qzero_);
	}
	/*
	 * Now reduce the result to the lowest common denominator.
	 */
	if (zisodd(tmp1)) {
		r->num = tmp1;
		zbitvalue(places, &r->den);
		return r;
	}
	twopow = zlowbit(tmp1);
	if (twopow > places)
		twopow = places;
	places -= twopow;
	zshift(tmp1, -twopow, &r->num);
	zfree(tmp1);
	zbitvalue(places, &r->den);
	return r;
}


/*
 * Approximate a number by using binary rounding with the minimum number
 * of binary places so that the resulting number is within the specified
 * epsilon of the original number.
 */
NUMBER *
qbappr(q, e)
	NUMBER *q, *e;
{
	long bits;

	if (qisneg(e) || qiszero(e))
		math_error("Bad epsilon value for qbappr");
	if (e == _epsilon_)
		bits = _epsilonprec_ + 1;
	else
		bits = qprecision(e) + 1;
	return qbround(q, bits);
}


/*
 * Approximate a number using continued fractions until the approximation
 * error is less than the specified value.  If a NULL pointer is given
 * for the error value, then the closest simpler fraction is returned.
 */
NUMBER *
qcfappr(q, e)
	NUMBER *q, *e;
{
	NUMBER qtest, *qtmp;
	ZVALUE u1, u2, u3, v1, v2, v3, t1, t2, t3, qq, tt;
	int i;
	BOOL haveeps;

	haveeps = TRUE;
	if (e == NULL) {
		haveeps = FALSE;
		e = &_qzero_;
	}
	if (qisneg(e))
		math_error("Negative epsilon for cfappr");
	if (qisint(q) || zisunit(q->num) || (haveeps && qiszero(e)))
		return qlink(q);
	u1 = _one_;
	u2 = _zero_;
	u3 = q->num;
	u3.sign = 0;
	v1 = _zero_;
	v2 = _one_;
	v3 = q->den;
	while (!ziszero(v3)) {
		if (!ziszero(u2) && !ziszero(u1)) {
			qtest.num = u2;
			qtest.den = u1;
			qtest.den.sign = 0;
			qtest.num.sign = q->num.sign;
			qtmp = qsub(q, &qtest);
			qtest = *qtmp;
			qtest.num.sign = 0;
			i = qrel(&qtest, e);
			qfree(qtmp);
			if (i <= 0)
				break;
		}
		zquo(u3, v3, &qq);
		zmul(qq, v1, &tt); zsub(u1, tt, &t1); zfree(tt);
		zmul(qq, v2, &tt); zsub(u2, tt, &t2); zfree(tt);
		zmul(qq, v3, &tt); zsub(u3, tt, &t3); zfree(tt);
		zfree(qq); zfree(u1); zfree(u2);
		if ((u3.v != q->num.v) && (u3.v != q->den.v))
			zfree(u3);
		u1 = v1; u2 = v2; u3 = v3;
		v1 = t1; v2 = t2; v3 = t3;
	}
	if (u3.v != q->den.v)
		zfree(u3);
	zfree(v1);
	zfree(v2);
	i = ziszero(v3);
	zfree(v3);
	if (i && haveeps) {
		zfree(u1);
		zfree(u2);
		return qlink(q);
	}
	qtest.num = u2;
	qtest.den = u1;
	qtest.den.sign = 0;
	qtest.num.sign = q->num.sign;
	qtmp = qcopy(&qtest);
	zfree(u1);
	zfree(u2);
	return qtmp;
}


/*
 * Return an indication on whether or not two fractions are approximately
 * equal within the specified epsilon. Returns negative if the absolute value
 * of the difference between the two numbers is less than epsilon, zero if
 * the difference is equal to epsilon, and positive if the difference is
 * greater than epsilon.
 */
FLAG
qnear(q1, q2, epsilon)
	NUMBER *q1, *q2, *epsilon;
{
	int res;
	NUMBER qtemp, *qq;

	if (qisneg(epsilon))
		math_error("Negative epsilon for near");
	if (q1 == q2) {
		if (qiszero(epsilon))
			return 0;
		return -1;
	}
	if (qiszero(epsilon))
		return qcmp(q1, q2);
	if (qiszero(q2)) {
		qtemp = *q1;
		qtemp.num.sign = 0;
		return qrel(&qtemp, epsilon);
	}
	if (qiszero(q1)) {
		qtemp = *q2;
		qtemp.num.sign = 0;
		return qrel(&qtemp, epsilon);
	}
	qq = qsub(q1, q2);
	qtemp = *qq;
	qtemp.num.sign = 0;
	res = qrel(&qtemp, epsilon);
	qfree(qq);
	return res;
}


/*
 * Compute the gcd (greatest common divisor) of two numbers.
 *	q3 = qgcd(q1, q2);
 */
NUMBER *
qgcd(q1, q2)
	register NUMBER *q1, *q2;
{
	ZVALUE z;
	NUMBER *q;

	if (q1 == q2)
		return qabs(q1);
	if (qisfrac(q1) || qisfrac(q2)) {
		q = qalloc();
		zgcd(q1->num, q2->num, &q->num);
		zlcm(q1->den, q2->den, &q->den);
		return q;
	}
	if (qiszero(q1))
		return qabs(q2);
	if (qiszero(q2))
		return qabs(q1);
	if (qisunit(q1) || qisunit(q2))
		return qlink(&_qone_);
	zgcd(q1->num, q2->num, &z);
	if (zisunit(z)) {
		zfree(z);
		return qlink(&_qone_);
	}
	q = qalloc();
	q->num = z;
	return q;
}


/*
 * Compute the lcm (least common multiple) of two numbers.
 *	q3 = qlcm(q1, q2);
 */
NUMBER *
qlcm(q1, q2)
	register NUMBER *q1, *q2;
{
	NUMBER *q;

	if (qiszero(q1) || qiszero(q2))
		return qlink(&_qzero_);
	if (q1 == q2)
		return qabs(q1);
	if (qisunit(q1))
		return qabs(q2);
	if (qisunit(q2))
		return qabs(q1);
	q = qalloc();
	zlcm(q1->num, q2->num, &q->num);
	if (qisfrac(q1) || qisfrac(q2))
		zgcd(q1->den, q2->den, &q->den);
	return q;
}


/*
 * Remove all occurances of the specified factor from a number.
 * Returned number is always positive.
 */
NUMBER *
qfacrem(q1, q2)
	NUMBER *q1, *q2;
{
	long count;
	ZVALUE tmp;
	NUMBER *r;

	if (qisfrac(q1) || qisfrac(q2))
		math_error("Non-integers for factor removal");
	count = zfacrem(q1->num, q2->num, &tmp);
	if (zisunit(tmp)) {
		zfree(tmp);
		return qlink(&_qone_);
	}
	if (count == 0) {
		zfree(tmp);
		return qlink(q1);
	}
	r = qalloc();
	r->num = tmp;
	return r;
}


/*
 * Divide one number by the gcd of it with another number repeatedly until
 * the number is relatively prime.
 */
NUMBER *
qgcdrem(q1, q2)
	NUMBER *q1, *q2;
{
	ZVALUE tmp;
	NUMBER *r;

	if (qisfrac(q1) || qisfrac(q2))
		math_error("Non-integers for gcdrem");
	zgcdrem(q1->num, q2->num, &tmp);
	if (zisunit(tmp)) {
		zfree(tmp);
		return qlink(&_qone_);
	}
	if (zcmp(q1->num, tmp) == 0) {
		zfree(tmp);
		return qlink(q1);
	}
	r = qalloc();
	r->num = tmp;
	return r;
}


/*
 * Return the lowest prime factor of a number.
 * Search is conducted for the specified number of primes.
 * Returns one if no factor was found.
 */
NUMBER *
qlowfactor(q1, q2)
	NUMBER *q1, *q2;
{
	if (qisfrac(q1) || qisfrac(q2))
		math_error("Non-integers for lowfactor");
	return itoq(zlowfactor(q1->num, ztoi(q2->num)));
}


/*
 * Return the number of places after the decimal point needed to exactly
 * represent the specified number as a real number.  Integers return zero,
 * and non-terminating decimals return minus one.  Examples:
 *	qplaces(1/7)=-1, qplaces(3/10)= 1, qplaces(1/8)=3, qplaces(4)=0.
 */
long
qplaces(q)
	NUMBER *q;
{
	long twopow, fivepow;
	HALF fiveval[2];
	ZVALUE five;
	ZVALUE tmp;

	if (qisint(q))	/* no decimal places if number is integer */
		return 0;
	/*
	 * The number of decimal places of a fraction in lowest terms is finite
	 * if an only if the denominator is of the form 2^A * 5^B, and then the
	 * number of decimal places is equal to MAX(A, B).
	 */
	five.sign = 0;
	five.len = 1;
	five.v = fiveval;
	fiveval[0] = 5;
	fivepow = zfacrem(q->den, five, &tmp);
	if (!zisonebit(tmp)) {
		zfree(tmp);
		return -1;
	}
	twopow = zlowbit(tmp);
	zfree(tmp);
	if (twopow < fivepow)
		twopow = fivepow;
	return twopow;
}


/*
 * Perform a probabilistic primality test (algorithm P in Knuth).
 * Returns FALSE if definitely not prime, or TRUE if probably prime.
 * Second arg determines how many times to check for primality.
 */
BOOL
qprimetest(q1, q2)
	NUMBER *q1, *q2;
{
	if (qisfrac(q1) || qisfrac(q2) || qisneg(q2))
		math_error("Bad arguments for qprimetest");
	return zprimetest(q1->num, qtoi(q2));
}


/*
 * Return a trivial hash value for a number.
 */
HASH
qhash(q)
	NUMBER *q;
{
	HASH hash;

	hash = zhash(q->num);
	if (qisfrac(q))
		hash += zhash(q->den) * 2000003;
	return hash;
}

/* END CODE */