matfunc.c

Calc Software Package for Number Calc

	long index;
	MATRIX *res;

	if ((vp->v_type == V_NUM) && qiszero(vp->v_num)) {
		math_error("Division by zero");
		/*NOTREACHED*/
	}
	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");
		/*NOTREACHED*/
	}
	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_MATTRACE2);
	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_MATTRACE3);
	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 elemenbs 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);
		/*NOTREACHED*/
	}
	for (;;) {
		if (dim <  mp->m_dim) {
			math_error(
			    "Indexing a %ldd matrix as a %ldd matrix",
			    mp->m_dim, dim);
			/*NOTREACHED*/
		}
		offset = 0;
		for (i = 0; i < mp->m_dim; i++) {
			if (indices->v_type != V_NUM) {
				math_error("Non-numeric index for matrix");
				/*NOTREACHED*/
			}
			q = indices->v_num;
			if (qisfrac(q)) {
				math_error("Non-integral index for matrix");
				/*NOTREACHED*/
			}
			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");
				/*NOTREACHED*/
			}
			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");
			/*NOTREACHED*/
		}
		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;
	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");
		/*NOTREACHED*/
	}
	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");
		/*NOTREACHED*/
	}
	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.
 */
static 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");
		/*NOTREACHED*/
	}
	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");
		/*NOTREACHED*/
	}
	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(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");
		/*NOTREACHED*/
	}
	if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])) {
		math_error("Inverting non-square matrix");
		/*NOTREACHED*/
	}
	/*
	 * 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);
				math_error("Matrix is not invertible");
				/*NOTREACHED*/
			}
			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 submatrix is k x k */
	long i, j;