calc/matfunc.c

/*
 * matfunc - extended precision rational arithmetic matrix functions
 *
 * Copyright (C) 1999-2007,2021-2023  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.
 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 *
 * Under source code control:	1990/02/15 01:48:18
 * File existed as early as:	before 1990
 *
 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
 */

/*
 * Extended precision rational arithmetic matrix functions.
 * Matrices can contain arbitrary types of elements.
 */

#include "alloc.h"
#include "value.h"
#include "zrand.h"

#include "have_unused.h"


#include "errtbl.h"
#include "banned.h"	/* include after system header <> includes */


E_FUNC long irand(long s);

S_FUNC void matswaprow(MATRIX *m, long r1, long r2);
S_FUNC void matsubrow(MATRIX *m, long oprow, long baserow, VALUE *mulval);
S_FUNC void matmulrow(MATRIX *m, long row, VALUE *mulval);
S_FUNC MATRIX *matident(MATRIX *m);



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

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

	if (m1->m_dim != m2->m_dim) {
		math_error("Incompatible matrix dimensions for add");
		not_reached();
	}
	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))) {
			math_error("Incompatible matrix bounds for add");
			not_reached();
		}
		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(MATRIX *m1, MATRIX *m2)
{
	int dim;
	long min1, min2, max1, max2, index;
	VALUE *v1, *v2, *vres;
	MATRIX *res;
	MATRIX tmp;

	if (m1->m_dim != m2->m_dim) {
		math_error("Incompatible matrix dimensions for sub");
		not_reached();
	}
	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))) {
			math_error("Incompatible matrix bounds for sub");
			not_reached();
		}
		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(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(MATRIX *m1, MATRIX *m2)
{
	register MATRIX *res;
	long i1, i2, max1, max2, index, maxindex;
	VALUE *v1, *v2, *vres;
	VALUE sum, tmp1, tmp2;

	if (m1->m_dim == 0) {
		i2 = m2->m_size;
		v2 = m2->m_table;
		res = matalloc(i2);
		*res = *m2;
		vres = res->m_table;
		while (i2-- > 0)
			mulvalue(m1->m_table, v2++, vres++);
		return res;
	}
	if (m2->m_dim == 0) {
		i1 = m1->m_size;
		v1 = m1->m_table;
		res = matalloc(i1);
		*res = *m1;
		vres = res->m_table;
		while (i1-- > 0)
			mulvalue(v1++, m2->m_table, vres++);
		return res;
	}
	if (m1->m_dim == 1 && m2->m_dim == 1) {
		if (m1->m_max[0]-m1->m_min[0] != m2->m_max[0]-m2->m_min[0]) {
			math_error("Incompatible bounds for 1D * 1D  matmul");
			not_reached();
		}
		res = matalloc(m1->m_size);
		*res = *m1;
		v1 = m1->m_table;
		v2 = m2->m_table;
		vres = res->m_table;
		for (index = m1->m_size; index > 0; index--)
			mulvalue(v1++, v2++, vres++);
		return res;
	}
	if (m1->m_dim == 1 && m2->m_dim == 2) {
		if (m1->m_max[0]-m1->m_min[0] != m2->m_max[0]-m2->m_min[0]) {
			math_error("Incompatible bounds for 1D * 2D matmul");
			not_reached();
		}
		res = matalloc(m2->m_size);
		*res = *m2;
		i1 = m1->m_max[0] - m1->m_min[0] + 1;
		max2 = m2->m_max[1] - m2->m_min[1] + 1;
		v1 = m1->m_table;
		v2 = m2->m_table;
		vres = res->m_table;
		while (i1-- > 0) {
			i2 = max2;
			while (i2-- > 0)
				mulvalue(v1, v2++, vres++);
			v1++;
		}
		return res;
	}
	if (m1->m_dim == 2 && m2->m_dim == 1) {
		if (m1->m_max[1]-m1->m_min[1] != m2->m_max[0]-m2->m_min[0]) {
			math_error("Incompatible bounds for 2D * 1D matmul");
			not_reached();
		}
		res = matalloc(m1->m_size);
		*res = *m1;
		i1 = m1->m_max[0] - m1->m_min[0] + 1;
		max1 = m1->m_max[1] - m1->m_min[1] + 1;
		v1 = m1->m_table;
		vres = res->m_table;
		while (i1-- > 0) {
			v2 = m2->m_table;
			i2 = max1;
			while (i2-- > 0)
				mulvalue(v1++, v2++, vres++);
		}
		return res;
	}

	if ((m1->m_dim != 2) || (m2->m_dim != 2)) {
		math_error("Matrix dimensions not compatible for mul");
		not_reached();
	}
	if ((m1->m_max[1]-m1->m_min[1]) != (m2->m_max[0]-m2->m_min[0])) {
		math_error("Incompatible bounds for 2D * 2D matrix mul");
		not_reached();
	}
	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_type = V_NULL;
			sum.v_subtype = V_NOSUBTYPE;
			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++;
				if (index+1 < maxindex)
					v2 += max2;
			}
			index = (i1 * max2) + i2;
			res->m_table[index] = sum;
		}
	}
	return res;
}


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

	if (m->m_dim < 2) {
		res = matalloc(m->m_size);
		*res = *m;
		v1 = m->m_table;
		v2 = res->m_table;
		for (index = m->m_size; index > 0; index--)
			squarevalue(v1++, v2++);
		return res;
	}
	if (m->m_dim != 2) {
		math_error("Matrix dimension exceeds two for square");
		not_reached();
	}
	if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])) {
		math_error("Squaring non-square matrix");
		not_reached();
	}
	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_type = V_NULL;
			sum.v_subtype = V_NOSUBTYPE;
			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 matrix to an integer power if
 * dimension <= 2 and for dimension == 2, the matrix is square.
 * 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.
 *
 * given:
 *	m		matrix to be raised
 *	q		power to raise it to
 */
