4.4BSD/usr/src/contrib/calc-1.26.4/matfunc.c

/*
 * Copyright (c) 1993 David I. Bell
 * Permission is granted to use, distribute, or modify this source,
 * provided that this copyright notice remains intact.
 *
 * Extended precision rational arithmetic matrix functions.
 * Matrices can contain arbitrary types of elements.
 */

#include "calc.h"


#if 0
static char *abortmsg = "Calculation aborted";
static char *memmsg = "Not enough memory";
#endif


static void matswaprow(), matsubrow(), matmulrow();
#if 0
static void mataddrow();
#endif

static MATRIX *matident();



/*
 * Add two compatible matrices.
 */
MATRIX *
matadd(m1, m2)
	MATRIX *m1, *m2;
{
	int dim;

	long min1, min2, max1, max2, index;
	VALUE *v1, *v2, *vres;
	MATRIX *res;
	MATRIX tmp;

	if (m1->m_dim != m2->m_dim)
		error("Incompatible matrix dimensions for add");
	tmp.m_dim = m1->m_dim;
	tmp.m_size = m1->m_size;
	for (dim = 0; dim < m1->m_dim; dim++) {
		min1 = m1->m_min[dim];
		max1 = m1->m_max[dim];
		min2 = m2->m_min[dim];
		max2 = m2->m_max[dim];
		if ((min1 && min2 && (min1 != min2)) || ((max1-min1) != (max2-min2)))
			error("Incompatible matrix bounds for add");
		tmp.m_min[dim] = (min1 ? min1 : min2);
		tmp.m_max[dim] = tmp.m_min[dim] + (max1 - min1);
	}
	res = matalloc(m1->m_size);
	*res = tmp;
	v1 = m1->m_table;
	v2 = m2->m_table;
	vres = res->m_table;
	for (index = m1->m_size; index > 0; index--)
		addvalue(v1++, v2++, vres++);
	return res;
}


/*
 * Subtract two compatible matrices.
 */
MATRIX *
matsub(m1, m2)
	MATRIX *m1, *m2;
{
	int dim;
	long min1, min2, max1, max2, index;
	VALUE *v1, *v2, *vres;
	MATRIX *res;
	MATRIX tmp;

	if (m1->m_dim != m2->m_dim)
		error("Incompatible matrix dimensions for sub");
	tmp.m_dim = m1->m_dim;
	tmp.m_size = m1->m_size;
	for (dim = 0; dim < m1->m_dim; dim++) {
		min1 = m1->m_min[dim];
		max1 = m1->m_max[dim];
		min2 = m2->m_min[dim];
		max2 = m2->m_max[dim];
		if ((min1 && min2 && (min1 != min2)) || ((max1-min1) != (max2-min2)))
			error("Incompatible matrix bounds for sub");
		tmp.m_min[dim] = (min1 ? min1 : min2);
		tmp.m_max[dim] = tmp.m_min[dim] + (max1 - min1);
	}
	res = matalloc(m1->m_size);
	*res = tmp;
	v1 = m1->m_table;
	v2 = m2->m_table;
	vres = res->m_table;
	for (index = m1->m_size; index > 0; index--)
		subvalue(v1++, v2++, vres++);
	return res;
}


/*
 * Produce the negative of a matrix.
 */
MATRIX *
matneg(m)
	MATRIX *m;
{
	register VALUE *val, *vres;
	long index;
	MATRIX *res;

	res = matalloc(m->m_size);
	*res = *m;
	val = m->m_table;
	vres = res->m_table;
	for (index = m->m_size; index > 0; index--)
		negvalue(val++, vres++);
	return res;
}


/*
 * Multiply two compatible matrices.
 */
