/*	$Id: matrici.c,v 1.3 1997/05/01 17:34:49 sandro Exp $	*/

/*
 * Copyright (c) 1997
 *	Sandro Sigala, Brescia, Italy.  All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "matrici.h"
#include "umatrix.h"

int stile = STILE_CLASSICO;

double (valore)(Matrice m, int riga, int colonna)
{
	return m->valori[(riga-1) * m->colonne + (colonna-1)];
}

Matrice matr_crea(char *nome, int righe, int colonne)
{
	Matrice m;
	assert(righe > 0 && colonne > 0);
	m = (Matrice)malloc(sizeof *m);
	m->righe = righe;
	m->colonne = colonne;
	m->nome = strdup(nome);
	m->valori = (double *)calloc(righe * colonne, sizeof(double));
	return m;
}

void matr_elimina(Matrice m)
{
	assert(m != NULL);
	free(m->nome);
	free(m->valori);
	free(m);
}

void matr_stampa(Matrice m)
{
	int r, i, c, sp;

	assert(m != NULL);

	sp = strlen(m->nome) + 3;

	for (r = 1; r <= m->righe; r++) {
		if (sp) {
			if (r-1 == m->righe / 2)
				printf("%s = ", m->nome);
			else
				for (i = 0; i < sp; i++)
					putchar(' ');
		}
		if (stile == STILE_ALTERNATIVO)
			printf("|| ");
		else {
			if (m->righe > 1) {
				if (r == 1) printf(" / ");
				else if (r == m->righe) printf(" \\ ");
				else printf("|  ");
			} else printf("| ");
		}

		for (c = 1; c <= m->colonne; c++)
			printf("%3.0f ", valore(m, r, c));
		if (stile == STILE_ALTERNATIVO)
			printf("||\n");
		else {
			if (m->righe > 1) {
				if (r == 1) printf(" \\\n");
				else if (r == m->righe) printf(" /\n");
				else printf("  |\n");
			} else printf("|\n");
		}
	}
}

Matrice matr_somma(char *nome, Matrice m1, Matrice m2)
{
	Matrice mr;
	int r, c;

	assert(m1 != NULL && m2 != NULL);

	if (m1->righe != m2->righe || m1->colonne != m2->colonne) {
		fprintf(stderr, getstr("Errore: somma di matrici di dimensioni differenti\n"));
		return NULL;
	}

	mr = matr_crea(nome, m1->righe, m1->colonne);

	for (r = 1; r <= m1->righe; r++)
	for (c = 1; c <= m1->colonne; c++)
		valore(mr, r, c) = valore(m1, r, c) + valore(m2, r, c);

	return mr;
}

Matrice matr_prodotto_per_numero(char *nome, Matrice m, double n)
{
	Matrice mr;
	int r, c;

	assert(m != NULL);

	mr = matr_crea(nome, m->righe, m->colonne);

	for (r = 1; r <= m->righe; r++)
	for (c = 1; c <= m->colonne; c++)
		valore(mr, r, c) = valore(m, r, c) * n;

	return mr;
}

Matrice matr_prodotto(char *nome, Matrice m1, Matrice m2)
{
	Matrice mr;
	int r, c, i;

	assert(m1 != NULL && m2 != NULL);

	if (m1->colonne != m2->righe) {
		fprintf(stderr, getstr("Errore: prodotto di matrici di dimensioni errate\n"));
		return NULL;
	}

	mr = matr_crea(nome, m1->righe, m2->colonne);

	for (r = 1; r <= m1->righe; r++)
	for (c = 1; c <= m2->colonne; c++)
	for (i = 1; i <= m1->colonne; i++)
		valore(mr, r, c) += valore(m1, r, i) * valore(m2, i, c);

	return mr;
}

Matrice matr_trasposta(char *nome, Matrice m)
{
	Matrice mr;
	int r, c;

	assert(m != NULL);

	mr = matr_crea(nome, m->colonne, m->righe);

	for (r = 1; r <= m->righe; r++)
	for (c = 1; c <= m->colonne; c++)
		valore(mr, c, r) = valore(m, r, c);

	return mr;
}

void matr_assegna_valori_casuali(Matrice m, int min, int max)
{
	int r, c;

	assert(m != NULL);

	for (r = 1; r <= m->righe; r++)
	for (c = 1; c <= m->colonne; c++)
		valore(m, r, c) = rand() % (max + 1 - min) + min;
}

Matrice matr_crea_valori_casuali(char *nome, int righe, int colonne, int min, int max)
{
	Matrice m;

	m = matr_crea(nome, righe, colonne);
	matr_assegna_valori_casuali(m, min, max);

	return m;
}

Matrice matr_crea_matrice_unita(char *nome, int ordine)
{
	Matrice m;
	int i;

	m = matr_crea(nome, ordine, ordine);
	for (i = 1; i <= ordine; i++)
		valore(m, i, i) = 1;

	return m;
}

static double determinante(Matrice m)
{
	double d;
	int r, c, i, k;
	Matrice mr;

	/*
	 * Versioni inline non ricorsive del calcolo del determinante di
	 * matrice di ordine 1, 2 e 3.
	 */
	switch (m->righe) {
	case 1:
		return valore(m, 1, 1);
	case 2:
		return valore(m, 1, 1) * valore(m, 2, 2) -
			valore(m, 1, 2) * valore(m, 2, 1);
	case 3:
		return valore(m, 1, 1) * valore(m, 2, 2) * valore(m, 3, 3) -
			valore(m, 1, 1) * valore(m, 2, 3) * valore(m, 3, 2) -
			valore(m, 1, 2) * valore(m, 2, 1) * valore(m, 3, 3) +
			valore(m, 1, 3) * valore(m, 2, 1) * valore(m, 3, 2) +
			valore(m, 1, 2) * valore(m, 2, 3) * valore(m, 3, 1) -
			valore(m, 1, 3) * valore(m, 2, 2) * valore(m, 3, 1);
	}

	/* Matrice di ordine > 3. */

	for (d = 0, k = 1; k <= m->righe; k++) {
		if (valore(m, k, 1) == 0)
			continue;
		mr = matr_crea("", m->righe-1, m->righe-1);
		for (r = i = 1; r <= m->righe; r++) {
			if (r == k)
				continue;
			for (c = 1; c + 1 <= m->colonne; c++)
				valore(mr, i, c) = valore(m, r, c+1);
			i++;
		}
		d += determinante(mr) * valore(m, k, 1) * pow(-1, k + 1);
		matr_elimina(mr);
	}

	return d;
}

