qfunc.c

Calc Software Package for Number Calc

/*
 * qfunc - extended precision rational arithmetic non-primitive functions
 *
 * Copyright (C) 1999-2004  David I. Bell and Ernest Bowen
 *
 * Primary author:  David I. Bell
 *
 * Calc is open software; you can redistribute it and/or modify it under
 * the terms of the version 2.1 of the GNU Lesser General Public License
 * as published by the Free Software Foundation.
 *
 * Calc is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
 * Public License for more details.
 *
 * A copy of version 2.1 of the GNU Lesser General Public License is
 * distributed with calc under the filename COPYING-LGPL.  You should have
 * received a copy with calc; if not, write to Free Software Foundation, Inc.
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 *
 * @(#) $Revision: 29.9 $
 * @(#) $Id: qfunc.c,v 29.9 2006/05/20 08:43:55 chongo Exp $
 * @(#) $Source: /usr/local/src/cmd/calc/RCS/qfunc.c,v $
 *
 * Under source code control:	1990/02/15 01:48:20
 * File existed as early as:	before 1990
 *
 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
 */


#include "qmath.h"
#include "config.h"
#include "prime.h"

static NUMBER **B_table;
static long B_num;
static long B_allocnum;
static NUMBER **E_table;
static long E_num;

#define QALLOCNUM 64

/*
 * Set the default epsilon for approximate calculations.
 * This must be greater than zero.
 *
 * given:
 *	q		number to be set as the new epsilon
 */
void
setepsilon(NUMBER *q)
{
	NUMBER *old;

	if (qisneg(q) || qiszero(q)) {
		math_error("Epsilon value must be greater than zero");
		/*NOTREACHED*/
	}
	old = conf->epsilon;
	conf->epsilonprec = qprecision(q);
	conf->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(NUMBER *q1, NUMBER *q2)
{
	NUMBER *r;
	ZVALUE z1, z2, tmp;
	int s, t;
	long rnd;

	if (qisfrac(q1) || qisfrac(q2)) {
		math_error("Non-integers for minv");
		/*NOTREACHED*/
	}
	if (qiszero(q2)) {
		if (qisunit(q1))
			return qlink(q1);
		return qlink(&_qzero_);
	}
	if (qisunit(q2))
		return qlink(&_qzero_);
	rnd = conf->mod;
	s = (rnd & 4) ? 0 : q2->num.sign;
	if (rnd & 1)
		s^= 1;

	tmp = q2->num;
	tmp.sign = 0;
	if (zmodinv(q1->num, tmp, &z1))
		return qlink(&_qzero_);
	zsub(tmp, z1, &z2);
	if (rnd & 16) {
		t = zrel(z1, z2);
		if (t < 0)
			s = 0;
		else if (t > 0)
			s = 1;
	}
	r = qalloc();
	if (s) {
		zfree(z1);
		z2.sign = TRUE;
		r->num = z2;
		return r;
	}
	zfree(z2);
	r->num = z1;
	return r;
}


/*
 * Return the residue modulo an integer (q3) of an integer (q1)
 * raised to a positive integer (q2) power.
 */
NUMBER *
qpowermod(NUMBER *q1, NUMBER *q2, NUMBER *q3)
{
	NUMBER *r;
	ZVALUE z1, z2, tmp;
	int s, t;
	long rnd;

	if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3)) {
		math_error("Non-integers for pmod");
		/*NOTREACHED*/
	}
	if (qisneg(q2)) {
		math_error("Negative power for pmod");
		/*NOTREACHED*/
	}
	if (qiszero(q3))
		return qpowi(q1, q2);
	if (qisunit(q3))
		return qlink(&_qzero_);
	rnd = conf->mod;
	s = (rnd & 4) ? 0 : q3->num.sign;
	if (rnd & 1)
		s^= 1;
	tmp = q3->num;
	tmp.sign = 0;
	zpowermod(q1->num, q2->num, tmp, &z1);
	if (ziszero(z1)) {
		zfree(z1);
		return qlink(&_qzero_);
	}
	zsub(tmp, z1, &z2);
	if (rnd & 16) {
		t = zrel(z1, z2);
		if (t < 0)
			s = 0;
		else if (t > 0)
			s = 1;
	}
	r = qalloc();
	if (s) {
		zfree(z1);
		z2.sign = TRUE;
		r->num = z2;
		return r;
	}
	zfree(z2);
	r->num = z1;
	return r;
}


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

	if (qisfrac(q2)) {
		math_error("Raising number to fractional power");
		/*NOTREACHED*/
	}
	num = q1->num;
	zden = 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");
			/*NOTREACHED*/
		}
		return qlink(&_qzero_);
	}
	if (zisunit(num) && zisunit(zden)) {	/* 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(zden))
		zpowi(zden, 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 hypotenuse within
 * the specified error.	 This is sqrt(a^2 + b^2).
 */
NUMBER *
qhypot(NUMBER *q1, NUMBER *q2, NUMBER *epsilon)
{
	NUMBER *tmp1, *tmp2, *tmp3;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for hypot");
		/*NOTREACHED*/
	}
	if (qiszero(q1))
		return qqabs(q2);
	if (qiszero(q2))
		return qqabs(q1);
	tmp1 = qsquare(q1);
	tmp2 = qsquare(q2);
	tmp3 = qqadd(tmp1, tmp2);
	qfree(tmp1);
	qfree(tmp2);
	tmp1 = qsqrt(tmp3, epsilon, 24L);
	qfree(tmp3);
	return tmp1;
}