MATRIX *
matmul(m1, m2)
	MATRIX *m1, *m2;
{
	register MATRIX *res;
	long i1, i2, max1, max2, index, maxindex;
	VALUE *v1, *v2;
	VALUE sum, tmp1, tmp2;

	if ((m1->m_dim != 2) || (m2->m_dim != 2))
		error("Matrix dimension must be two for mul");
	if ((m1->m_max[1] - m1->m_min[1]) != (m2->m_max[0] - m2->m_min[0]))
		error("Incompatible bounds for matrix mul");
	max1 = (m1->m_max[0] - m1->m_min[0] + 1);
	max2 = (m2->m_max[1] - m2->m_min[1] + 1);
	maxindex = (m1->m_max[1] - m1->m_min[1] + 1);
	res = matalloc(max1 * max2);
	res->m_dim = 2;
	res->m_min[0] = m1->m_min[0];
	res->m_max[0] = m1->m_max[0];
	res->m_min[1] = m2->m_min[1];
	res->m_max[1] = m2->m_max[1];
	for (i1 = 0; i1 < max1; i1++) {
		for (i2 = 0; i2 < max2; i2++) {
			sum.v_num = qlink(&_qzero_);
			sum.v_type = V_NUM;
			v1 = &m1->m_table[i1 * maxindex];
			v2 = &m2->m_table[i2];
			for (index = 0; index < maxindex; index++) {
				mulvalue(v1, v2, &tmp1);
				addvalue(&sum, &tmp1, &tmp2);
				freevalue(&tmp1);
				freevalue(&sum);
				sum = tmp2;
				v1++;
				v2 += max2;
			}
			index = (i1 * max2) + i2;
			res->m_table[index] = sum;
		}
	}
	return res;
}


/*
 * Square a matrix.
 */
MATRIX *
matsquare(m)
	MATRIX *m;
{
	register MATRIX *res;
	long i1, i2, max, index;
	VALUE *v1, *v2;
	VALUE sum, tmp1, tmp2;

	if (m->m_dim != 2)
		error("Matrix dimension must be two for square");
	if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1]))
		error("Squaring non-square matrix");
	max = (m->m_max[0] - m->m_min[0] + 1);
	res = matalloc(max * max);
	res->m_dim = 2;
	res->m_min[0] = m->m_min[0];
	res->m_max[0] = m->m_max[0];
	res->m_min[1] = m->m_min[1];
	res->m_max[1] = m->m_max[1];
	for (i1 = 0; i1 < max; i1++) {
		for (i2 = 0; i2 < max; i2++) {
			sum.v_num = qlink(&_qzero_);
			sum.v_type = V_NUM;
			v1 = &m->m_table[i1 * max];
			v2 = &m->m_table[i2];
			for (index = 0; index < max; index++) {
				mulvalue(v1, v2, &tmp1);
				addvalue(&sum, &tmp1, &tmp2);
				freevalue(&tmp1);
				freevalue(&sum);
				sum = tmp2;
				v1++;
				v2 += max;
			}
			index = (i1 * max) + i2;
			res->m_table[index] = sum;
		}
	}
	return res;
}


/*
 * Compute the result of raising a square matrix to an integer power.
 * Negative powers mean the positive power of the inverse.
 * Note: This calculation could someday be improved for large powers
 * by using the characteristic polynomial of the matrix.
 */
MATRIX *
matpowi(m, q)
	MATRIX *m;		/* matrix to be raised */
	NUMBER *q;		/* power to raise it to */
{
	MATRIX *res, *tmp;
	long power;		/* power to raise to */
	unsigned long bit;	/* current bit value */

	if (m->m_dim != 2)
		error("Matrix dimension must be two for power");
	if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1]))
		error("Raising non-square matrix to a power");
	if (qisfrac(q))
		error("Raising matrix to non-integral power");
	if (isbig(q->num))
		error("Raising matrix to very large power");
	power = (istiny(q->num) ? z1tol(q->num) : z2tol(q->num));
	if (qisneg(q))
		power = -power;
	/*
	 * Handle some low powers specially
	 */
	if ((power <= 4) && (power >= -2)) {
		switch ((int) power) {
			case 0:
				return matident(m);
			case 1:
				return matcopy(m);
			case -1:
				return matinv(m);
			case 2:
				return matsquare(m);
			case -2:
				tmp = matinv(m);
				res = matsquare(tmp);
				matfree(tmp);
				return res;
			case 3:
				tmp = matsquare(m);
				res = matmul(m, tmp);
				matfree(tmp);
				return res;
			case 4:
				tmp = matsquare(m);
				res = matsquare(tmp);
				matfree(tmp);
				return res;
		}
	}
	if (power < 0) {
		m = matinv(m);
		power = -power;
	}
	/*
	 * Compute the power by squaring and multiplying.
	 * This uses the left to right method of power raising.
	 */
	bit = TOPFULL;
	while ((bit & power) == 0)
		bit >>= 1L;
	bit >>= 1L;
	res = matsquare(m);
	if (bit & power) {
		tmp = matmul(res, m);
		matfree(res);
		res = tmp;
	}
	bit >>= 1L;
	while (bit) {
		tmp = matsquare(res);
		matfree(res);
		res = tmp;
		if (bit & power) {
			tmp = matmul(res, m);
			matfree(res);
			res = tmp;
		}
		bit >>= 1L;
	}
	if (qisneg(q))
		matfree(m);
	return res;
}