/*
 * matr_determinante(): funzione.
 * La funzione restitusce il determinante della matrice.
 */
double matr_determinante(Matrice m)
{
	if (m->righe != m->colonne) {
		fprintf(stderr, getstr("Errore: calcolo determinante di matrice non quadrata\n"));
		return 0;
	}

	return determinante(m);
}

/*
 * Predicati.
 */

/*
 * matr_quadrata_p(): predicato.
 * La funzione e` vera se la matrice e` quadrata.
 */
int (matr_quadrata_p)(Matrice m)
{
	assert(m != NULL);

	return m->righe == m->colonne;
}

/*
 * matr_rettangolare_p(): predicato.
 * La funzione e` vera se la matrice e` rettangolare.
 */
int (matr_rettangolare_p)(Matrice m)
{
	assert(m != NULL);

	return m->righe != m->colonne;
}

/*
 * matr_unita_p(): predicato.
 * La funzione e` vera se la matrice e` una matrice unita`.
 *
 * Una matrice quadrata A di ordine n si dice matrice unita` se i suoi
 * elementi della diagonale principale sono eguali a 1 e tutti gli altri
 * elementi eguali a zero, cioe` se:
 *   per i,k=1..n --> a[i,k] = 0 con i != k e
 *   per i=1..n --> a[i,k] = 1
 */
int matr_unita_p(Matrice m)
{
	int r, c;

	assert(m != NULL);

	if (!matr_quadrata_p(m))
		return 0;

	for (r = 1; r <= m->righe; r++)
	for (c = 1; c <= m->colonne; c++)
		if (r == c) {
			/* Diagonale principale. */
			if (valore(m, r, c) != 1)
				return 0;
		} else if (valore(m, r, c) != 0)
				return 0;

	return 1;
}

/*
 * matr_simmetrica_p(): predicato.
 * La funzione e` vera se la matrice e` simmetrica.
 *
 * Una matrice quadrata A di ordine n si dice simmetrica se sono eguali le
 * coppie di elementi simmetrici rispetto alla diagonale principale, cioe` se:
 *   a[i,k] = a[k,i] con i,k=1..n e i != k
 *
 * In altre parole, una matrice e` simmetrica solo e solo se coincide con la
 * sua trasposta.
 */
int matr_simmetrica_p(Matrice m)
{
	int r, c;

	assert(m != NULL);

	if (!matr_quadrata_p(m))
		return 0;

	for (r = 1; r <= m->righe; r++)
	for (c = 1; c <= m->colonne; c++)
		if (valore(m, r, c) != valore(m, c, r))
			return 0;

	return 1;
}

/*
 * matr_stesso_tipo_p(): predicato.
 * La funzione e` vera se le due matrici sono dello stesso tipo, ovvero
 * hanno la stessa dimensione.
 */
int (matr_stesso_tipo_p)(Matrice m1, Matrice m2)
{
	assert(m1 != NULL && m2 != NULL);

	return m1->righe == m2->righe && m1->colonne == m2->colonne;
}

/*
 * matr_eguali_p(): predicato.
 * La funzione e` vera se le due matrici sono eguali.
 */
int matr_eguali_p(Matrice m1, Matrice m2)
{
	assert(m1 != NULL && m2 != NULL);

	if (!matr_stesso_tipo_p(m1, m2))
		return 0;

#if 1
	/* Versione classica. */
	{
		int r, c;
		for (r = 1; r <= m1->righe; r++)
		for (c = 1; c <= m1->colonne; c++)
			if (valore(m1, r, c) != valore(m2, r, c))
				return 0;
		return 1;
	}
#else
	/* Versione ottimizata.  In alcuni casi non puo` funzionare.  */
	return memcmp(m1->valori, m2->valori,
		      m1->righe * m1->colonne * sizeof(double)) == 0;
#endif
}