/*
 * mtrxmath - inverse.c
 * 
 * Jeff Craig (luserjeff@linuxfreemail.com)
 *
 * This file contains the following functions:
 *   reduce_matrix - makes a matrix a smaller size
 *   determinant - computes the determinant of a matrix
 *   inverse - computes the inverse of a matrix
 */

#include "mtrxmath.h"

/*
 * reduce_matrix
 * Reduces a matrix in size
 * Based on the "Cross Product of Two Vectors in Space"
 * Calc Theorem
 * I believe it has been fixed
 */

MATRIX_PTR reduce_matrix (int cut_row, int cut_column, MATRIX_PTR matrix)
{
  MATRIX_PTR reduced;
  int y, x, rc, rr;
  rc=rr=0;
  
  reduced = alloc_matrix(matrix->rows-1, matrix->rows-1);

  /* This function performs a fairly simple function. It
     reduces a matrix in size by one element on each dimesion
     around the coordinates sent in the cut_row and cut_column
     values.  For example:
        3x3 Matrix:  1 2 3
	             4 5 6
		     7 8 9
	is sent to this function with the cut_row == 2 and 
	cut_column == 1, the function returns the 
	2x2 Matrix: 2 3
	            8 9
  */

  for ( x=0 ; x < matrix->rows ; x++) {
    for ( y=0 ; y < matrix->columns; y++) {
      if( x == cut_row || y == cut_column ) { }
      else {
	reduced->matrix[rr][rc] = matrix->matrix[x][y];
	rc++;
      }
    }
    if (rr != cut_row || x > cut_row ) { rr++; rc=0; }
  }
  return reduced;
}

/* 
 * Determinant
 * Solve the determinant of a matrix
 *
 * Yes, I know this uses the Cramer's Method of finding a
 * determinate, and that Gaussian would be better.  I'm
 * looking into implementing a Gaussian function, but for
 * now this works. 
 */

float determinant ( MATRIX_PTR matrix ) 
{
  float det=0;
  MATRIX_PTR tmp;
  int i, sign=1;
  
  /* Do I need to explain this? */
  if ( matrix->rows != matrix->columns )
    die("Matrix must be square");

  /* This may never be used, but it's necessary for error
     checking */
  if (matrix->rows == 1)
    return matrix->matrix[0][0];

  /* This may not be necessary, but it keeps the recursive function
     down to one less call, so it speeds things up just a bit */
  if (matrix->rows == 2)
    return (matrix->matrix[0][0]*matrix->matrix[1][1]-matrix->matrix[0][1]*matrix->matrix[1][0]);

  /* This is the recursive algorithm that does most of the computation. */
  else {
    for ( i=0 ; i < matrix->rows ; i++, sign*=-1 ) {
      tmp = reduce_matrix( 0,i,matrix );
      det += sign * matrix->matrix[0][i] * determinant(tmp);
      free_matrix(tmp);
    }
    return det;
  }
  return 0;
}

/*
 * inverse
 * 
 * Function to find the inverse of a matrix
 */

MATRIX_PTR inverse ( MATRIX_PTR matrix )
{
  MATRIX_PTR inverse;
  MATRIX_PTR tmp;
  float det=0;
  int row, col, sign=1;
  
  if ( matrix->rows != matrix->columns ) { die("Matrix must be Square"); }

  inverse = alloc_matrix(matrix->rows,matrix->columns);

  det = determinant(matrix);
  if ( det == 0 ) { die("Matrix not Invertable\n"); }
  
  for ( row=0; row<matrix->rows; row++ ) {
    for ( col=0; col<matrix->columns; col++ ) {
      /* 
	 This looks kind of confusing.  All it does is take the 
	 inverse and multiply each square by the determinant that
	 is reduced around the spot the matrix is being reduced by
	 For Instance:

	 In a 3x3 Matrix: A B C
	                  D E F
			  G H I

	 When computing with Element B, Element B is multiplied by
	 the determinant of the 2x2 matrix: D F
	                                    G I
      */
      tmp = reduce_matrix(row,col,matrix);
	if ((row+col)%2 == 0) sign = 1;
	else sign = -1;
      inverse->matrix[col][row] = sign * determinant(tmp);
      free_matrix(tmp);
      
      }
    }

    scalar_mult(inverse,1/det);
    return inverse;
  
  return NULL;
}