MATRIX *
matpowi(MATRIX *m, NUMBER *q)
{
	MATRIX *res, *tmp;
	long power;		/* power to raise to */
	FULL bit;		/* current bit value */

	if (m->m_dim > 2) {
		math_error("Matrix dimension greater than 2 for power");
		not_reached();
	}
	if (m->m_dim == 2 && (m->m_max[0] - m->m_min[0] !=
			 m->m_max[1] - m->m_min[1])) {
		math_error("Raising non-square 2D matrix to a power");
		not_reached();
	}
	if (qisfrac(q)) {
		math_error("Raising matrix to non-integral power");
		not_reached();
	}
	if (zge31b(q->num)) {
		math_error("Raising matrix to very large power");
		not_reached();
	}
	power = ztolong(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(MATRIX *m1, MATRIX *m2)
{
	MATRIX *res;
	VALUE *v1, *v2, *vr;
	VALUE tmp1, tmp2;

	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(MATRIX *m1, MATRIX *m2)
{
	VALUE *v1, *v2;
	VALUE result, tmp1, tmp2;
	long len;

	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.
 *
 * given:
 *	m		matrix to be scaled
 *	n		scale factor
 */
MATRIX *
matscale(MATRIX *m, long n)
{
	register VALUE *val, *vres;
	VALUE temp;
	long index;
	MATRIX *res;		/* resulting matrix */

	if (n == 0)
		return matcopy(m);
	temp.v_type = V_NUM;
	temp.v_subtype = V_NOSUBTYPE;
	temp.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++, &temp, vres++);
	qfree(temp.v_num);
	return res;
}


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

	if (n == 0)
		return matcopy(m);
	temp.v_type = V_NUM;
	temp.v_subtype = V_NOSUBTYPE;
	temp.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++, &temp, false, vres++);
	qfree(temp.v_num);
	return res;
}


/*
 * Multiply the elements of a matrix by a specified value.
 *
 * given:
 *	m		matrix to be multiplied
 *	vp		value to multiply by
 */
MATRIX *
matmulval(MATRIX *m, VALUE *vp)
{
	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.
 *
 * given:
 *	m		matrix to be divided
 *	vp		value to divide by
 *	v3		rounding type parameter
 */
MATRIX *
matquoval(MATRIX *m, VALUE *vp, VALUE *v3)
{
	register VALUE *val, *vres;
	long index;
	MATRIX *res;

	if ((vp->v_type == V_NUM) && qiszero(vp->v_num)) {
		math_error("Division by zero");
		not_reached();
	}
	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, v3, vres++);
	return res;
}


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

	if ((vp->v_type == V_NUM) && qiszero(vp->v_num)) {
		math_error("Division by zero");
		not_reached();
	}
	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, v3, vres++);
	return res;
}


