matfunc.c

Calc Software Package for Number Calc

/*
 * matfunc - extended precision rational arithmetic matrix functions
 *
 * Copyright (C) 1999-2004  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.
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 *
 * @(#) $Revision: 29.7 $
 * @(#) $Id: matfunc.c,v 29.7 2006/06/02 10:24:09 chongo Exp $
 * @(#) $Source: /usr/local/src/cmd/calc/RCS/matfunc.c,v $
 *
 * 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 "value.h"
#include "zrand.h"
#include "calcerr.h"

#include "have_unused.h"

extern long irand(long s);

static void matswaprow(MATRIX *m, long r1, long r2);
static void matsubrow(MATRIX *m, long oprow, long baserow, VALUE *mulval);
static void matmulrow(MATRIX *m, long row, VALUE *mulval);
static 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");
		/*NOTREACHED*/
	}
	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");
			/*NOTREACHED*/
		}
		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");
		/*NOTREACHED*/
	}
	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");
			/*NOTREACHED*/
		}
		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");
			/*NOTREACHED*/
		}
		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");
			/*NOTREACHED*/
		}
		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");
			/*NOTREACHED*/
		}
		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");
		/*NOTREACHED*/
	}
	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");
		/*NOTREACHED*/
	}
	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");
		/*NOTREACHED*/
	}
	if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])) {
		math_error("Squaring non-square matrix");
		/*NOTREACHED*/
	}
	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");
		/*NOTREACHED*/
	}
	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");
		/*NOTREACHED*/
	}
	if (qisfrac(q)) {
		math_error("Raising matrix to non-integral power");
		/*NOTREACHED*/
	}
	if (zge31b(q->num)) {
		math_error("Raising matrix to very large power");
		/*NOTREACHED*/
	}
	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;