matfunc.c

Calc Software Package for Number Calc

	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_DET2);
	if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1]))
		return error_value(E_DET3);

	/*
	 * 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 submatrix, 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.
 */
static 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.
 */
static 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.
 */
static 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");
			/*NOTREACHED*/
	}

	/*
	 * 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 diaginal 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 %d",
			   size);
		/*NOTREACHED*/
	}
	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");
		/*NOTREACHED*/
	}
	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!");
		/*NOTREACHED*/
	}
}

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 unambigous 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");
}