qfunc.c

Calc Software Package for Number Calc

qbtrunc(NUMBER *q1, NUMBER *q2)
{
	long places;
	NUMBER *r, *e;

	if (qisfrac(q2) || !zistiny(q2->num)) {
		math_error("Bad number of places for qtrunc");
		/*NOTREACHED*/
	}
	places = qtoi(q2);
	e = qbitvalue(-places);
	r = qmappr(q1, e, 2);
	qfree(e);
	return r;
}


/*
 * Round a number to a specified number of binary places.
 */
NUMBER *
qbround(NUMBER *q, long places, long rnd)
{
	NUMBER *e, *r;

	if (qiszero(q))
		return qlink(&_qzero_);
	if (rnd & 32)
		places -= qilog2(q) + 1;
	e = qbitvalue(-places);
	r = qmappr(q, e, rnd & 31);
	qfree(e);
	return r;
}

/*
 * Round a number to a specified number of decimal digits.
 */
NUMBER *
qround(NUMBER *q, long places, long rnd)
{
	NUMBER *e, *r;

	if (qiszero(q))
		return qlink(&_qzero_);
	if (rnd & 32)
		places -= qilog10(q) + 1;
	e = qtenpow(-places);
	r = qmappr(q, e, rnd & 31);
	qfree(e);
	return r;
}

/*
 * Approximate a number to nearest multiple of a given number.	Whether
 * rounding is down, up, etc. is determined by rnd.
 */
NUMBER *
qmappr(NUMBER *q, NUMBER *e, long rnd)
{
	NUMBER *r;
	ZVALUE tmp1, tmp2, mul;

	if (qiszero(e))
		return qlink(q);
	if (qiszero(q))
		return qlink(&_qzero_);
	zmul(q->num, e->den, &tmp1);
	zmul(q->den, e->num, &tmp2);
	zquo(tmp1, tmp2, &mul, rnd);
	zfree(tmp1);
	zfree(tmp2);
	if (ziszero(mul)) {
		zfree(mul);
		return qlink(&_qzero_);
	}
	r = qalloc();
	zreduce(mul, e->den, &tmp1, &r->den);
	zmul(tmp1, e->num, &r->num);
	zfree(tmp1);
	zfree(mul);
	return r;
}


/*
 * Determine the smallest-denominator rational number in the interval of
 * length abs(epsilon) (< 1) with centre or one end at q, or to determine
 * the number nearest above or nearest below q with denominator
 * not exceeding abs(epsilon).
 * Whether the approximation is nearest above or nearest below is
 * determined by rnd and the signs of epsilon and q.
 */

