qfunc.c

Calc Software Package for Number Calc

	zfree(tmp2);
	return itoq(n);
}


/*
 * Return the number of digits in the representation to a specified
 * base of the integral part of a number.
 *
 * Examples: qdigits(3456,10) = 4, qdigits(-23.45, 10) = 2.
 *
 * One should remember these special cases:
 *
 *	digits(12.3456) == 2    computes with the integer part only
 *	digits(-1234) == 4      computes with the absolute value only
 *	digits(0) == 1          specical case
 *	digits(-0.123) == 1     combination of all of the above
 *
 * given:
 *	q		number to count digits of
 */
long
qdigits(NUMBER *q, ZVALUE base)
{
	long n;			/* number of digits */
	ZVALUE temp;		/* temporary value */

	if (zabsrel(q->num, q->den) < 0)
		return 1;
	if (qisint(q))
		return 1 + zlog(q->num, base);
	zquo(q->num, q->den, &temp, 2);
	n = 1 + zlog(temp, base);
	zfree(temp);
	return n;
}


/*
 * Return the digit at the specified place in the expansion to specified
 * base of a rational number.  The places specified by dpos are numbered from
 * the "units" place just before the "decimal" point, so that negative
 * dpos indicates the (-dpos)th place to the right of the point.
 * Examples: qdigit(1234.5678, 1, 10) = 3, qdigit(1234.5678, -3, 10) = 7.
 * The signs of the number and the base are ignored.
 */