/*
 * Given one leg of a right triangle with unit hypotenuse, 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(NUMBER *q, NUMBER *epsilon, BOOL wantneg)
{
	NUMBER *res, *qtmp1, *qtmp2;
	ZVALUE num;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for legtoleg");
		/*NOTREACHED*/
	}
	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");
		/*NOTREACHED*/
	}
	qtmp1 = qsquare(q);
	qtmp2 = qsub(&_qone_, qtmp1);
	qfree(qtmp1);
	res = qsqrt(qtmp2, epsilon, 24L);
	qfree(qtmp2);
	if (wantneg) {
		qtmp1 = qneg(res);
		qfree(res);
		res = qtmp1;
	}
	return res;
}


/*
 * Compute the square root of a real number.
 * Type of rounding if any depends on rnd.
 * If rnd & 32 is nonzero, result is exact for square numbers;
 * If rnd & 64 is nonzero, the negative square root is returned;
 * If rnd < 32, result is rounded to a multiple of epsilon
 *	up, down, etc. depending on bits 0, 2, 4 of v.
 */

NUMBER *
qsqrt(NUMBER *q1, NUMBER *epsilon, long rnd)
{
	NUMBER *r, etemp;
	ZVALUE tmp1, tmp2, quo, mul, divisor;
	long  s1, s2, up, RR, RS;
	int sign;

	if (qisneg(q1)) {
		math_error("Square root of negative number");
		/*NOTREACHED*/
	}
	if (qiszero(q1))
		return qlink(&_qzero_);
	sign = (rnd & 64) != 0;
	if (qiszero(epsilon)) {
		math_error("Zero epsilon for qsqrt");
		/*NOTREACHED*/
	}

	etemp = *epsilon;
	etemp.num.sign = 0;
	RS = rnd & 25;
	if (!(RS & 8))
		RS ^= epsilon->num.sign;
	if (rnd & 2)
		RS ^= sign ^ epsilon->num.sign;
	if (rnd & 4)
		RS ^= epsilon->num.sign;
	RR = zisunit(q1->den) && qisunit(epsilon);
	if (rnd & 32 || RR) {
		s1 = zsqrt(q1->num, &tmp1, RS);
		if (RR) {
			if (ziszero(tmp1)) {
				zfree(tmp1);
				return qlink(&_qzero_);
			}
			r = qalloc();
			tmp1.sign = sign;
			r->num = tmp1;
			return r;
		}
		if (!s1) {
			s2 = zsqrt(q1->den, &tmp2, 0);
			if (!s2) {
				r = qalloc();
				tmp1.sign = sign;
				r->num = tmp1;
				r->den = tmp2;
				return r;
			}
			zfree(tmp2);
		}
		zfree(tmp1);
	}
	up = 0;
	zsquare(epsilon->den, &tmp1);
	zmul(tmp1, q1->num, &tmp2);
	zfree(tmp1);
	zsquare(epsilon->num, &tmp1);
	zmul(tmp1, q1->den, &divisor);
	zfree(tmp1);
	if (rnd & 16) {
		zshift(tmp2, 2, &tmp1);
		zfree(tmp2);
		s1 = zquo(tmp1, divisor, &quo, 16);
		zfree(tmp1);
		s2 = zsqrt(quo, &tmp1, s1 ? s1 < 0 : 16);
		zshift(tmp1, -1, &mul);
		up = (*tmp1.v & 1) ? s1 + s2 : -1;
		zfree(tmp1);
	} else {
		s1 = zquo(tmp2, divisor, &quo, 0);
		zfree(tmp2);
		s2 = zsqrt(quo, &mul, 0);
		up = (s1 + s2) ? 0 : -1;
	}
	if (up == 0) {
		if (rnd & 8)
			up = (long)((RS ^ *mul.v) & 1);
		else
			up = RS ^ sign;
	}
	if (up > 0) {
		zadd(mul, _one_, &tmp2);
		zfree(mul);
		mul = tmp2;
	}
	zfree(divisor);
	zfree(quo);
	if (ziszero(mul)) {
		zfree(mul);
		return qlink(&_qzero_);
	}
	r = qalloc();
	zreduce(mul, etemp.den, &tmp1, &r->den);
	zfree(mul);
	tmp1.sign = sign;
	zmul(tmp1, etemp.num, &r->num);
	zfree(tmp1);
	return r;
}


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

	if (qisneg(q)) {
		math_error("Square root of negative number");
		/*NOTREACHED*/
	}
	if (qiszero(q))
		return qlink(&_qzero_);
	r = qalloc();
	if (qisint(q)) {
		(void) zsqrt(q->num, &r->num,0);
		return r;
	}
	zquo(q->num, q->den, &tmp, 0);
	(void) zsqrt(tmp, &r->num,0);
	freeh(tmp.v);
	return r;
}

