qfunc.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 non-primitive functions
 */

#include "qmath.h"


NUMBER *_epsilon_;	/* default precision for real functions */
long _epsilonprec_;	/* bits of precision for epsilon */


/*
 * Set the default precision for real calculations.
 * The precision must be between zero and one.
 */
void
setepsilon(q)
	NUMBER *q;		/* number to be set as the new epsilon */
{
	NUMBER *old;

	if (qisneg(q) || qiszero(q) || (qreli(q, 1L) >= 0))
		math_error("Epsilon value must be between zero and one");
	old = _epsilon_;
	_epsilonprec_ = qprecision(q);
	_epsilon_ = qlink(q);
	if (old)
		qfree(old);
}


/*
 * Return the inverse of one number modulo another.
 * That is, find x such that:
 *	Ax = 1 (mod B)
 * Returns zero if the numbers are not relatively prime (temporary hack).
 */
NUMBER *
qminv(q1, q2)
	NUMBER *q1, *q2;
{
	NUMBER *r;

	if (qisfrac(q1) || qisfrac(q2))
		math_error("Non-integers for minv");
	r = qalloc();
	if (zmodinv(q1->num, q2->num, &r->num)) {
		qfree(r);
		return qlink(&_qzero_);
	}
	return r;
}


/*
 * Return the residue modulo an integer (q3) of an integer (q1) 
 * raised to a positive integer (q2) power.  
 */
NUMBER *
qpowermod(q1, q2, q3)
	NUMBER *q1, *q2, *q3;
{
	NUMBER *r, *t;
	BOOL s;

	if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3)) 
		math_error("Non-integers for pmod");
	if (qisneg(q2))
		math_error("Negative power for pmod");
	if (qiszero(q3)) return qpowi(q1, q2);
	s = q3->num.sign;
	q3->num.sign = 0;
	r = qalloc();
	zpowermod(q1->num, q2->num, q3->num, &r->num);
	if (!s || qiszero(r)) return r;
	q3->num.sign = 1;
	t = qadd(r, q3);
	qfree(r);
	return t;
}


/*
 * Return the power function of one number with another.
 * The power must be integral.
 *	q3 = qpowi(q1, q2);
 */
NUMBER *
qpowi(q1, q2)
	NUMBER *q1, *q2;
{
	register NUMBER *r;
	BOOL invert, sign;
	ZVALUE num, den, z2;

	if (qisfrac(q2))
		math_error("Raising number to fractional power");
	num = q1->num;
	den = q1->den;
	z2 = q2->num;
	sign = (num.sign && zisodd(z2));
	invert = z2.sign;
	num.sign = 0;
	z2.sign = 0;
	/*
	* Check for trivial cases first.
	*/
	if (ziszero(num) && !ziszero(z2)) {	/* zero raised to a power */
		if (invert)
			math_error("Zero raised to negative power");
		return qlink(&_qzero_);
	}
	if (zisunit(num) && zisunit(den)) {	/* 1 or -1 raised to a power */
		r = (sign ? q1 : &_qone_);
		r->links++;
		return r;
	}
	if (ziszero(z2))	/* raising to zeroth power */
		return qlink(&_qone_);
	if (zisunit(z2)) {	/* raising to power 1 or -1 */
		if (invert)
			return qinv(q1);
		return qlink(q1);
	}
	/*
	 * Not a trivial case.  Do the real work.
	 */
	r = qalloc();
	if (!zisunit(num))
		zpowi(num, z2, &r->num);
	if (!zisunit(den))
		zpowi(den, z2, &r->den);
	if (invert) {
		z2 = r->num;
		r->num = r->den;
		r->den = z2;
	}
	r->num.sign = sign;
	return r;
}


/*
 * Given the legs of a right triangle, compute its hypothenuse within
 * the specified error.  This is sqrt(a^2 + b^2).
 */