VALUE
mattrace(MATRIX *m)
{
	VALUE *vp;
	VALUE sum;
	VALUE tmp;
	long i, j;

	if (m->m_dim < 2) {
		matsum(m, &sum);
		return sum;
	}
	if (m->m_dim != 2)
		return error_value(E_MATTRACE_2);
	i = (m->m_max[0] - m->m_min[0] + 1);
	j = (m->m_max[1] - m->m_min[1] + 1);
	if (i != j)
		return error_value(E_MATTRACE_3);
	vp = m->m_table;
	copyvalue(vp, &sum);
	j++;
	while (--i > 0) {
		vp += j;
		addvalue(&sum, vp, &tmp);
		freevalue(&sum);
		sum = tmp;
	}
	return sum;
}


/*
 * Transpose a 2-dimensional matrix
 */
MATRIX *
mattrans(MATRIX *m)
{
	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)
		return matcopy(m);
	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(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;
}


/*
 * Round elements of a matrix to specified number of decimal digits
 */
MATRIX *
matround(MATRIX *m, VALUE *v2, VALUE *v3)
{
	VALUE *p, *q;
	long s;
	MATRIX *res;

	s = m->m_size;
	res = matalloc(s);
	*res = *m;
	p = m->m_table;
	q = res->m_table;
	while (s-- > 0)
		roundvalue(p++, v2, v3, q++);
	return res;
}


/*
 * Round elements of a matrix to specified number of binary digits
 */
MATRIX *
matbround(MATRIX *m, VALUE *v2, VALUE *v3)
{
	VALUE *p, *q;
	long s;
	MATRIX *res;

	s = m->m_size;
	res = matalloc(s);
	*res = *m;
	p = m->m_table;
	q = res->m_table;
	while (s-- > 0)
		broundvalue(p++, v2, v3, q++);
	return res;
}

/*
 * Approximate a matrix by approximating elements to be multiples of
 * v2, rounding type determined by v3.
 */
MATRIX *
matappr(MATRIX *m, VALUE *v2, VALUE *v3)
{
	VALUE *p, *q;
	long s;
	MATRIX *res;

	s = m->m_size;
	res = matalloc(s);
	*res = *m;
	p = m->m_table;
	q = res->m_table;
	while (s-- > 0)
		apprvalue(p++, v2, v3, q++);
	return res;
}




/*
 * Produce a matrix with values all of which have been truncated to integers.
 */
