/*
 * Copyright 1995,96 Thierry Bousch
 * Licensed under the Gnu Public License, Version 2
 *
 * $Id: Matrix.c,v 2.9 1996/09/14 09:40:04 bousch Exp $
 *
 * Rectangular matrices, on an arbitrary ring. Matrices are represented
 * as arrays of "Lines" (one-dimensional arrays)
 */

#include <limits.h>
#include <stdlib.h>
#include "saml.h"
#include "saml-errno.h"
#include "saml-util.h"
#include "mnode.h"
#include "builtin.h"

static gr_string* array_stringify (std_mnode*);
static s_mnode* array_add (std_mnode*, std_mnode*);
static s_mnode* array_sub (std_mnode*, std_mnode*);
static int array_notzero (std_mnode*);
static int array_differ (std_mnode*, std_mnode*);
static s_mnode* array_zero (std_mnode*);
static s_mnode* array_negate (std_mnode*);
static s_mnode* array2array (std_mnode*, std_mnode*);

static s_mnode* matrix_mul (std_mnode*, std_mnode*);
static s_mnode* matrix_one (std_mnode*);
static s_mnode* matrix_invert (std_mnode*);
static s_mnode* matrix2matrix (std_mnode*, std_mnode*);
s_mnode* matrix_determinant (std_mnode*);

static unsafe_s_mtype MathType_Line = {
	"Line",
	mstd_free, NULL, array_stringify,
	NULL, NULL,
	array_add, array_sub, NULL, NULL, NULL,
	array_notzero, NULL, NULL, array_differ, NULL,
	array_zero, array_negate, NULL, NULL, NULL
};

static unsafe_s_mtype MathType_Matrix = {
	"Matrix",
	mstd_free, NULL, array_stringify,
	NULL, NULL,
	array_add, array_sub, matrix_mul, NULL, NULL,
	array_notzero, NULL, NULL, array_differ, NULL,
	array_zero, array_negate, matrix_one, matrix_invert, NULL
};

void init_MathType_Matrix (void)
{
	register_mtype(ST_LINE, &MathType_Line);
	register_mtype(ST_MATRIX, &MathType_Matrix);
	register_CV_routine(ST_LINE, ST_LINE, array2array);
	register_CV_routine(ST_MATRIX, ST_MATRIX, matrix2matrix);
}

static gr_string* array_stringify (std_mnode* ar)
{
	gr_string *grs = new_gr_string(0);
	return grs_append(grs, "array", 5);
}

static s_mnode* array2array (std_mnode* ar, std_mnode* model)
{
	std_mnode *mn;
	s_mnode *cmodel;
	int i, len;

	if (!model)
		return copy_mnode((mn_ptr)ar);
	len = model->length;
	if (ar->length != len)
		return mnode_error(SE_ICAST, "array2array");

	mn = mstd_alloc(ST_LINE, len);
	cmodel = model->x[0];
	for (i = 0; i < len; i++)
		mn->x[i] = mnode_promote(ar->x[i], cmodel);
	return (mn_ptr)mn;
}

static s_mnode* matrix2matrix (std_mnode* ar, std_mnode* model)
{
	return mnode_error(SE_NOTRDY, "matrix2matrix");
}


static s_mnode* array_add (std_mnode *ar1, std_mnode *ar2)
{
	int i, len;
	std_mnode *mn;

	if ((len = ar1->length) != (ar2->length))
		return mnode_error(SE_SCONFL, "array_add");
	mn = mstd_alloc(ar1->hdr.type, len);
	for (i = 0; i < len; i++)
		mn->x[i] = mnode_add(ar1->x[i], ar2->x[i]);
	return (mn_ptr)mn;
}

static s_mnode* array_sub (std_mnode* ar1, std_mnode* ar2)
{
	int i, len;
	std_mnode *mn;

	if ((len = ar1->length) != (ar2->length))
		return mnode_error(SE_SCONFL, "array_sub");
	mn = mstd_alloc(ar1->hdr.type, len);
	for (i = 0; i < len; i++)
		mn->x[i] = mnode_sub(ar1->x[i], ar2->x[i]);
	return (mn_ptr)mn;
}

static int array_notzero (std_mnode* ar)
{
	int i, len = ar->length;

	for (i = 0; i < len; i++)
	    if (mnode_notzero(ar->x[i]))
	    	return 1;
	return 0;
}

static int array_differ (std_mnode* ar1, std_mnode* ar2)
{
	int i, len1=ar1->length, len2=ar2->length;

	if (len1 != len2)
	    return 1;
	for (i = 0; i < len1; i++)
	    if (mnode_differ(ar1->x[i],ar2->x[i]))
		return 1;
	return 0;
}

static s_mnode* array_zero (std_mnode* ar1)
{
	int i, len = ar1->length;
	s_mnode *zero;
	std_mnode *mn;

	if (len == 0)
		return copy_mnode((mn_ptr)ar1);
	zero = mnode_zero(ar1->x[0]);
	mn = mstd_alloc(ar1->hdr.type, len);
	for (i = 0; i < len; i++)
		mn->x[i] = zero;
	zero->refs += len - 1;
	return (mn_ptr)mn;
}

static s_mnode* array_negate (std_mnode* ar1)
{
	int i, len = ar1->length;
	std_mnode *mn;

	mn = mstd_alloc(ar1->hdr.type, len);
	for (i = 0; i < len; i++)
		mn->x[i] = mnode_negate(ar1->x[i]);
	return (mn_ptr)mn;
}