/*
 * Calculate the cross product of two one dimensional matrices each
 * with three components.
 *	m3 = matcross(m1, m2);
 */
MATRIX *
matcross(m1, m2)
	MATRIX *m1, *m2;
{
	MATRIX *res;
	VALUE *v1, *v2, *vr;
	VALUE tmp1, tmp2;

	if ((m1->m_dim != 1) || (m2->m_dim != 1))
		error("Matrix not 1d for cross product");
	if ((m1->m_size != 3) || (m2->m_size != 3))
		error("Matrix not size 3 for cross product");
	res = matalloc(3L);
	res->m_dim = 1;
	res->m_min[0] = 0;
	res->m_max[0] = 2;
	v1 = m1->m_table;
	v2 = m2->m_table;
	vr = res->m_table;
	mulvalue(v1 + 1, v2 + 2, &tmp1);
	mulvalue(v1 + 2, v2 + 1, &tmp2);
	subvalue(&tmp1, &tmp2, vr + 0);
	freevalue(&tmp1);
	freevalue(&tmp2);
	mulvalue(v1 + 2, v2 + 0, &tmp1);
	mulvalue(v1 + 0, v2 + 2, &tmp2);
	subvalue(&tmp1, &tmp2, vr + 1);
	freevalue(&tmp1);
	freevalue(&tmp2);
	mulvalue(v1 + 0, v2 + 1, &tmp1);
	mulvalue(v1 + 1, v2 + 0, &tmp2);
	subvalue(&tmp1, &tmp2, vr + 2);
	freevalue(&tmp1);
	freevalue(&tmp2);
	return res;
}


/*
 * Return the dot product of two matrices.
 *	result = matdot(m1, m2);
 */
VALUE
matdot(m1, m2)
	MATRIX *m1, *m2;
{
	VALUE *v1, *v2;
	VALUE result, tmp1, tmp2;
	long len;

	if ((m1->m_dim != 1) || (m2->m_dim != 1))
		error("Matrix not 1d for dot product");
	if (m1->m_size != m2->m_size)
		error("Incompatible matrix sizes for dot product");
	v1 = m1->m_table;
	v2 = m2->m_table;
	mulvalue(v1, v2, &result);
	len = m1->m_size;
	while (--len > 0) {
		mulvalue(++v1, ++v2, &tmp1);
		addvalue(&result, &tmp1, &tmp2);
		freevalue(&tmp1);
		freevalue(&result);
		result = tmp2;
	}
	return result;
}


/*
 * Scale the elements of a matrix by a specified power of two.
 */
MATRIX *
matscale(m, n)
	MATRIX *m;		/* matrix to be scaled */
	long n;
{
	register VALUE *val, *vres;
	VALUE num;
	long index;
	MATRIX *res;		/* resulting matrix */

	if (n == 0)
		return matcopy(m);
	num.v_type = V_NUM;
	num.v_num = itoq(n);
	res = matalloc(m->m_size);
	*res = *m;
	val = m->m_table;
	vres = res->m_table;
	for (index = m->m_size; index > 0; index--)
		scalevalue(val++, &num, vres++);
	qfree(num.v_num);
	return res;
}


/*
 * Shift the elements of a matrix by the specified number of bits.
 * Positive shift means leftwards, negative shift rightwards.
 */
MATRIX *
matshift(m, n)
	MATRIX *m;		/* matrix to be scaled */
	long n;
{
	register VALUE *val, *vres;
	VALUE num;
	long index;
	MATRIX *res;		/* resulting matrix */

	if (n == 0)
		return matcopy(m);
	num.v_type = V_NUM;
	num.v_num = itoq(n);
	res = matalloc(m->m_size);
	*res = *m;
	val = m->m_table;
	vres = res->m_table;
	for (index = m->m_size; index > 0; index--)
		shiftvalue(val++, &num, FALSE, vres++);
	qfree(num.v_num);
	return res;
}


/*
 * Multiply the elements of a matrix by a specified value.
 */
MATRIX *
matmulval(m, vp)
	MATRIX *m;		/* matrix to be multiplied */
	VALUE *vp;		/* value to multiply by */
{
	register VALUE *val, *vres;
	long index;
	MATRIX *res;

	res = matalloc(m->m_size);
	*res = *m;
	val = m->m_table;
	vres = res->m_table;
	for (index = m->m_size; index > 0; index--)
		mulvalue(val++, vp, vres++);
	return res;
}