MATRIX *
matint(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(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;
}


/*
 * Index a matrix normally by the specified set of index values.
 * Returns the address of the matrix element if it is valid, or generates
 * an error if the index values are out of range.  The create flag is true
 * if the element is to be written, but this is ignored here.
 *
 * given:
 *	mp		matrix to operate on
 *	create		true => create if element does not exist
 *	dim		dimension of the indexing
 *	indices		table of values being indexed by
 */
/*ARGSUSED*/
VALUE *
matindex(MATRIX *mp, bool UNUSED(create), long dim, VALUE *indices)
{
	NUMBER *q;		/* index value */
	VALUE *vp;
	long index;		/* index value as an integer */
	long offset;		/* current offset into array */
	int i;			/* loop counter */

	if (dim < 0) {
		math_error("Negative dimension %ld for matrix", dim);
		not_reached();
	}
	for (;;) {
		if (dim <  mp->m_dim) {
			math_error(
			    "Indexing a %ldd matrix as a %ldd matrix",
			    mp->m_dim, dim);
			not_reached();
		}
		offset = 0;
		for (i = 0; i < mp->m_dim; i++) {
			if (indices->v_type != V_NUM) {
				math_error("Non-numeric index for matrix");
				not_reached();
			}
			q = indices->v_num;
			if (qisfrac(q)) {
				math_error("Non-integral index for matrix");
				not_reached();
			}
			index = qtoi(q);
			if (zge31b(q->num) || (index < mp->m_min[i]) ||
			    (index > mp->m_max[i])) {
				math_error("Index out of bounds for matrix");
				not_reached();
			}
			offset *= (mp->m_max[i] - mp->m_min[i] + 1);
			offset += (index - mp->m_min[i]);
			indices++;
		}
		vp = mp->m_table + offset;
		dim -= mp->m_dim;
		if (dim == 0)
			break;
		if (vp->v_type != V_MAT) {
			math_error("Non-matrix argument for matindex");
			not_reached();
		}
		mp = vp->v_mat;
	}
	return vp;
}


/*
 * Returns the list of indices for a matrix element with specified
 * double-bracket index.
 */
LIST *
matindices(MATRIX *mp, long index)
{
	LIST *lp;
	int j;
	long d;
	VALUE val;

	if (index < 0 || index >= mp->m_size)
		return NULL;

	lp = listalloc();
	val.v_type = V_NUM;
	val.v_subtype = V_NOSUBTYPE;
	j = mp->m_dim;

	while (--j >= 0) {
		d = mp->m_max[j] - mp->m_min[j] + 1;
		val.v_num = itoq(index % d + mp->m_min[j]);
		insertlistfirst(lp, &val);
		qfree(val.v_num);
		index /= d;
	}
	return lp;
}


/*
 * Search a matrix for the specified value, starting with the specified index.
 * Returns 0 and stores index if value found; otherwise returns 1.
 */
int
matsearch(MATRIX *m, VALUE *vp, long i, long j, ZVALUE *index)
{
	register VALUE *val;

	val = &m->m_table[i];
	if (i < 0 || j > m->m_size) {
		math_error("This should not happen in call to matsearch");
		not_reached();
	}
	while (i < j) {
		if (acceptvalue(val++, vp)) {
			utoz(i, index);
			return 0;
		}
		i++;
	}
	return 1;
}


/*
 * Search a matrix backwards for the specified value, starting with the
 * specified index.  Returns 0 and stores index if value found; otherwise
 * returns 1.
 */
int
matrsearch(MATRIX *m, VALUE *vp, long i, long j, ZVALUE *index)
{
	register VALUE *val;

	if (i < 0 || j > m->m_size) {
		math_error("This should not happen in call to matrsearch");
		not_reached();
	}
	val = &m->m_table[--j];
	while (j >= i) {
		if (acceptvalue(val--, vp)) {
			utoz(j, index);
			return 0;
		}
		j--;
	}
	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 w-dimensional 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.
 *
 * given:
 *	m		matrix to be filled
 *	v1		value to fill most of matrix with
 *	v2		value for diagonal entries or null
 */
void
matfill(MATRIX *m, VALUE *v1, VALUE *v2)
{
	register VALUE *val;
	VALUE temp1, temp2;
	long rows, cols;
	long i, j;

	copyvalue(v1, &temp1);

	val = m->m_table;
	if (m->m_dim != 2 || v2 == NULL) {
		for (i = m->m_size; i > 0; i--) {
			freevalue(val);
			copyvalue(&temp1, val++);
		}
		freevalue(&temp1);
		return;
	}

	copyvalue(v2, &temp2);
	rows = m->m_max[0] - m->m_min[0] + 1;
	cols = m->m_max[1] - m->m_min[1] + 1;

	for (i = 0; i < rows; i++) {
		for (j = 0; j < cols; j++) {
			freevalue(val);
			if (i == j)
				copyvalue(&temp2, val++);
			else
				copyvalue(&temp1, val++);
		}
	}
	freevalue(&temp1);
	freevalue(&temp2);
}



/*
 * Set a copy of a square matrix to the identity matrix.
 */
S_FUNC MATRIX *
matident(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) {
		math_error(
		    "Matrix dimension must be two for setting to identity");
		not_reached();
	}
	if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])) {
		math_error("Matrix must be square for setting to identity");
		not_reached();
	}
	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_subtype = V_NOSUBTYPE;
			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(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 *vres;		/* current value in result for dim < 2 */
	VALUE mulval;		/* value to multiply rows by */
	VALUE tmpval;		/* temporary value */

	if (m->m_dim < 2) {
		res = matalloc(m->m_size);
		*res = *m;
		val = m->m_table;
		vres = res->m_table;
		for (cur = m->m_size; cur > 0; cur--)
			invertvalue(val++, vres++);
		return res;
	}
	if (m->m_dim != 2) {
		math_error("Matrix dimension exceeds two for inverse");
		not_reached();
	}
	if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])) {
		math_error("Inverting non-square matrix");
		not_reached();
	}
	/*
	 * 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->v_subtype = V_NOSUBTYPE;
			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);
				math_error("Matrix is not invertible");
				not_reached();
			}
			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++) {
			if ((row == cur) || (testvalue(val) == 0)) {
				if (row+1 < rows)
					val += rows;
				continue;
			}
			mulvalue(val, &mulval, &tmpval);
			matsubrow(m, row, cur, &tmpval);
			matsubrow(res, row, cur, &tmpval);
			freevalue(&tmpval);
			if (row+1 < rows)
				val += rows;
		}
		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);
		}
		if (row+1 < rows)
			val += (rows + 1);
	}
	matfree(m);
	return res;
}


/*
 * Calculate the determinant of a square matrix.
 * This uses the fraction-free Gauss-Bareiss algorithm.
 */