/*
 * Return whether or not a number is an exact square.
 */
BOOL
qissquare(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(NUMBER *q1, NUMBER *q2)
{
	NUMBER *r;
	ZVALUE tmp;

	if (qisneg(q2) || qiszero(q2) || qisfrac(q2)) {
		math_error("Taking number to bad root value");
		/*NOTREACHED*/
	}
	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, 0);
	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.
 *
 * given:
 *	q		number to take log of
 */
long
qilog2(NUMBER *q)
{
	long n;			/* power of two */
	int c;			/* result of comparison */
	ZVALUE tmp1, tmp2;	/* temporary values */

	if (qiszero(q)) {
		math_error("Zero argument for ilog2");
		/*NOTREACHED*/
	}
	if (qisint(q))
		return zhighbit(q->num);
	tmp1 = q->num;
	tmp1.sign = 0;
	n = zhighbit(tmp1) - zhighbit(q->den);
	if (n == 0)
		c = zrel(tmp1, q->den);
	else if (n > 0) {
		zshift(q->den, n, &tmp2);
		c = zrel(tmp1, tmp2);
		zfree(tmp2);
	} else {
		zshift(tmp1, -n, &tmp2);
		c = zrel(tmp2, q->den);
		zfree(tmp2);
	}
	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.
 *
 * given:
 *	q		number to take log of
 */
long
qilog10(NUMBER *q)
{
	long n;			/* log value */
	ZVALUE tmp1, tmp2;	/* temporary values */

	if (qiszero(q)) {
		math_error("Zero argument for ilog10");
		/*NOTREACHED*/
	}
	tmp1 = q->num;
	tmp1.sign = 0;
	if (qisint(q))
		return zlog10(tmp1, NULL);
	/*
	 * The number is not an integer.
	 * Compute the result if the number is greater than one.
	 */
	if (zrel(tmp1, q->den) > 0) {
		zquo(tmp1, q->den, &tmp2, 0);
		n = zlog10(tmp2, NULL);
		zfree(tmp2);
		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(tmp1))
		zsub(q->den, _one_, &tmp2);
	else
		zquo(q->den, tmp1, &tmp2, 0);
	n = -zlog10(tmp2, NULL) - 1;
	zfree(tmp2);
	return n;
}

/*
 * Return the integer floor of the logarithm of a number relative to
 * a specified integral base.
 */
NUMBER *
qilog(NUMBER *q, ZVALUE base)
{
	long n;
	ZVALUE tmp1, tmp2;

	if (qiszero(q))
		return NULL;

	if (qisunit(q))
		return qlink(&_qzero_);
	if (qisint(q))
		return itoq(zlog(q->num, base));
	tmp1 = q->num;
	tmp1.sign = 0;
	if (zrel(tmp1, q->den) > 0) {
		zquo(tmp1, q->den, &tmp2, 0);
		n = zlog(tmp2, base);
		zfree(tmp2);
		return itoq(n);
	}
	if (zisunit(tmp1))
		zsub(q->den, _one_, &tmp2);
	else
		zquo(q->den, tmp1, &tmp2, 0);
	n = -zlog(tmp2, base) - 1;