/*
 * Divide the elements of a matrix by a specified value, keeping
 * only the integer quotient.
 */
MATRIX *
matquoval(m, vp)
	MATRIX *m;		/* matrix to be divided */
	VALUE *vp;		/* value to divide by */
{
	register VALUE *val, *vres;
	long index;
	MATRIX *res;

	if ((vp->v_type == V_NUM) && qiszero(vp->v_num))
		error("Division by zero");
	res = matalloc(m->m_size);
	*res = *m;
	val = m->m_table;
	vres = res->m_table;
	for (index = m->m_size; index > 0; index--)
		quovalue(val++, vp, vres++);
	return res;
}


/*
 * Divide the elements of a matrix by a specified value, keeping
 * only the remainder of the division.
 */
MATRIX *
matmodval(m, vp)
	MATRIX *m;		/* matrix to be divided */
	VALUE *vp;		/* value to divide by */
{
	register VALUE *val, *vres;
	long index;
	MATRIX *res;

	if ((vp->v_type == V_NUM) && qiszero(vp->v_num))
		error("Division by zero");
	res = matalloc(m->m_size);
	*res = *m;
	val = m->m_table;
	vres = res->m_table;
	for (index = m->m_size; index > 0; index--)
		modvalue(val++, vp, vres++);
	return res;
}


MATRIX *
mattrans(m)
	MATRIX *m;			/* matrix to be transposed */
{
	register VALUE *v1, *v2;	/* current values */
	long rows, cols;		/* rows and columns in new matrix */
	long row, col;			/* current row and column */
	MATRIX *res;

	if (m->m_dim != 2)
		error("Matrix dimension must be two for transpose");
	res = matalloc(m->m_size);
	res->m_dim = 2;
	res->m_min[0] = m->m_min[1];
	res->m_max[0] = m->m_max[1];
	res->m_min[1] = m->m_min[0];
	res->m_max[1] = m->m_max[0];
	rows = (m->m_max[1] - m->m_min[1] + 1);
	cols = (m->m_max[0] - m->m_min[0] + 1);
	v1 = res->m_table;
	for (row = 0; row < rows; row++) {
		v2 = &m->m_table[row];
		for (col = 0; col < cols; col++) {
			copyvalue(v2, v1);
			v1++;
			v2 += rows;
		}
	}
	return res;
}


/*
 * Produce a matrix with values all of which are conjugated.
 */
MATRIX *
matconj(m)
	MATRIX *m;
{
	register VALUE *val, *vres;
	long index;
	MATRIX *res;

	res = matalloc(m->m_size);
	*res = *m;
	val = m->m_table;
	vres = res->m_table;
	for (index = m->m_size; index > 0; index--)
		conjvalue(val++, vres++);
	return res;
}


/*
 * Produce a matrix with values all of which have been rounded to the
 * specified number of decimal places.
 */
MATRIX *
matround(m, places)
	MATRIX *m;
	long places;
{
	register VALUE *val, *vres;
	VALUE tmp;
	long index;
	MATRIX *res;

	if (places < 0)
		error("Negative number of places for matround");
	res = matalloc(m->m_size);
	*res = *m;
	tmp.v_type = V_INT;
	tmp.v_int = places;
	val = m->m_table;
	vres = res->m_table;
	for (index = m->m_size; index > 0; index--)
		roundvalue(val++, &tmp, vres++);
	return res;
}


/*
 * Produce a matrix with values all of which have been rounded to the
 * specified number of binary places.
 */
MATRIX *
matbround(m, places)
	MATRIX *m;
	long places;
{
	register VALUE *val, *vres;
	VALUE tmp;
	long index;
	MATRIX *res;

	if (places < 0)
		error("Negative number of places for matbround");
	res = matalloc(m->m_size);
	*res = *m;
	tmp.v_type = V_INT;
	tmp.v_int = places;
	val = m->m_table;
	vres = res->m_table;
	for (index = m->m_size; index > 0; index--)
		broundvalue(val++, &tmp, vres++);
	return res;
}


/*
 * Produce a matrix with values all of which have been truncated to integers.
 */
MATRIX *
matint(m)
	MATRIX *m;
{
	register VALUE *val, *vres;
	long index;
	MATRIX *res;

	res = matalloc(m->m_size);
	*res = *m;
	val = m->m_table;
	vres = res->m_table;
	for (index = m->m_size; index > 0; index--)
		intvalue(val++, vres++);
	return res;
}