VALUE
matdet(MATRIX *m)
{
	long n;			/* original matrix is n x n */
	long k;			/* working sub-matrix is k x k */
	long i, j;
	VALUE *pivot, *div, *val;
	VALUE *vp, *vv;
	VALUE tmp1, tmp2, tmp3;
	bool neg;		/* whether to negate determinant */

	if (m->m_dim < 2) {
		vp = m->m_table;
		i = m->m_size;
		copyvalue(vp, &tmp1);

		while (--i > 0) {
			mulvalue(&tmp1, ++vp, &tmp2);
			freevalue(&tmp1);
			tmp1 = tmp2;
		}
		return tmp1;
	}

	if (m->m_dim != 2)
		return error_value(E_DET_2);
	if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1]))
		return error_value(E_DET_3);

	/*
	 * 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 = false;
	m = matcopy(m);
	n = (m->m_max[0] - m->m_min[0] + 1);
	pivot = div = m->m_table;
	for (k = n; k > 0; k--) {
		/*
		 * Find the first nonzero value in the rest of the column
		 * downwards from pivot.  If there is no such value, then
		 * the determinant is zero.  If the first nonzero entry is not
		 * the pivot, then swap rows in the k * k sub-matrix, and
		 * remember that the determinant changes sign.
		 */
		val = pivot;
		i = k;
		while (!testvalue(val)) {
			if (--i <= 0) {
				tmp1.v_type = V_NUM;
				tmp1.v_subtype = V_NOSUBTYPE;
				tmp1.v_num = qlink(&_qzero_);
				matfree(m);
				return tmp1;
			}
			val += n;
		}
		if (i < k) {
			vp = pivot;
			vv = val;
			j = k;
			while (j-- > 0) {
				tmp1 = *vp;
				*vp++ = *vv;
				*vv++ = tmp1;
			}
			neg = !neg;
		}
		/*
		 * Now for every val below the pivot, for each entry to
		 * the right of val, calculate the 2 x 2 determinant
		 * with corners at the pivot and the entry.  If
		 * k < n, divide by div (the previous pivot value).
		 */
		val = pivot;
		i = k;
		while (--i > 0) {
			val += n;
			vp = pivot;
			vv = val;
			j = k;
			while (--j > 0) {
				mulvalue(pivot, ++vv, &tmp1);
				mulvalue(val, ++vp, &tmp2);
				subvalue(&tmp1, &tmp2, &tmp3);
				freevalue(&tmp1);
				freevalue(&tmp2);
				freevalue(vv);
				if (k < n) {
					divvalue(&tmp3, div, vv);
					freevalue(&tmp3);
				}
				else
					*vv = tmp3;
			}
		}
		div = pivot;
		if (k > 0)
			pivot += n + 1;
	}
	if (neg)
		negvalue(div, &tmp1);
	else
		copyvalue(div, &tmp1);
	matfree(m);
	return tmp1;
}