NUMBER *
qhypot(q1, q2, epsilon)
	NUMBER *q1, *q2, *epsilon;
{
	NUMBER *tmp1, *tmp2, *tmp3;

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Bad epsilon value for hypot");
	if (qiszero(q1))
		return qabs(q2);
	if (qiszero(q2))
		return qabs(q1);
	tmp1 = qsquare(q1);
	tmp2 = qsquare(q2);
	tmp3 = qadd(tmp1, tmp2);
	qfree(tmp1);
	qfree(tmp2);
	tmp1 = qsqrt(tmp3, epsilon);
	qfree(tmp3);
	return tmp1;
}


/*
 * Given one leg of a right triangle with unit hypothenuse, calculate
 * the other leg within the specified error.  This is sqrt(1 - a^2).
 * If the wantneg flag is nonzero, then negative square root is returned.
 */
NUMBER *
qlegtoleg(q, epsilon, wantneg)
	NUMBER *q, *epsilon;
	BOOL wantneg;
{
	NUMBER *res, *qtmp1, *qtmp2;
	ZVALUE num;

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Bad epsilon value for legtoleg");
	if (qisunit(q))
		return qlink(&_qzero_);
	if (qiszero(q)) {
		if (wantneg)
			return qlink(&_qnegone_);
		return qlink(&_qone_);
	}
	num = q->num;
	num.sign = 0;
	if (zrel(num, q->den) >= 0)
		math_error("Leg too large in legtoleg");
	qtmp1 = qsquare(q);
	qtmp2 = qsub(&_qone_, qtmp1);
	qfree(qtmp1);
	res = qsqrt(qtmp2, epsilon);
	qfree(qtmp2);
	if (wantneg) {
		qtmp1 = qneg(res);
		qfree(res);
		res = qtmp1;
	}
	return res;
}

/*
 * Compute the square root of a number to within the specified error.
 * If the number is an exact square, the exact result is returned.
 *	q3 = qsqrt(q1, q2);
 */
NUMBER *
qsqrt(q1, epsilon)
	register NUMBER *q1, *epsilon;
{
	long bits, bits2;	/* number of bits of precision */
	int exact;
	NUMBER *r;
	ZVALUE t1, t2;

	if (qisneg(q1))
		math_error("Square root of negative number");
	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Bad epsilon value for sqrt");
	if (qiszero(q1))
		return qlink(&_qzero_);
	if (qisunit(q1))
		return qlink(&_qone_);
	if (qiszero(epsilon) && qisint(q1) && zistiny(q1->num) && (*q1->num.v <= 3))
		return qlink(&_qone_);
	bits = zhighbit(epsilon->den) - zhighbit(epsilon->num) + 1;
	if (bits < 0)
		bits = 0;
	bits2 = zhighbit(q1->num) - zhighbit(q1->den) + 1;
	if (bits2 > 0)
		bits += bits2;
	r = qalloc();
	zshift(q1->num, bits * 2, &t2);
	zmul(q1->den, t2, &t1);
	zfree(t2);
	exact = zsqrt(t1, &t2);
	zfree(t1);
	if (exact) {
		zshift(q1->den, bits, &t1);
		zreduce(t2, t1, &r->num, &r->den);
	} else {
		zquo(t2, q1->den, &t1);
		zfree(t2);
		zbitvalue(bits, &t2);
		zreduce(t1, t2, &r->num, &r->den);
	}
	zfree(t1);
	zfree(t2);
	if (qiszero(r)) {
		qfree(r);
		r = qlink(&_qzero_);
	}
	return r;
}


/*
 * Calculate the integral part of the square root of a number.
 * Example:  qisqrt(13) = 3.
 */
NUMBER *
qisqrt(q)
	NUMBER *q;
{
	NUMBER *r;
	ZVALUE tmp;

	if (qisneg(q))
		math_error("Square root of negative number");
	if (qiszero(q))
		return qlink(&_qzero_);
	if (qisint(q) && zistiny(q->num) && (z1tol(q->num) < 4))
		return qlink(&_qone_);
	r = qalloc();
	if (qisint(q)) {
		(void) zsqrt(q->num, &r->num);
		return r;
	}
	zquo(q->num, q->den, &tmp);
	(void) zsqrt(tmp, &r->num);
	zfree(tmp);
	return r;
}