/*
 * Produce a matrix with values all of which have only the fraction part left.
 */
MATRIX *
matfrac(m)
	MATRIX *m;
{
	register VALUE *val, *vres;
	long index;
	MATRIX *res;

	res = matalloc(m->m_size);
	*res = *m;
	val = m->m_table;
	vres = res->m_table;
	for (index = m->m_size; index > 0; index--)
		fracvalue(val++, vres++);
	return res;
}


/*
 * Search a matrix for the specified value, starting with the specified index.
 * Returns the index of the found value, or -1 if the value was not found.
 */
long
matsearch(m, vp, index)
	MATRIX *m;
	VALUE *vp;
	long index;
{
	register VALUE *val;

	if (index < 0)
		index = 0;
	val = &m->m_table[index];
	while (index < m->m_size) {
		if (!comparevalue(vp, val))
			return index;
		index++;
		val++;
	}
	return -1;
}


/*
 * Search a matrix backwards for the specified value, starting with the
 * specified index.  Returns the index of the found value, or -1 if the
 * value was not found.
 */
long
matrsearch(m, vp, index)
	MATRIX *m;
	VALUE *vp;
	long index;
{
	register VALUE *val;

	if (index >= m->m_size)
		index = m->m_size - 1;
	val = &m->m_table[index];
	while (index >= 0) {
		if (!comparevalue(vp, val))
			return index;
		index--;
		val--;
	}
	return -1;
}


/*
 * Fill all of the elements of a matrix with one of two specified values.
 * All entries are filled with the first specified value, except that if
 * the matrix is square and the second value pointer is non-NULL, then
 * all diagonal entries are filled with the second value.  This routine
 * affects the supplied matrix directly, and doesn't return a copy.
 */
void
matfill(m, v1, v2)
	MATRIX *m;		/* matrix to be filled */
	VALUE *v1;		/* value to fill most of matrix with */
	VALUE *v2;		/* value for diagonal entries (or NULL) */
{
	register VALUE *val;
	long row, col;
	long rows;
	long index;

	if (v2 && ((m->m_dim != 2) ||
		((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1]))))
			error("Filling diagonals of non-square matrix");
	val = m->m_table;
	for (index = m->m_size; index > 0; index--)
		freevalue(val++);
	val = m->m_table;
	if (v2 == NULL) {
		for (index = m->m_size; index > 0; index--)
			copyvalue(v1, val++);
		return;
	}
	rows = m->m_max[0] - m->m_min[0] + 1;
	for (row = 0; row < rows; row++) {
		for (col = 0; col < rows; col++) {
			copyvalue(((row != col) ? v1 : v2), val++);
		}
	}
}


/*
 * Set a copy of a square matrix to the identity matrix.
 */
static MATRIX *
matident(m)
	MATRIX *m;
{
	register VALUE *val;	/* current value */
	long row, col;		/* current row and column */
	long rows;		/* number of rows */
	MATRIX *res;		/* resulting matrix */

	if (m->m_dim != 2)
		error("Matrix dimension must be two for setting to identity");
	if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1]))
		error("Matrix must be square for setting to identity");
	res = matalloc(m->m_size);
	*res = *m;
	val = res->m_table;
	rows = (res->m_max[0] - res->m_min[0] + 1);
	for (row = 0; row < rows; row++) {
		for (col = 0; col < rows; col++) {
			val->v_type = V_NUM;
			val->v_num = ((row == col) ? qlink(&_qone_) : qlink(&_qzero_));
			val++;
		}
	}
	return res;
}


/*
 * Calculate the inverse of a matrix if it exists.
 * This is done by using transformations on the supplied matrix to convert
 * it to the identity matrix, and simultaneously applying the same set of
 * transformations to the identity matrix.
 */