/*
 * Local utility routine to swap two rows of a square matrix.
 * No checks are made to verify the legality of the arguments.
 */
S_FUNC void
matswaprow(MATRIX *m, long r1, long 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.
 */
S_FUNC void
matsubrow(MATRIX *m, long oprow, long 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++;
	}
}


/*
 * Local utility routine to multiply a row by a specified number.
 * No checks are made to verify the legality of the arguments.
 */
S_FUNC void
matmulrow(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(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) {
		copyvalue(v1, v2);
		v1++;
		v2++;
	}
	return res;
}


/*
 * Make a matrix the same size as another and filled with a fixed value.
 *
 * given:
 *	m		matrix to initialize
 *	v1		value to fill most of matrix with
 *	v2		value for diagonal entries (or NULL)
 */
MATRIX *
matinit(MATRIX *m, VALUE *v1, VALUE *v2)
{
	MATRIX *res;
	register VALUE *v;
	register long i;
	long row;
	long rows;

	/*
	 * clone matrix size
	 */
	res = matalloc(m->m_size);
	*res = *m;

	/*
	 * firewall
	 */
	if (v2 && ((res->m_dim != 2) ||
		((res->m_max[0] - res->m_min[0]) !=
		 (res->m_max[1] - res->m_min[1])))) {
			math_error("Filling diagonals of non-square matrix");
			not_reached();
	}

	/*
	 * fill the bulk of the matrix
	 */
	v = res->m_table;
	if (v2 == NULL) {
		i = m->m_size;
		while (i-- > 0) {
			copyvalue(v1, v++);
		}
		return res;
	}

	/*
	 * fill the diagonal of a square matrix if requested
	 */
	rows = res->m_max[0] - res->m_min[0] + 1;
	v = res->m_table;
	for (row = 0; row < rows; row++) {
		copyvalue(v2, v+row);
		v += rows;
	}
	return res;
}


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

	m = (MATRIX *) malloc(matsize(size));
	if (m == NULL) {
		math_error("Cannot get memory to allocate matrix of size %ld",
			   size);
		not_reached();
	}
	m->m_size = size;
	for (i = size, vp = m->m_table; i > 0; i--, vp++)
		vp->v_subtype = V_NOSUBTYPE;
	return m;
}


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

	vp = m->m_table;
	i = m->m_size;
	while (i-- > 0)
		freevalue(vp++);
	free(m);
}


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

	vp = m->m_table;
	i = m->m_size;
	while (i-- > 0) {
		if (testvalue(vp++))
			return true;
	}
	return false;
}

/*
 * Sum the elements in a matrix.
 */
void
matsum(MATRIX *m, VALUE *vres)
{
	VALUE *vp;
	VALUE tmp;		/* first sum value */
	VALUE sum;		/* final sum value */
	long i;

	vp = m->m_table;
	i = m->m_size;
	copyvalue(vp, &sum);

	while (--i > 0) {
		addvalue(&sum, ++vp, &tmp);
		freevalue(&sum);
		sum = tmp;
	}
	*vres = sum;
}