NUMBER *
qcfappr(NUMBER *q, NUMBER *epsilon, long rnd)
{
	NUMBER *res, etemp, *epsilon1;
	ZVALUE num, zden, oldnum, oldden;
	ZVALUE rem, oldrem, quot;
	ZVALUE tmp1, tmp2, tmp3, tmp4;
	ZVALUE denbnd;
	ZVALUE f, g, k;
	BOOL esign;
	int s;
	BOOL bnddencase;
	BOOL useold = FALSE;

	if (qiszero(epsilon) || qisint(q))
		return qlink(q);

	esign = epsilon->num.sign;
	etemp = *epsilon;
	etemp.num.sign = 0;
	bnddencase = (zrel(etemp.num, etemp.den) >= 0);
	if (bnddencase) {
		zquo(etemp.num, etemp.den, &denbnd, 0);
		if (zrel(q->den, denbnd) <= 0) {
			zfree(denbnd);
			return qlink(q);
		}
	} else {
		if (rnd & 16)
			epsilon1 = qscale(epsilon, -1);
		else
			epsilon1 = qlink(epsilon);
		zreduce(q->den, epsilon1->den, &tmp1, &g);
		zmul(epsilon1->num, tmp1, &f);
		f.sign = 0;
		zfree(tmp1);
		qfree(epsilon1);
	}
	if (rnd & 16 && !zistwo(q->den)) {
		s = 0;
	} else {
		s = esign ? -1 : 1;
		if (rnd & 1)
			s = -s;
		if (rnd & 2 && q->num.sign ^ esign)
			s = -s;
		if (rnd & 4 && esign)
			s = -s;
	}
	oldnum = _one_;
	oldden = _zero_;
	zcopy(q->den, &oldrem);
	zdiv(q->num, q->den, &num, &rem, 0);
	zden = _one_;
	for (;;) {
		if (!bnddencase) {
			zmul(f, zden, &tmp1);
			zmul(g, rem, &tmp2);
			if (ziszero(rem) || (s >= 0 && zrel(tmp1,tmp2) >= 0))
				break;
			zfree(tmp1);
			zfree(tmp2);
		}
		zdiv(oldrem, rem, &quot, &tmp1, 0);
		zfree(oldrem);
		oldrem = rem;
		rem = tmp1;
		zmul(quot, zden, &tmp1);
		zadd(tmp1, oldden, &tmp2);
		zfree(tmp1);
		zfree(oldden);
		oldden = zden;
		zden = tmp2;
		zmul(quot, num, &tmp1);
		zadd(tmp1, oldnum, &tmp2);
		zfree(tmp1);
		zfree(oldnum);
		oldnum = num;
		num = tmp2;
		zfree(quot);
		if (bnddencase) {
			if (zrel(zden, denbnd) >= 0)
				break;
		}
		s = -s;
	}
	if (bnddencase) {
		if (s > 0) {
			useold = TRUE;
		} else {
			zsub(zden, denbnd, &tmp1);
			zquo(tmp1, oldden, &k, 1);
			zfree(tmp1);
		}
		zfree(denbnd);
	} else {
		if (s < 0) {
			zfree(tmp1);
			zfree(tmp2);
			zfree(f);
			zfree(g);
			zfree(oldnum);
			zfree(oldden);
			zfree(num);
			zfree(zden);
			zfree(oldrem);
			zfree(rem);
			return qlink(q);
		}
		zsub(tmp1, tmp2, &tmp3);
		zfree(tmp1);
		zfree(tmp2);
		zmul(f, oldden, &tmp1);
		zmul(g, oldrem, &tmp2);
		zfree(f);
		zfree(g);
		zadd(tmp1, tmp2, &tmp4);
		zfree(tmp1);
		zfree(tmp2);
		zquo(tmp3, tmp4, &k, 0);
		zfree(tmp3);
		zfree(tmp4);
	}
	if (!useold && !ziszero(k)) {
		zmul(k, oldnum, &tmp1);
		zsub(num, tmp1, &tmp2);
		zfree(tmp1);
		zfree(num);
		num = tmp2;
		zmul(k, oldden, &tmp1);
		zsub(zden, tmp1, &tmp2);
		zfree(tmp1);
		zfree(zden);
		zden = tmp2;
	}
	if (bnddencase && s == 0) {
		zmul(k, oldrem, &tmp1);
		zadd(rem, tmp1, &tmp2);
		zfree(tmp1);
		zfree(rem);
		rem = tmp2;
		zmul(rem, oldden, &tmp1);
		zmul(zden, oldrem, &tmp2);
		useold = (zrel(tmp1, tmp2) >= 0);
		zfree(tmp1);
		zfree(tmp2);
	}
	if (!bnddencase || s <= 0)
		zfree(k);
	zfree(rem);
	zfree(oldrem);
	res = qalloc();
	if (useold) {
		zfree(num);
		zfree(zden);
		res->num = oldnum;
		res->den = oldden;
		return res;
	}
	zfree(oldnum);
	zfree(oldden);
	res->num = num;
	res->den = zden;
	return res;
}


/*
 * Calculate the nearest-above, or nearest-below, or nearest, number
 * with denominator less than the given number, the choice between
 * possibilities being determined by the parameter rnd.
 */
NUMBER *
qcfsim(NUMBER *q, long rnd)
{
	NUMBER *res;
	ZVALUE tmp1, tmp2, den1, den2;
	int s;

	if (qiszero(q) && rnd & 26)
		return qlink(&_qzero_);
	if (rnd & 24) {
		s = q->num.sign;
	} else {
		s = rnd & 1;
		if (rnd & 2)
			s ^= q->num.sign;
	}
	if (qisint(q)) {
		if ((rnd & 8) && !(rnd & 16))
			return qlink(&_qzero_);
		if (s)
			return qinc(q);
		return qdec(q);
	}
	if (zistwo(q->den)) {
		if (rnd & 16)
			s ^= 1;
		if (s)
			zadd(q->num, _one_, &tmp1);
		else
			zsub(q->num, _one_, &tmp1);
		res = qalloc();
		zshift(tmp1, -1, &res->num);
		zfree(tmp1);
		return res;
	}
	s = s ? 1 : -1;
	if (rnd & 24)
		s = 0;
	res = qalloc();
	zmodinv(q->num, q->den, &den1);
	if (s >= 0) {
		zsub(q->den, den1, &den2);
		if (s > 0 || ((zrel(den1, den2) < 0) ^ (!(rnd & 16)))) {
			zfree(den1);
			res->den = den2;
			zmul(den2, q->num, &tmp1);
			zadd(tmp1, _one_, &tmp2);
			zfree(tmp1);
			zequo(tmp2, q->den, &res->num);
			zfree(tmp2);
			return res;
		}
		zfree(den2);
	}
	res->den = den1;
	zmul(den1, q->num, &tmp1);
	zsub(tmp1, _one_, &tmp2);
	zfree(tmp1);
	zequo(tmp2, q->den, &res->num);
	zfree(tmp2);
	return res;
}