MATRIX *
matinv(m)
	MATRIX *m;
{
	MATRIX *res;		/* matrix to become the inverse */
	long rows;		/* number of rows */
	long cur;		/* current row being worked on */
	long row, col;		/* temp row and column values */
	VALUE *val;		/* current value in matrix*/
	VALUE mulval;		/* value to multiply rows by */
	VALUE tmpval;		/* temporary value */

	if (m->m_dim != 2)
		error("Matrix dimension must be two for inverse");
	if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1]))
		error("Inverting non-square matrix");
	/*
	 * Begin by creating the identity matrix with the same attributes.
	 */
	res = matalloc(m->m_size);
	*res = *m;
	rows = (m->m_max[0] - m->m_min[0] + 1);
	val = res->m_table;
	for (row = 0; row < rows; row++) {
		for (col = 0; col < rows; col++) {
			if (row == col)
				val->v_num = qlink(&_qone_);
			else
				val->v_num = qlink(&_qzero_);
			val->v_type = V_NUM;
			val++;
		}
	}
	/*
	 * Now loop over each row, and eliminate all entries in the
	 * corresponding column by using row operations.  Do the same
	 * operations on the resulting matrix.  Copy the original matrix
	 * so that we don't destroy it.
	 */
	m = matcopy(m);
	for (cur = 0; cur < rows; cur++) {
		/*
		 * Find the first nonzero value in the rest of the column
		 * downwards from [cur,cur].  If there is no such value, then
		 * the matrix is not invertible.  If the first nonzero entry
		 * is not the current row, then swap the two rows to make the
		 * current one nonzero.
		 */
		row = cur;
		val = &m->m_table[(row * rows) + row];
		while (testvalue(val) == 0) {
			if (++row >= rows) {
				matfree(m);
				matfree(res);
				error("Matrix is not invertible");
			}
			val += rows;
		}
		invertvalue(val, &mulval);
		if (row != cur) {
			matswaprow(m, row, cur);
			matswaprow(res, row, cur);
		}
		/*
		 * Now for every other nonzero entry in the current column, subtract
		 * the appropriate multiple of the current row to force that entry
		 * to become zero.
		 */
		val = &m->m_table[cur];
		for (row = 0; row < rows; row++, val += rows) {
			if ((row == cur) || (testvalue(val) == 0))
				continue;
			mulvalue(val, &mulval, &tmpval);
			matsubrow(m, row, cur, &tmpval);
			matsubrow(res, row, cur, &tmpval);
			freevalue(&tmpval);
		}
		freevalue(&mulval);
	}
	/*
	 * Now the original matrix has nonzero entries only on its main diagonal.
	 * Scale the rows of the result matrix by the inverse of those entries.
	 */
	val = m->m_table;
	for (row = 0; row < rows; row++) {
		if ((val->v_type != V_NUM) || !qisone(val->v_num)) {
			invertvalue(val, &mulval);
			matmulrow(res, row, &mulval);
			freevalue(&mulval);
		}
		val += (rows + 1);
	}
	matfree(m);
	return res;
}


/*
 * Calculate the determinant of a square matrix.
 * This is done using row operations to create an upper-diagonal matrix.
 */
VALUE
matdet(m)
	MATRIX *m;
{
	long rows;		/* number of rows */
	long cur;		/* current row being worked on */
	long row;		/* temp row values */
	int neg;		/* whether to negate determinant */
	VALUE *val;		/* current value */
	VALUE mulval, tmpval;	/* other values */

	if (m->m_dim != 2)
		error("Matrix dimension must be two for determinant");
	if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1]))
		error("Non-square matrix for determinant");
	/*
	 * Loop over each row, and eliminate all lower entries in the
	 * corresponding column by using row operations.  Copy the original
	 * matrix so that we don't destroy it.
	 */
	neg = 0;
	m = matcopy(m);
	rows = (m->m_max[0] - m->m_min[0] + 1);
	for (cur = 0; cur < rows; cur++) {
		/*
		 * Find the first nonzero value in the rest of the column
		 * downwards from [cur,cur].  If there is no such value, then
		 * the determinant is zero.  If the first nonzero entry is not
		 * the current row, then swap the two rows to make the current
		 * one nonzero, and remember that the determinant changes sign.
		 */
		row = cur;
		val = &m->m_table[(row * rows) + row];
		while (testvalue(val) == 0) {
			if (++row >= rows) {
				matfree(m);
				mulval.v_type = V_NUM;
				mulval.v_num = qlink(&_qzero_);
				return mulval;
			}
			val += rows;
		}
		invertvalue(val, &mulval);
		if (row != cur) {
			matswaprow(m, row, cur);
			neg = !neg;
		}
		/*
		 * Now for every other nonzero entry lower down in the current column,
		 * subtract the appropriate multiple of the current row to force that
		 * entry to become zero.
		 */
		row = cur + 1;
		val = &m->m_table[(row * rows) + cur];
		for (; row < rows; row++, val += rows) {
			if (testvalue(val) == 0)
				continue;
			mulvalue(val, &mulval, &tmpval);
			matsubrow(m, row, cur, &tmpval);
			freevalue(&tmpval);
		}
		freevalue(&mulval);
	}
	/*
	 * Now the matrix is upper-diagonal, and the determinant is the
	 * product of the main diagonal entries, and is possibly negated.
	 */
	val = m->m_table;
	mulval.v_type = V_NUM;
	mulval.v_num = qlink(&_qone_);
	for (row = 0; row < rows; row++) {
		mulvalue(&mulval, val, &tmpval);
		freevalue(&mulval);
		mulval = tmpval;
		val += (rows + 1);
	}
	matfree(m);
	if (neg) {
		negvalue(&mulval, &tmpval);
		freevalue(&mulval);
		return tmpval;
	}
	return mulval;
}