/*
 * Return whether or not a number is an exact square.
 */
BOOL
qissquare(q)
	NUMBER *q;
{
	BOOL flag;

	flag = zissquare(q->num);
	if (qisint(q) || !flag)
		return flag;
	return zissquare(q->den);
}


/*
 * Compute the greatest integer of the Kth root of a number.
 * Example:  qiroot(85, 3) = 4.
 */
NUMBER *
qiroot(q1, q2)
	register NUMBER *q1, *q2;
{
	NUMBER *r;
	ZVALUE tmp;

	if (qisneg(q2) || qiszero(q2) || qisfrac(q2))
		math_error("Taking number to bad root value");
	if (qiszero(q1))
		return qlink(&_qzero_);
	if (qisone(q1) || qisone(q2))
		return qlink(q1);
	if (qistwo(q2))
		return qisqrt(q1);
	r = qalloc();
	if (qisint(q1)) {
		zroot(q1->num, q2->num, &r->num);
		return r;
	}
	zquo(q1->num, q1->den, &tmp);
	zroot(tmp, q2->num, &r->num);
	zfree(tmp);
	return r;
}


/*
 * Return the greatest integer of the base 2 log of a number.
 * This is the number such that  1 <= x / log2(x) < 2.
 * Examples:  qilog2(8) = 3, qilog2(1.3) = 1, qilog2(1/7) = -3.
 */
long
qilog2(q)
	NUMBER *q;		/* number to take log of */
{
	long n;			/* power of two */
	int c;			/* result of comparison */
	ZVALUE tmp;		/* temporary value */

	if (qisneg(q) || qiszero(q))
		math_error("Non-positive number for log2");
	if (qisint(q))
		return zhighbit(q->num);
	n = zhighbit(q->num) - zhighbit(q->den);
	if (n == 0)
		c = zrel(q->num, q->den);
	else if (n > 0) {
		zshift(q->den, n, &tmp);
		c = zrel(q->num, tmp);
	} else {
		zshift(q->num, n, &tmp);
		c = zrel(tmp, q->den);
	}
	if (n)
		zfree(tmp);
	if (c < 0)
		n--;
	return n;
}


/*
 * Return the greatest integer of the base 10 log of a number.
 * This is the number such that  1 <= x / log10(x) < 10.
 * Examples:  qilog10(100) = 2, qilog10(12.3) = 1, qilog10(.023) = -2.
 */
long
qilog10(q)
	NUMBER *q;		/* number to take log of */
{
	long n;			/* log value */
	ZVALUE temp;		/* temporary value */

	if (qisneg(q) || qiszero(q))
		math_error("Non-positive number for log10");
	if (qisint(q))
		return zlog10(q->num);
	/*
	 * The number is not an integer.
	 * Compute the result if the number is greater than one.
	 */
	if ((q->num.len > q->den.len) ||
		((q->num.len == q->den.len) && (zrel(q->num, q->den) > 0))) {
			zquo(q->num, q->den, &temp);
			n = zlog10(temp);
			zfree(temp);
			return n;
	}
	/*
	 * Here if the number is less than one.
	 * If the number is the inverse of a power of ten, then the obvious answer
	 * will be off by one.  Subtracting one if the number is the inverse of an
	 * integer will fix it.
	 */
	if (zisunit(q->num))
		zsub(q->den, _one_, &temp);
	else
		zquo(q->den, q->num, &temp);
	n = -zlog10(temp) - 1;
	zfree(temp);
	return n;
}


/*
 * Return the number of digits in a number, ignoring the sign.
 * For fractions, this is the number of digits of its greatest integer.
 * Examples: qdigits(3456) = 4, qdigits(-23.45) = 2, qdigits(.0120) = 1.
 */
