#include "determinant.h"

#include "number.h"
#include "rational.h"
#include "definitions.cpp"

// determinant of a matrix
// replaces the matrix by its LU decomposition
int determinant_inplace(Number_Matrix M, Number* out){
  int i;
  int sign_correction;

  LU_dcmp_inplace(M, &sign_correction);

  if(sign_correction==0){
    *out=number_zero();
    return(0);
  }

  *out=number_one();
  if(sign_correction==-1){
    number_Qprod_chain(quot(-1,1), out);
  }

  for(i=0;i<M.length;i++){
    number_prod_chain(M.matrix[i][i], out);
  }

  return(0);
}

// LU decomposition
// uses pivoting to avoid dividing by 0
// the sign_correction should be multiplied to the determinant to obtain the right value
// if dividing by 0 is unavoidable, then the determinant is 0, and sign_correction is set to 0
int LU_dcmp_inplace(Number_Matrix M, int* sign_correction){
  int i,j,k,pivot;
  Number tmp;

  *sign_correction=1;

  for(j=0;j<M.length;j++){
    for(i=0;i<=j;i++){
      for(k=0;k<i;k++){
	// -M[i][k]*M[k][j]
	number_prod(M.matrix[i][k], M.matrix[k][j], &tmp);
	number_Qprod_chain(quot(-1,1), &tmp);
	number_add_chain(tmp, M.matrix[i]+j);
	free_Number(tmp);
      }
    }
    for(i=j+1;i<M.length;i++){
      for(k=0;k<j;k++){
	// -M[i][k]*M[k][j]
	number_prod(M.matrix[i][k], M.matrix[k][j], &tmp);
	number_Qprod_chain(quot(-1,1), &tmp);
	number_add_chain(tmp, M.matrix[i]+j);
	free_Number(tmp);
      }
    }

    // pivot if M[j][j]==0
    // find first M[j][j] that is not 0
    for(pivot=j;pivot<M.length && number_is_zero(M.matrix[pivot][j])==1;pivot++){}

    // no non-zero M[j][j] left: return
    if(pivot>=M.length){
      *sign_correction=0;
      return(0);
    }
    // pivot if needed
    if(pivot!=j){
      for(k=0;k<M.length;k++){
	tmp=M.matrix[j][k];
	M.matrix[j][k]=M.matrix[pivot][k];
	M.matrix[pivot][k]=tmp;
      }
      *sign_correction*=-1;

    }

    for(i=j+1;i<M.length;i++){
      // do not use the inplace algorithm if M[j][j] has more than one terms, since it would be modified by the inplace function
      if(M.matrix[j][j].length<=1){
	number_quot_inplace(M.matrix[i]+j, M.matrix[j]+j);
      }
      else{
	number_quot_chain(M.matrix[i]+j, M.matrix[j][j]);
      }
    }
  }
  return(0);
}