/*
 * Local utility routine to swap two rows of a square matrix.
 * No checks are made to verify the legality of the arguments.
 */
static void
matswaprow(m, r1, r2)
	MATRIX *m;
	long r1, r2;
{
	register VALUE *v1, *v2;
	register long rows;
	VALUE tmp;

	if (r1 == r2)
		return;
	rows = (m->m_max[0] - m->m_min[0] + 1);
	v1 = &m->m_table[r1 * rows];
	v2 = &m->m_table[r2 * rows];
	while (rows-- > 0) {
		tmp = *v1;
		*v1 = *v2;
		*v2 = tmp;
		v1++;
		v2++;
	}
}


/*
 * Local utility routine to subtract a multiple of one row to another one.
 * The row to be changed is oprow, the row to be subtracted is baserow.
 * No checks are made to verify the legality of the arguments.
 */
static void
matsubrow(m, oprow, baserow, mulval)
	MATRIX *m;
	long oprow, baserow;
	VALUE *mulval;
{
	register VALUE *vop, *vbase;
	register long entries;
	VALUE tmp1, tmp2;

	entries = (m->m_max[0] - m->m_min[0] + 1);
	vop = &m->m_table[oprow * entries];
	vbase = &m->m_table[baserow * entries];
	while (entries-- > 0) {
		mulvalue(vbase, mulval, &tmp1);
		subvalue(vop, &tmp1, &tmp2);
		freevalue(&tmp1);
		freevalue(vop);
		*vop = tmp2;
		vop++;
		vbase++;
	}
}


#if 0
/*
 * Local utility routine to add one row to another one.
 * No checks are made to verify the legality of the arguments.
 */
static void
mataddrow(m, r1, r2)
	MATRIX *m;
	long r1;	/* row to be added into */
	long r2;	/* row to add */
{
	register VALUE *v1, *v2;
	register long rows;
	VALUE tmp;

	rows = (m->m_max[0] - m->m_min[0] + 1);
	v1 = &m->m_table[r1 * rows];
	v2 = &m->m_table[r2 * rows];
	while (rows-- > 0) {
		addvalue(v1, v2, &tmp);
		freevalue(v1);
		*v1 = tmp;
		v1++;
		v2++;
	}
}
#endif


/*
 * Local utility routine to multiply a row by a specified number.
 * No checks are made to verify the legality of the arguments.
 */
static void
matmulrow(m, row, mulval)
	MATRIX *m;
	long row;
	VALUE *mulval;
{
	register VALUE *val;
	register long rows;
	VALUE tmp;

	rows = (m->m_max[0] - m->m_min[0] + 1);
	val = &m->m_table[row * rows];
	while (rows-- > 0) {
		mulvalue(val, mulval, &tmp);
		freevalue(val);
		*val = tmp;
		val++;
	}
}


/*
 * Make a full copy of a matrix.
 */
MATRIX *
matcopy(m)
	MATRIX *m;
{
	MATRIX *res;
	register VALUE *v1, *v2;
	register long i;

	res = matalloc(m->m_size);
	*res = *m;
	v1 = m->m_table;
	v2 = res->m_table;
	i = m->m_size;
	while (i-- > 0) {
		if (v1->v_type == V_NUM) {
			v2->v_type = V_NUM;
			v2->v_num = qlink(v1->v_num);
		} else
			copyvalue(v1, v2);
		v1++;
		v2++;
	}
	return res;
}


/*
 * Allocate a matrix with the specified number of elements.
 */
MATRIX *
matalloc(size)
	long size;
{
	MATRIX *m;

	m = (MATRIX *) malloc(matsize(size));
	if (m == NULL)
		error("Cannot get memory to allocate matrix of size %d", size);
	m->m_size = size;
	return m;
}


/*
 * Free a matrix, along with all of its element values.
 */
