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

typedef double elemento;

/* In questo modo se dovremo riusare il codice per allocare ad esempio delle matrici di 
long double bastera' cambiare la definizione di "elemento" come long double e tutto il 
resto del codice sara' riutilizzabile */

typedef elemento ** Matrice;

Matrice CreaMatrice(const int nrighe, const int ncol);
Matrice RiempiMatrice(FILE *ifp,const int nrighe,const int ncol);
elemento Determinante(const Matrice mat, const int nrighe, const int ncol);
void DeallocaMatrice(Matrice a, const int nrighe);
 
int main(int argc, char *argv[])
{

  int n=0;
  FILE *ifp;
  Matrice M;
  elemento det;

  if (argc<2)
    {
	fprintf(stderr,"Uso: >> determinante nome_file_dati \n");
	exit(1);
    }


  ifp=fopen(argv[1],"r");
  if(ifp == NULL){
	fprintf(stderr,"Impossibile aprire il file %s \n",argv[1]);
	exit(1);
  }

  /* In vista di un'estensione a matrici n X n chiede in inpt la dimensione della matrice */
  printf("Inserisci la dimensione della matrice\n");
  scanf("%d",&n);
  

  /* Riempie la matrice M con i punti letti dal file in input;
     la matrice viene allocata e inizializzata in RiempiMatrice
   */
  M = RiempiMatrice(ifp,n,n);
  fclose(ifp);

  /* Calcola il determinante */
  det = Determinante(M,n,n);

  /* Lo stampa in output */
  printf("Il determinante della matrice inserita in %s vale %lf\n",argv[1],det);

  /* E' buona abitudine deallocare le cose allocate quando non le utilizziamo piu' */
  DeallocaMatrice(M,n);

  return(0);
}

Matrice CreaMatrice(const int nrighe,const int ncol){
  Matrice p;
  int i;

  /* Alloca il puntatore alla matrice, cioe' ad un vettore di puntatori, uno per ciascuna riga */

  p=(Matrice)malloc(nrighe*sizeof(elemento*));

  /* Controlla che l'allocazione abbia avuto successo */
  
  if(p==NULL) {
    fprintf(stderr,"Non ho abbastanza spazio per allocare una matrice\n");
    exit(1);
  }
  for(i=0; i<nrighe;i++){

    /* Per ogni riga allochiamo spazio per gli elementi della matrice, inizializzandoli a zero */

    p[i]=(elemento*)calloc(ncol,sizeof(elemento));

    /*...e controlliamo che lo spazio ci sia effettivamente ...*/

    if(p[i]==NULL) {
      fprintf(stderr,"Non ho abbastanza spazio per allocare la riga %d di una matrice\n",i);
      exit(1);
    }
  }
  return(p);
} 


Matrice RiempiMatrice(FILE *ifp,const int nrighe,const int ncol){
  int i,j;
  Matrice ma;
  /* Allocazione dinamica della matrice */
  ma=CreaMatrice(nrighe,ncol);

  /* Copia dati dai file */

  for(i=0;i<nrighe;i++)
    {
      for(j=0;j<ncol;j++)
	{
	  fscanf(ifp,"%lf", &ma[i][j]);
	  if(feof(ifp)){ 
	    fprintf(stderr,"Il file non contiene abbastanza elementi !\n");
	    exit(1);

	  }
	}
    }
  /* IMPORTANTE: "riavvolgiamo" il file in modo da poter riutilizzare il puntatore */
  rewind(ifp);
  return(ma);
}
void DeallocaMatrice(Matrice a, const int nrighe){
  int i=0;
  /* Dealloca ***PRIMA*** i puntatori alle righe...*/
  for(i=0;i<nrighe;i++){
    free(a[i]);
  }
  /* .. e ***POI*** quello a tutta la matrice...*/ 
  free(a);
}

elemento Determinante(const Matrice M, const int nrighe, const int ncol){
  int i=0,j=0,k=0;
  elemento det=0;
  Matrice * Minori;
  /* Questa funzione usa la ricorsione per calcolare il determinante di una  
     matrice di dimensione generica n X n. Per lo svolgimento del compito di esame
     era sufficiente implementare l'algoritmo per il determinante 3 X 3.
     Il caso generale (default) e' inserito in questa soluzione per chi volesse
     approfondire problemi del genere.
  */
  if(nrighe != ncol){
    fprintf(stderr,"Si puo' calcolare il determinante solo di una matrice quadrata!\n");
    exit(1);
  }
  
  switch(nrighe){
  case 0:
    fprintf(stderr,"La matrice deve contenere almeno un elemento!\n");
    exit(1);
    break;
  case 1:
    det = M[0][0];
    break;
  case 2:
    det = M[0][0]*M[1][1]-M[0][1]*M[1][0];
    break;
  case 3:
    det = M[0][0]*(M[1][1]*M[2][2]-M[1][2]*M[2][1]) +
      M[0][1]*(M[1][2]*M[2][0]-M[1][0]*M[2][2]) +
      M[0][2]*(M[1][0]*M[2][1]-M[1][1]*M[2][0]);
    break;
  default:
    /* Prima di tutto alloco un vettore di matrici (sic!)
       una per ogni elemento della prima riga di M, per contenere i 
       minori
    */
    Minori = (Matrice *)malloc(ncol*sizeof(Matrice));

    /* Poi creo e riempio ciascuna matrice */

    for(i=0;i<ncol;i++){
      Minori[i] = CreaMatrice(nrighe-1,ncol-1);

      /* Sto creando i minori ripsetto alla prima riga, quindi devo
	 riempire la prima riga del minore i-simo con la seconda della
	 matrice originale etc. percio' parto da j=1 (seconda riga della
	 matrice originale, e metto gli elementi nella j-1 riga del minore.
       */
      for(j=1;j<nrighe;j++){
	/*
	  Per calcolare l'i-simo minore devo saltare la i-sima colonna, 
	  cosa che realizzo facendo un loop su k che varia su tutte le 
	  colonne della matrice originaria, ma riempendo il minore solo
	  se k e' diverso da i. 
	 */
	for(k=0;k<ncol;k++){
	  if(k<i){
	    Minori[i][j-1][k] = M[j][k];
	      }
	  else if(k>i){
	    Minori[i][j-1][k-1] = M[j][k];
	  }
	}
      }
    }
    /* Quindi calcolo il determinante ricorsivamente: 
       invoco ancora la funzione Determinante, ma con dimensione minore;
       non appena ho raggiunto dei minori di dimensione pari a 3 la 
       ricorsione termina, essendomi ricondotto ad un "caso base" che 
       non necessita di invocare ulteriormente Determinante ma puo' essere
       calcolato direttamente. 
     */

    for(i=0;i<ncol;i++){
      det += M[0][i]*pow(-1,i)*Determinante(Minori[i],nrighe-1,ncol-1);
    }

    /* Se non ripuliamo la memoria adesso perdiamo le "chiavi di accesso" all'area 
       da noi allocata, che restara' bloccata fino al termine del programma 
       (memory leak). Quindi, subito dopo aver calcolato det, che e' l'unica cosa
       che ci interessa, deallochiamo prima ciascuna matrice del vettore minori,
       e poi il puntatore a minori stesso, con la stessa logica con cui deallochiamo
       le matrici (semplice no ? ;-))
     */
    for(i=0;i<ncol;i++){
      DeallocaMatrice(Minori[i],nrighe-1);
    }
    free(Minori);
    break;
  }
  return(det);
}