NUMBER *
qdigit(NUMBER *q, ZVALUE dpos, ZVALUE base)
{
	ZVALUE N, D;
	ZVALUE K;
	long k;
	ZVALUE A, B, C;		/* temporary integers */
	NUMBER *res;

	/*
	 * In the first stage, q is expressed as base^k * N/D where
	 *	gcd(D, base) = 1
  	 * K is k as a ZVALUE
	 */
	base.sign = 0;
	if (ziszero(base) || zisunit(base))
			return NULL;
	if (qiszero(q) || (qisint(q) && zisneg(dpos)) ||
			(zge31b(dpos) && !zisneg(dpos)))
		return qlink(&_qzero_);
	k = zfacrem(q->num, base, &N);
	if (k == 0) {
		k = zgcdrem(q->den, base, &D);
		if (k > 0) {
			zequo(q->den, D, &A);
			itoz(k, &K);
			zpowi(base, K, &B);
			zfree(K);
			zequo(B, A, &C);
			zfree(A);
			zfree(B);
			zmul(C, q->num, &N);
			zfree(C);
			k = -k;
		}
		else
			N = q->num;
	}
	if (k >= 0)
		D = q->den;

	itoz(k, &K);
	if (zrel(dpos, K) >= 0) {
		zsub(dpos, K, &A);
		zfree(K);
		zpowi(base, A, &B);
		zfree(A);
		zmul(D, B, &A);
		zfree(B);
		zquo(N, A, &B, 0);
	} else {
		if (zisunit(D)) {
			if (k != 0)
				zfree(N);
			zfree(K);
			if (k < 0)
				zfree(D);
			return qlink(&_qzero_);
		}
		zsub(K, dpos, &A);
		zfree(K);
		zpowermod(base, A, D, &C);
		zfree(A);
		zmod(N, D, &A, 0);
		zmul(C, A, &B);
		zfree(A);
		zfree(C);
		zmod(B, D, &A, 0);
		zfree(B);
		zmodinv(D, base, &B);
		zsub(base, B, &C);
		zfree(B);
		zmul(C, A, &B);
		zfree(C);
	}
	zfree(A);
	if (k != 0)
		zfree(N);
	if (k < 0)
		zfree(D);
	zmod(B, base, &A, 0);
	zfree(B);
	if (ziszero(A)) {
		zfree(A);
		return qlink(&_qzero_);
	}
	if (zisone(A)) {
		zfree(A);
		return qlink(&_qone_);
	}
	res = qalloc();
	res->num = A;
	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(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, 2);
		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(NUMBER *q)
{
	register NUMBER *r;

	if (qisfrac(q)) {
		math_error("Non-integral factorial");
		/*NOTREACHED*/
	}
	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(NUMBER *q)
{
	NUMBER *r;

	if (qisfrac(q)) {
		math_error("Non-integral factorial");
		/*NOTREACHED*/
	}
	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(NUMBER *q)
{
	NUMBER *r;

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


/*
 * Compute the permutation function  q1 * (q1-1) * ... * (q1-q2+1).
 */
NUMBER *
qperm(NUMBER *q1, NUMBER *q2)
{
	NUMBER *r;
	NUMBER *qtmp1, *qtmp2;
	long i;

	if (qisfrac(q2)) {
		math_error("Non-integral second arg for permutation");
		/*NOTREACHED*/
	}
	if (qiszero(q2))
		return qlink(&_qone_);
	if (qisone(q2))
		return qlink(q1);
	if (qisint(q1) && !qisneg(q1)) {
		if (qrel(q2, q1) > 0)
			return qlink(&_qzero_);
		if (qispos(q2)) {
			r = qalloc();
			zperm(q1->num, q2->num, &r->num);
			return r;
		}
	}
	if (zge31b(q2->num)) {
		math_error("Too large arg2 for permutation");
		/*NOTREACHED*/
	}
	i = qtoi(q2);
	if (i > 0) {
		q1 = qlink(q1);
		r = qlink(q1);
		while (--i > 0) {
			qtmp1 = qdec(q1);
			qtmp2 = qmul(r, qtmp1);
			qfree(q1);
			q1 = qtmp1;
			qfree(r);
			r = qtmp2;
		}
		qfree(q1);
		return r;
	}
	i = -i;
	qtmp1 = qinc(q1);
	r = qinv(qtmp1);
	while (--i > 0) {
		qtmp2 = qinc(qtmp1);
		qfree(qtmp1);
		qtmp1 = qqdiv(r, qtmp2);
		qfree(r);
		r = qtmp1;
		qtmp1 = qtmp2;
	}
	qfree(qtmp1);
	return r;
}


/*
 * Compute the combinatorial function  q(q - 1) ...(q - n + 1)/n!
 * n is to be a nonnegative integer
 */
NUMBER *
qcomb(NUMBER *q, NUMBER *n)
{
	NUMBER *r;
	NUMBER *q1, *q2;
	long i, j;
	ZVALUE z;

	if (!qisint(n) || qisneg(n)) {
		math_error("Bad second arg in call to qcomb!");
		/*NOTREACHED*/
	}
	if (qisint(q)) {
		switch (zcomb(q->num, n->num, &z)) {
		case 0:
			return qlink(&_qzero_);
		case 1:
			return qlink(&_qone_);
		case -1:
			return qlink(&_qnegone_);
		case 2:
			return qlink(q);
		case -2:
			return NULL;
		default:
			r = qalloc();
			r->num = z;
			return r;
		}
	}
	if (zge31b(n->num))
		return NULL;
	i = ztoi(n->num);
	q = qlink(q);
	r = qlink(q);
	j = 1;
	while (--i > 0) {
		q1 = qdec(q);
		qfree(q);
		q = q1;
		q2 = qmul(r, q);
		qfree(r);
		r = qdivi(q2, ++j);
		qfree(q2);
	}
	qfree(q);
	return r;
}


/*
 * Compute the Bernoulli number with index n
 * For even positive n, newly calculated values for all even indices up
 * to n are stored in table at B_table.
 */
NUMBER *
qbern(ZVALUE z)
{
	long n, i, k, m, nn, dd;
	NUMBER **p;
	NUMBER *s, *s1, *c, *c1, *t;
	size_t sz;

	if (zisone(z))
		return qlink(&_qneghalf_);

	if (zisodd(z) || z.sign)
		return qlink(&_qzero_);

	if (zge31b(z))
		return NULL;

	n = ztoi(z);

	if (n == 0)

		return qlink(&_qone_);

	m = (n >> 1) - 1;

	if (m < B_num)
		return qlink(B_table[m]);

	if (m >= B_allocnum) {
		k = (m/QALLOCNUM + 1) * QALLOCNUM;
		sz = k * sizeof(NUMBER *);
		if (sz < (size_t) k)
			return NULL;
		if (B_allocnum == 0)
			p = (NUMBER **) malloc(sz);
		else
			p = (NUMBER **) realloc(B_table, sz);
		if (p == NULL)
			return NULL;
		B_allocnum = k;
		B_table = p;
	}
	for (k = B_num; k <= m; k++) {
		nn = 2 * k + 3;
		dd = 1;
		c1 = itoq(nn);
		c = qinv(c1);
		qfree(c1);
		s = qsub(&_qonehalf_, c);
		i = k;
		for (i = 0; i < k; i++) {
			c1 = qmuli(c, nn--);
			qfree(c);
			c = qdivi(c1, dd++);
			qfree(c1);
			c1 = qmuli(c, nn--);
			qfree(c);
			c = qdivi(c1, dd++);
			qfree(c1);
			t = qmul(c, B_table[i]);
			s1 = qsub(s, t);
			qfree(t);
			qfree(s);
			s = s1;
		}
		B_table[k] = s;
		qfree(c);
	}
	B_num = k;
	return qlink(B_table[m]);
}


void
qfreebern(void)
{
	long k;

	if (B_num > 0) {
		for (k = 0; k < B_num; k++)
			qfree(B_table[k]);
		free(B_table);
	}
	B_table = NULL;
	B_allocnum = 0;
	B_num = 0;
}


/*
 * Compute the Euler number with index n = z
 * For even positive n, newly calculated values with all even indices up
 * to n are stored in E_table for later quick access.
 */
NUMBER *
qeuler(ZVALUE z)
{
	long i, k, m, n, nn, dd;
	NUMBER **p;
	NUMBER *s, *s1, *c, *c1, *t;
	size_t sz;


	if (ziszero(z))
		return qlink(&_qone_);
	if (zisodd(z) || zisneg(z))
		return qlink(&_qzero_);
	if (zge31b(z))
		return NULL;
	n = ztoi(z);
	m = (n >> 1) - 1;
	if (m < E_num)
		return qlink(E_table[m]);
	sz = (m + 1) * sizeof(NUMBER *);
	if (sz < (size_t) m + 1)
		return NULL;
	if (E_num)
		p = (NUMBER **) realloc(E_table, sz);
	else
		p = (NUMBER **) malloc(sz);
	if (p == NULL)
		return NULL;
	E_table = p;
	for (k = E_num; k <= m; k++) {
		nn = 2 * k + 2;
		dd = 1;
		c = qlink(&_qone_);
		s = qlink(&_qnegone_);
		i = k;
		for (i = 0; i < k; i++) {
			c1 = qmuli(c, nn--);
			qfree(c);
			c = qdivi(c1, dd++);
			qfree(c1);
			c1 = qmuli(c, nn--);
			qfree(c);
			c = qdivi(c1, dd++);
			qfree(c1);
			t = qmul(c, E_table[i]);
			s1 = qsub(s, t);
			qfree(t);
			qfree(s);
			s = s1;
		}
		E_table[k] = s;
		qfree(c);
	}
	E_num = k;
	return qlink(E_table[m]);
}


void
qfreeeuler(void)
{
	long k;

	if (E_num > 0) {
		for (k = 0; k < E_num; k++)
			qfree(E_table[k]);
		free(E_table);
	}
	E_table = NULL;
	E_num = 0;
}


/*
 * Catalan numbers: catalan(n) = comb(2*n, n)/(n+1)
 * To be called only with integer q
 */
NUMBER *
qcatalan(NUMBER *q)
{
	NUMBER *A, *B;
	NUMBER *res;

	if (qisneg(q))
		return qlink(&_qzero_);
	A = qscale(q, 1);
	B = qcomb(A, q);
	if (B == NULL)
		return NULL;
	qfree(A);
	A = qinc(q);
	res = qqdiv(B, A);
	qfree(A);
	qfree(B);
	return res;
}

/*
 * 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(NUMBER *q1, NUMBER *q2)
{
	if (qisfrac(q1) || qisfrac(q2)) {
		math_error("Non-integral arguments for jacobi");
		/*NOTREACHED*/
	}
	return itoq((long) zjacobi(q1->num, q2->num));
}


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

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


/*
 * Truncate a number to the specified number of decimal places.
 */
NUMBER *
qtrunc(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 = qtenpow(-places);
	r = qmappr(q1, e, 2);
	qfree(e);
	return r;
}




/*
 * Truncate a number to the specified number of binary places.
 * Specifying zero places makes the result identical to qint.
 */
NUMBER *