/*
 * 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(MATRIX *m1, MATRIX *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;
}


void
matreverse(MATRIX *m)
{
	VALUE *p, *q;
	VALUE tmp;

	p = m->m_table;
	q = m->m_table + m->m_size - 1;
	while (q > p) {
		tmp = *p;
		*p++ = *q;
		*q-- = tmp;
	}
}


void
matsort(MATRIX *m)
{
	VALUE *a, *b, *next, *end;
	VALUE *buf, *p;
	VALUE *S[LONG_BITS];
	long len[LONG_BITS];
	long i, j, k;

	buf = (VALUE *) malloc(m->m_size * sizeof(VALUE));
	if (buf == NULL) {
		math_error("Not enough memory for matsort");
		not_reached();
	}
	next = m->m_table;
	end = next + m->m_size;
	for (k = 0; next && k < LONG_BITS; k++) {
		S[k] = next++;			/* S[k] is start of a run */
		len[k] = 1;
		if (next == end)
			next = NULL;
		while (k > 0 && (!next || len[k] >= len[k - 1])) {/* merging */
			j = len[k];
			b = S[k--];
			i = len[k];
			a = S[k];
			len[k] += j;
			p = buf;
			if (precvalue(b, a)) {
				do {
					*p++ = *b++;
					j--;
				} while (j > 0 && precvalue(b,a));
				if (j == 0) {
					memcpy(p, a, i * sizeof(VALUE));
					memcpy(S[k], buf,
					       len[k] * sizeof(VALUE));
					continue;
				}
			}

			do {
				do {
					*p++ = *a++;
					i--;
				} while (i > 0 && !precvalue(b,a));
				if (i == 0) {
					break;
				}
				do {
					*p++ = *b++;
					j--;
				} while (j > 0 && precvalue(b,a));
			} while (j != 0);

			if (i == 0) {
				memcpy(S[k], buf, (p - buf) * sizeof(VALUE));
			} else if (j == 0) {
				memcpy(p, a, i * sizeof(VALUE));
				memcpy(S[k], buf, len[k] * sizeof(VALUE));
			}
		}
	}
	free(buf);
	if (k >= LONG_BITS) {
		/* this should never happen */
		math_error("impossible k overflow in matsort!");
		not_reached();
	}
}

void
matrandperm(MATRIX *m)
{
	VALUE *vp;
	long s, i;
	VALUE val;

	s = m->m_size;
	for (vp = m->m_table; s > 1; vp++, s--) {
		i = irand(s);
		if (i > 0) {
			val = *vp;
			*vp = vp[i];
			vp[i] = val;
		}
	}
}


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

	val = m->m_table;
	if (m->m_dim == 0) {
		return (val->v_type == V_NUM && qisone(val->v_num));
	}
	if (m->m_dim == 1) {
		for (row = m->m_min[0]; row <= m->m_max[0]; row++, val++) {
			if (val->v_type != V_NUM || !qisone(val->v_num))
				return false;
		}
		return true;
	}
	if ((m->m_dim != 2) ||
		((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])))
			return false;
	for (row = m->m_min[0]; row <= m->m_max[0]; row++) {
		/*
		 * We could use col = m->m_min[1]; col < m->m_max[1]
		 * but if m->m_min[0] != m->m_min[1] this won't work.
		 * We know that we have a square 2-dimensional matrix
		 * so we will pretend that m->m_min[0] == m->m_min[1].
		 */
		for (col = m->m_min[0]; col <= m->m_max[0]; 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;
}


/*
 * 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(MATRIX *m, long max_print)
{
	VALUE *vp;
	long fullsize, count, index;
	long dim, i, j;
	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 [");
	if (dim) {
		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 = ",";
		}
	} else {
		math_str("mat [");
	}
	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 unambiguous format.
	 */
	math_str(":\n");
	vp = m->m_table;
	for (index = 0; index < max_print; index++) {
		msg = "	 [";
		j = index;
		if (dim) {
			for (i = 0; i < dim; i++) {
				math_fmt("%s%ld", msg,
					 m->m_min[i] + (j / sizes[i]));
				j %= sizes[i];
				msg = ",";
			}
		} else {
			math_str(msg);
		}
		math_str("] = ");
		printvalue(vp++, PRINT_SHORT | PRINT_UNAMBIG);
		math_str("\n");
	}
	if (max_print < fullsize)
		math_str("  ...\n");
}