
#include <stdio.h>
#include "matrix.h"

/* This file defines functions to calculate determinant of a matrix:
 * 	A function that calculates determinant of a permutation matrix
 * 	A function that calculates determinant of an arbitrary matrix
 */



/* Computes the determinant of permutation matrix P.
 * It is +1, or -1 since every row and column of P has exactly one non-zero entry,
 * and that is equal to 1.
 */
float compute_det_permutation(Matrix P)
{
	Matrix Q; // stores a submatrix of P
	int i; // P[0][i] = 1
	float det; // determinant value

	if (P.rows == 1) // P = [1]
		return 1;

	i = find_nonzero_column(P);

	// Drop first row and ith column of P and copy into Q
	Q = allocate_matrix(P.rows-1, P.cols-1);

	for (int j = 0; j < P.rows-1; j++)
		for (int k = 0; k < P.cols-1; k++)
			if (k < i)
				Q.elements[j][k] = P.elements[j+1][k];
			else
				Q.elements[j][k] = P.elements[j+1][k+1];
	
	// determinant of P equals determinant of Q times (-1)^i
	if (i % 2 == 1) // i is odd
		det =  -compute_det_permutation(Q);
	else
		det =  compute_det_permutation(Q);

	// free space and return
	free_matrix(Q);
	return det;
}


/* Computes the determinant of square matrix A */
float determinant(Matrix A)
{
	int i;
	Matrix B; // stores a submatrix of A
	Matrix T; // copy of A
	float det; // determinant value
	
	if (A.rows == 1) // 1 x 1 matrix
		return A.elements[0][0];

	T = allocate_matrix(A.rows, A.cols);
	copy_matrix(T, 0, A, 0, A.rows, A.cols); // copy A to T; now work on T

	if (T.elements[0][0] == 0) {
		i = find_nonzero_row(T);
		if (i >= T.rows) { // no non-zero row
			free_matrix(T);
			return 0; // determinant is 0
		}
		add_row(T.elements[0], T.elements[i], 1, T.cols); // add row #i to row #0
	}

	for (int t = 1; t < T.rows; t++) // make first column of T[t] zero
		add_row(T.elements[t], T.elements[0], -T.elements[t][0] / T.elements[0][0], T.cols);

	// Drop the first row and column of T and store in B
	B = allocate_matrix(T.rows-1, T.cols-1);
	copy_matrix(B, 0, T, 1, T.rows-1, T.cols-1);

	/* The determinant of A equals the determinant of B times
	 * T[0][0].
	 */
	det = T.elements[0][0] * determinant(B);

	// free space and return
	free_matrix(T);
	free_matrix(B);
	return det;
}