static s_mnode* matrix_mul (std_mnode* ar1, std_mnode* ar2)
{
	return mnode_error(SE_NOTRDY, "matrix_mul");
}

static s_mnode* matrix_invert (std_mnode* ar1)
{
	return mnode_error(SE_NOTRDY, "matrix_invert");
}

static s_mnode* matrix_one (std_mnode* model)
{
	int i, j, len = model->length;
	s_mnode *eltmodel, *zero, *one;
	std_mnode *mn, *mnl;

	if (len == 0) {
		/* No rows? */
		return copy_mnode((mn_ptr)model);
	}
	mnl = (smn_ptr)(model->x[0]);
	if (mnl->length != len) {
		/* Not a square matrix */
		return mnode_error(SE_SCONFL, "matrix_one");
	}
	eltmodel = mnl->x[0];
	zero = mnode_zero(eltmodel);
	one  = mnode_one(eltmodel);
	mn = mstd_alloc(ST_MATRIX, len);
	for (i = 0; i < len; i++) {
		mnl = mstd_alloc(ST_LINE, len);
		mn->x[i] = (mn_ptr)mnl;
		for (j = 0; j < len; j++)
			mnl->x[j] = zero;
		mnl->x[i] = one;
	}
	/* Adjust the reference counters */
	zero->refs += len * (len - 1) - 1;
	one ->refs += len - 1;
	return (mn_ptr)mn;
}

s_mnode* matrix_determinant (std_mnode* mat)
{
	int i, j, k, chsign = 0, len = mat->length;
	std_mnode *line;
	s_mnode ***m, *npiv, *det;

	if (!len) {
		/* Huh? The determinant of a void matrix? */
		return mnode_error(SE_WR_SIZE, "matrix_determinant");
	}
	line = (smn_ptr)(mat->x[0]);
	if (line->length != len) {
		/* Not a square matrix */
		return mnode_error(SE_SCONFL, "matrix_one");
	}
	if (len == 1) {
		/* A one-by-one determinant: that's easy! */
		return copy_mnode(line->x[0]);
	}
	/*
	 * For the general case, we use Bareiss' algorithm. We try to
	 * choose a pivot of minimal complexity, using the ->length field
	 * as an approximate measure of complexity -- this will work well
	 * for univariate and multivariate polynomials.
	 *
	 *   Bareiss, E.H., "Sylvester's identity and multistep
	 *   integer-preserving gaussian elimination",
	 *   Math. Comp. 22, pp 565-578 (1968)
	 *
	 * Temporary pointers-to-mnodes are put on the stack.
	 */
	m = (s_mnode ***)alloca(len * sizeof(s_mnode**));
	for (i = 0; i < len; i++) {
		m[i] = (s_mnode **)alloca(len * sizeof(s_mnode*));
		line = (smn_ptr)(mat->x[i]);
		for (j = 0; j < len; j++)
			m[i][j] = copy_mnode(line->x[j]);
	}
	/* Now the main loop */
	for (k = 0; k < len-1; k++) {
		int best_len = INT_MAX, cur_len, best_idx = -1;
		extern int apoly_length(s_mnode*);

		/* Find an appropriate pivot */
		for (i = k; i < len; i++) {
			npiv = m[i][k];
			if (!mnode_notzero(npiv))
				continue;

			if (npiv->type == ST_POLY)
				cur_len = ((smn_ptr)npiv)->length;
			else if (npiv->type == ST_APOLY)
				cur_len = apoly_length(npiv);
			else
				cur_len = 0;	/* unknown complexity */

			/* Too complicated? */
			if (cur_len >= best_len)
				continue;
			/* Found it */
			best_idx = i;
			best_len = cur_len;
			if (best_len <= 1)
				break;
		}
		if (best_idx < 0) {
			/*
			 * The k-th column is entirely zero, so we can free
			 * everything and report that the determinant
			 * is zero.
			 */
			for (i = 0; i < len; i++)
			    for (j = 0; j < len; j++)
				unlink_mnode(m[i][j]);
			line = (smn_ptr)(mat->x[0]);
			return mnode_zero(line->x[0]);
		}
		if ((i = best_idx) != k) {
			/* Swap lines i and k */
			s_mnode **tmp;
			tmp = m[i]; m[i] = m[k]; m[k] = tmp;
			++chsign;
		}
		npiv = m[k][k];
		for (i = k+1; i < len; i++)
		    for (j = k+1; j < len; j++) {
		    	s_mnode *t1, *t2, *t3;
		    	t1 = mnode_mul(npiv, m[i][j]);
		    	unlink_mnode(m[i][j]);
		    	t2 = mnode_mul(m[i][k], m[k][j]);
		    	t3 = mnode_sub(t1, t2);
		    	unlink_mnode(t1); unlink_mnode(t2);
		    	if (k == 0)
		    		m[i][j] = t3;
		    	else {
		    		m[i][j] = mnode_div(t3, m[k-1][k-1]);
		    		unlink_mnode(t3);
		    	}
		    }
	}
	/* The result is now in m[len-1][len-1], up to the sign */
	if (chsign & 1)
		det = mnode_negate(m[len-1][len-1]);
	else
		det = copy_mnode(m[len-1][len-1]);

	/* Free our temporary matrix and return the result */
	for (i = 0; i < len; i++)
	    for (j = 0; j < len; j++)
	    	unlink_mnode(m[i][j]);
	return det;
}