/*
 * 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(NUMBER *q1, NUMBER *q2, NUMBER *epsilon)
{
	int res;
	NUMBER qtemp, etemp, *qq;

	etemp = *epsilon;
	etemp.num.sign = 0;
	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, &etemp);
	}
	if (qiszero(q1)) {
		qtemp = *q2;
		qtemp.num.sign = 0;
		return qrel(&qtemp, &etemp);
	}
	qq = qsub(q1, q2);
	qtemp = *qq;
	qtemp.num.sign = 0;
	res = qrel(&qtemp, &etemp);
	qfree(qq);
	return res;
}


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

	if (q1 == q2)
		return qqabs(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 qqabs(q2);
	if (qiszero(q2))
		return qqabs(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(NUMBER *q1, NUMBER *q2)
{
	NUMBER *q;

	if (qiszero(q1) || qiszero(q2))
		return qlink(&_qzero_);
	if (q1 == q2)
		return qqabs(q1);
	if (qisunit(q1))
		return qqabs(q2);
	if (qisunit(q2))
		return qqabs(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 occurrences of the specified factor from a number.
 * Returned number is always positive or zero.
 */
NUMBER *
qfacrem(NUMBER *q1, NUMBER *q2)
{
	long count;
	ZVALUE tmp;
	NUMBER *r;

	if (qisfrac(q1) || qisfrac(q2)) {
		math_error("Non-integers for factor removal");
		/*NOTREACHED*/
	}
	if (qiszero(q2))
		return qqabs(q1);
	if (qiszero(q1))
		return qlink(&_qzero_);
	count = zfacrem(q1->num, q2->num, &tmp);
	if (zisunit(tmp)) {
		zfree(tmp);
		return qlink(&_qone_);
	}
	if (count == 0 && !qisneg(q1)) {
		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(NUMBER *q1, NUMBER *q2)
{
	ZVALUE tmp;
	NUMBER *r;

	if (qisfrac(q1) || qisfrac(q2)) {
		math_error("Non-integers for gcdrem");
		/*NOTREACHED*/
	}
	if (qiszero(q2))
		return qlink(&_qone_);
	if (qiszero(q1))
		return qlink(&_qzero_);
	if (zgcdrem(q1->num, q2->num, &tmp) == 0)
		return qqabs(q1);
	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(NUMBER *q1, NUMBER *q2)
{
	unsigned long count;

	if (qisfrac(q1) || qisfrac(q2)) {
		math_error("Non-integers for lowfactor");
		/*NOTREACHED*/
	}
	count = ztoi(q2->num);
	if (count > PIX_32B) {
		math_error("lowfactor count is too large");
		/*NOTREACHED*/
	}
	return utoq(zlowfactor(q1->num, count));
}

/*
 * 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:
 * qdecplaces(1.23) = 2, qdecplaces(3) = 0, qdecplaces(1/7) = -1.
 */
long
qdecplaces(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;
}


/*
 * Return, if possible, the minimum number of places after the decimal
 * point needed to exactly represent q with the specified base.
 * Integers return 0 and numbers with non-terminating expansions -1.
 * Returns -2 if base inadmissible
 */
long
qplaces(NUMBER *q, ZVALUE base)
{
	long count;
	ZVALUE tmp;

	if (base.len == 1 && base.v[0] == 10)
		return qdecplaces(q);
	if (ziszero(base) || zisunit(base))
			return -2;
	if (qisint(q))
		return 0;
	if (zisonebit(base)) {
		if (!zisonebit(q->den))
			return -1;
		return 1 + (zlowbit(q->den) - 1)/zlowbit(base);
	}
	count = zgcdrem(q->den, base, &tmp);
	if (count == 0)
		return -1;
	if (!zisunit(tmp))
		count = -1;
	zfree(tmp);
	return count;
}


/*
 * Perform a probabilistic primality test (algorithm P in Knuth).
 * Returns FALSE if definitely not prime, or TRUE if probably prime.
 *
 * The absolute value of the 2nd arg determines how many times
 * to check for primality.  If 2nd arg < 0, then the trivial
 * check is omitted.  The 3rd arg determines how tests to
 * initially skip.
 */
BOOL
qprimetest(NUMBER *q1, NUMBER *q2, NUMBER *q3)
{
	if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3)) {
		math_error("Bad arguments for ptest");
		/*NOTREACHED*/
	}
	if (zge24b(q2->num)) {
		math_error("ptest count >= 2^24");
		/*NOTREACHED*/
	}
	return zprimetest(q1->num, ztoi(q2->num), q3->num);
}