void
matfree(m)
	MATRIX *m;
{
	register VALUE *vp;
	register long i;

	vp = m->m_table;
	i = m->m_size;
	while (i-- > 0) {
		if (vp->v_type == V_NUM) {
			vp->v_type = V_NULL;
			qfree(vp->v_num);
		} else
			freevalue(vp);
		vp++;
	}
	free(m);
}


/*
 * Test whether a matrix has any nonzero values.
 * Returns TRUE if so.
 */
BOOL
mattest(m)
	MATRIX *m;
{
	register VALUE *vp;
	register long i;

	vp = m->m_table;
	i = m->m_size;
	while (i-- > 0) {
		if ((vp->v_type != V_NUM) || (!qiszero(vp->v_num)))
			return TRUE;
		vp++;
	}
	return FALSE;
}


/*
 * Test whether or not two matrices are equal.
 * Equality is determined by the shape and values of the matrices,
 * but not by their index bounds.  Returns TRUE if they differ.
 */
BOOL
matcmp(m1, m2)
	MATRIX *m1, *m2;
{
	VALUE *v1, *v2;
	long i;

	if (m1 == m2)
		return FALSE;
	if ((m1->m_dim != m2->m_dim) || (m1->m_size != m2->m_size))
		return TRUE;
	for (i = 0; i < m1->m_dim; i++) {
		if ((m1->m_max[i] - m1->m_min[i]) != (m2->m_max[i] - m2->m_min[i]))
		return TRUE;
	}
	v1 = m1->m_table;
	v2 = m2->m_table;
	i = m1->m_size;
	while (i-- > 0) {
		if (comparevalue(v1++, v2++))
			return TRUE;
	}
	return FALSE;
}


#if 0
/*
 * Test whether or not a matrix is the identity matrix.
 * Returns TRUE if so.
 */
BOOL
matisident(m)
	MATRIX *m;
{
	register VALUE *val;	/* current value */
	long row, col;		/* row and column numbers */

	if ((m->m_dim != 2) ||
		((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])))
			return FALSE;
	val = m->m_table;
	for (row = 0; row < m->m_size; row++) {
		for (col = 0; col < m->m_size; col++) {
			if (val->v_type != V_NUM)
				return FALSE;
			if (row == col) {
				if (!qisone(val->v_num))
					return FALSE;
			} else {
				if (!qiszero(val->v_num))
					return FALSE;
			}
			val++;
		}
	}
	return TRUE;
}
#endif


/*
 * Print a matrix and possibly few of its elements.
 * The argument supplied specifies how many elements to allow printing.
 * If any elements are printed, they are printed in short form.
 */
void
matprint(m, max_print)
	MATRIX *m;
	long max_print;
{
	VALUE *vp;
	long fullsize, count, index, num;
	int dim, i;
	char *msg;
	long sizes[MAXDIM];

	dim = m->m_dim;
	fullsize = 1;
	for (i = dim - 1; i >= 0; i--) {
		sizes[i] = fullsize;
		fullsize *= (m->m_max[i] - m->m_min[i] + 1);
	}
	msg = ((max_print > 0) ? "\nmat [" : "mat [");
	for (i = 0; i < dim; i++) {
		if (m->m_min[i])
			math_fmt("%s%ld:%ld", msg, m->m_min[i], m->m_max[i]);
		else
			math_fmt("%s%ld", msg, m->m_max[i] + 1);
		msg = ",";
	}
	if (max_print > fullsize)
		max_print = fullsize;
	vp = m->m_table;
	count = 0;
	for (index = 0; index < fullsize; index++) {
		if ((vp->v_type != V_NUM) || !qiszero(vp->v_num))
			count++;
		vp++;
	}
	math_fmt("] (%ld element%s, %ld nonzero)",
		fullsize, (fullsize == 1) ? "" : "s", count);
	if (max_print <= 0)
		return;

	/*
	 * Now print the first few elements of the matrix in short
	 * and unambigous format.
	 */
	math_str(":\n");
	vp = m->m_table;
	for (index = 0; index < max_print; index++) {
		msg = "  [";
		num = index;
		for (i = 0; i < dim; i++) {
			math_fmt("%s%ld", msg, m->m_min[i] + (num / sizes[i]));
			num %= sizes[i];
			msg = ",";
		}
		math_str("] = ");
		printvalue(vp++, PRINT_SHORT | PRINT_UNAMBIG);
		math_str("\n");
	}
	if (max_print < fullsize)
		math_str("  ...\n");
}

/* END CODE */