long
qdigits(q)
	NUMBER *q;		/* number to count digits of */
{
	long n;			/* number of digits */
	ZVALUE temp;		/* temporary value */

	if (qisint(q))
		return zdigits(q->num);
	zquo(q->num, q->den, &temp);
	n = zdigits(temp);
	zfree(temp);
	return n;
}


/*
 * Return the digit at the specified decimal place of a number represented
 * in floating point.  The lowest digit of the integral part of a number
 * is the zeroth decimal place.  Negative decimal places indicate digits
 * to the right of the decimal point.  Examples: qdigit(1234.5678, 1) = 3,
 * qdigit(1234.5678, -3) = 7.
 */
FLAG
qdigit(q, n)
	NUMBER *q;
	long n;
{
	ZVALUE tenpow, tmp1, tmp2;
	FLAG res;

	/*
	 * Zero number or negative decimal place of integer is trivial.
	 */
	if (qiszero(q) || (qisint(q) && (n < 0)))
		return 0;
	/*
	 * For non-negative decimal places, answer is easy.
	 */
	if (n >= 0) {
		if (qisint(q))
			return zdigit(q->num, n);
		zquo(q->num, q->den, &tmp1);
		res = zdigit(tmp1, n);
		zfree(tmp1);
		return res;
	}
	/*
	 * Fractional value and want negative digit, must work harder.
	 */
	ztenpow(-n, &tenpow);
	zmul(q->num, tenpow, &tmp1);
	zfree(tenpow);
	zquo(tmp1, q->den, &tmp2);
	res = zmodi(tmp2, 10L);
	zfree(tmp1);
	zfree(tmp2);
	return res;
}


/*
 * Return whether or not a bit is set at the specified bit position in a
 * number.  The lowest bit of the integral part of a number is the zeroth
 * bit position.  Negative bit positions indicate bits to the right of the
 * binary decimal point.  Examples: qdigit(17.1, 0) = 1, qdigit(17.1, -1) = 0.
 */
BOOL
qisset(q, n)
	NUMBER *q;
	long n;
{
	NUMBER *qtmp1, *qtmp2;
	ZVALUE ztmp;
	BOOL res;

	/*
	 * Zero number or negative bit position place of integer is trivial.
	 */
	if (qiszero(q) || (qisint(q) && (n < 0)))
		return FALSE;
	/*
	 * For non-negative bit positions, answer is easy.
	 */
	if (n >= 0) {
		if (qisint(q))
			return zisset(q->num, n);
		zquo(q->num, q->den, &ztmp);
		res = zisset(ztmp, n);
		zfree(ztmp);
		return res;
	}
	/*
	 * Fractional value and want negative bit position, must work harder.
	 */
	qtmp1 = qscale(q, -n);
	qtmp2 = qint(qtmp1);
	qfree(qtmp1);
	res = ((qtmp2->num.v[0] & 0x01) != 0);
	qfree(qtmp2);
	return res;
}


/*
 * Compute the factorial of an integer.
 *	q2 = qfact(q1);
 */
NUMBER *
qfact(q)
	register NUMBER *q;
{
	register NUMBER *r;

	if (qisfrac(q))
		math_error("Non-integral factorial");
	if (qiszero(q) || zisone(q->num))
		return qlink(&_qone_);
	r = qalloc();
	zfact(q->num, &r->num);
	return r;
}


/*
 * Compute the product of the primes less than or equal to a number.
 *	q2 = qpfact(q1);
 */
NUMBER *
qpfact(q)
	register NUMBER *q;
{
	NUMBER *r;

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


/*
 * Compute the lcm of all the numbers less than or equal to a number.
 *	q2 = qlcmfact(q1);
 */
NUMBER *
qlcmfact(q)
	register NUMBER *q;
{
	NUMBER *r;

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


/*
 * Compute the permutation function  M! / (M - N)!.
 */
NUMBER *
qperm(q1, q2)
	register NUMBER *q1, *q2;
{
	register NUMBER *r;

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


/*
 * Compute the combinatorial function  M! / (N! * (M - N)!).
 */
NUMBER *
qcomb(q1, q2)