DETERMINANTS (GAUSS & COFACTEURS)

//-------------------------------------------------
// calcul d'un determinant
//-------------------------------------------------

#ifndef _UTIL_H_
	#include "util.h"
#endif // _UTIL_H_

#ifndef _MATH_H_
	#include "math.h"
#endif // _MATH_H_

#include <conio.h>

//-------------------------------------------------
typedef double FLT;
//-------------------------------------------------
FLT **CreateMat(int n)
{
FLT **mat; // la matrice
int   y;   // <y> pour les lignes

mat = Malloc(FLT*,n);
for(y=0;y<n;y++)
  {
  mat[y] = Malloc(FLT,n);
  }

return mat;
} // CreateMat()


//-------------------------------------------------
FLT **RandomMat(int n,FLT coeffMin,FLT coeffMax)
{
FLT **mat;
int   x,y;  // <x> pour les colonnes
            // <y> pour les lignes

mat = CreateMat(n);

for(x=0;x<n;x++)
  {
  for(y=0;y<n;y++)
    {
    mat[y][x] = (FLT)RandomDouble(coeffMin,coeffMax);
    }
  }

return mat;
} // RandomMat()


//-------------------------------------------------
void DeleteMat(FLT **mat,int n)
{
int y; // <y> pour les lignes
for(y=0;y<n;y++)
  {
  Free(mat[y]);
  }
Free(mat);
} // DeleteMat()


//-------------------------------------------------
// calcul d'un determinant avec le pivot de Gauss
// la complexite est en n^3
FLT FastDet(FLT **_mat,int n)
{
FLT **mat;  // la matrice "tampon", que l'on pourra modifier
FLT   det;  // la valeur du determinant
int   x,y;  // <x> pour les colonnes
            // <y> pour les lignes
int   i,j;  // <i> pour les colonnes
            // <j> pour les lignes

// on copie <_mat> dans <mat> car on n'a pas le droit
// de modifier les valeurs de la matrice <_mat>
mat = CreateMat(n);
for(y=0;y<n;y++)
  {
  for(x=0;x<n;x++)
    {
    mat[y][x] = _mat[y][x];
    }
  }

// on calcule le determinant par la methode des pivots de Gauss
// on rend la matrice triangulaire superieur tout en conservant
// le determinant de la matrice a chaque iteration
det = 1.;

// on balaye les lignes
for(j=0;j<n-1;j++)
  {
  // Le principe est simple, montrons(le sur une exemple 6x6
  // La matrice est sous cette forme :
  // 
  // ( a1 *  *  *  *  * )  [L1]
  // ( 0  a2 *  *  *  * )  [L2]  <-- noms des lignes
  // ( 0  0  a3 *  *  * )  [L3]
  // ( 0  0  0  #1 *  * )  [L4]
  // ( 0  0  0  #2 *  * )  [L5]
  // ( 0  0  0  #3 *  * )  [L6]
  // 
  // Le determinant courant <det> vaut a1.a2.a3
  // Si la matrice est non-singuliere, (i.e. que sont determinant est non-nulle,
  // autrement dit elle est inversible) alors parmis les trois nombres #1 #2 et #3
  // il y en a un qui est non-nul, on prend le plus grand de ceux la en module pour
  // une raison de stabilite numerique.
  // Si #1=#2=#3=0 alors <det> est nul car la matrice est singuliere !
  // Par exemple si le plus grand est #2 alors on echange [L4] et [L5], ce qui va modifier
  // le determinant d'un signe.
  // 
  // La matrice sera donc sous cette forme :
  // 
  // ( a1 *  *  *  *  *  )  [L1]
  // ( 0  a2 *  *  *  *  )  [L2]
  // ( 0  0  a3 *  *  *  )  [L3]
  // ( 0  0  0  a4 *  *  )  [L4] (avec a4=#2)
  // ( 0  0  0  #1 *  *  )  [L5]
  // ( 0  0  0  #3 *  *  )  [L6]
  // 
  // Ensuite il suffit de retrancher #1/a4.[L4] a [L5] et #3/a4.[L4] a [L6] pour obtenir :
  // 
  // ( a1 *  *  *  *  *  )  [L1]
  // ( 0  a2 *  *  *  *  )  [L2]
  // ( 0  0  a3 *  *  *  )  [L3]
  // ( 0  0  0  a4 *  *  )  [L4]
  // ( 0  0  0  0  *  *  )  [L5]
  // ( 0  0  0  0  *  *  )  [L6]
  // 
  // Le determinant est inchange par cette derniere operation.
  // Le determinant courant vaut maintenant a1.a2.a3.a4
  // A la fin de l'operation on aura <mat> :
  // 
  // ( a1 *  *  *  *  *  )  [L1]
  // ( 0  a2 *  *  *  *  )  [L2]
  // ( 0  0  a3 *  *  *  )  [L3]
  // ( 0  0  0  a4 *  *  )  [L4]
  // ( 0  0  0  0  a5 *  )  [L5]
  // ( 0  0  0  0  0  a6 )  [L6]
  // 
  // Remarquons qu'il suffit de faire l'operation <n-1> fois
  // car a la fin on aura a6 directement
  // 
  // En resume voici les differentes etapes :
  //   ( etape 1 ) : on recherche du plus grand coefficient (en module)
  //                 si il est nul, alors le determinant est 0 et l'algorithme est fini !
  //   ( etape 2 ) : on echange des lignes avec eventuellement changement du signe de <det>
  //                 et le determinant courant est multiplie par ledit coefficient
  //   ( etape 3 ) : on retranche la lignes a toutes les lignes du dessous
  //                 pour annuler les coefficients
  int  rankMax,rank;
  FLT  coeffMax;

  // ( etape 1 )
  rankMax = j;
  for(rank=j+1;rank<n;rank++)
    {
    if(fabs(mat[rankMax][j]) < fabs(mat[rank][j]))
      {
      rankMax = rank;
      }
    }

  coeffMax = mat[rankMax][j];
  if(fabs(coeffMax) <= MATH_EPSILON)
    {
    det = 0.;
    goto label_end;
    }
  // ( etape 2 )
  if(rankMax != j)
    {
    for(i=j;i<n;i++)
      {
      FLT tmp;
      tmp             = mat[j][i];
      mat[j][i]       = mat[rankMax][i];
      mat[rankMax][i] = tmp;
      }
    det *= -1.;
    }

  det *= coeffMax;
  // ( etape 3 )
  for(rank=j+1;rank<n;rank++)
    {
    FLT coeff; // c'est dans l'exemple le rapport #3/a4
    coeff = mat[rank][j]/coeffMax;
    for(i=j;i<n;i++)
      {
      mat[rank][i] -= coeff*mat[j][i];
      }
    }
    
  }

det *= mat[n-1][n-1];

// on libere <mat>
label_end:
DeleteMat(mat,n);
return det;
} // FastDet()


//-------------------------------------------------
// calcul d'un determinant avec les cofacteurs
// la complexite est en n!
FLT StupidDet(FLT **mat,int n)
{
// cas trivial d'une matrice 1x1
if(1 == n)
  {
  return mat[0][0];
  }
// cas trivial d'une matrice 2x2
// 
// | a b |
// | c d |  = a.d - c.d
// 
else if(2 == n)
  {
  return mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0];
  }
else
  {
  FLT    **cof; // matrice qui represente un cofacteur
  FLT      det; // le determinant
  int      x,y; // <x> pour les colonnes
                // <y> pour les lignes
  int      i,j; // <i> pour les colonnes
                // <j> pour les lignes
  FLT      sgn;

  // on cree une matrice temporaire pour le calcul des cofacteurs
  cof = CreateMat(n-1);
  for(y=0;y<n-1;y++)
    {
    for(x=0;x<n-1;x++)
      {
      cof[y][x] = mat[y][x];
      }
    }

  // quelques initialisations
  det =  0.0;
  sgn = +1.0;
  // on balaye la derniere colonne (de bas en haut)
  for(j=n-1;j>=0;j--)
    {
    // Le principe est tres simple, il s'agit de developper le determinant
    // suivant la derniere colonne. Cet algorithme est evidemment recursif
    // puisqu'il faut les sous-determinants
    // 
    // Exemple sur une matrice 4x4 :
    // 
    // | *  *  *  a1 |
    // | *  *  *  a2 |
    // | *  *  *  a3 |
    // | *  *  *  a4 |
    // 
    // A chaque coefficient de la derniere colonne, on associe une sous-dertermiant ("sd") :
    // (notation : '#')
    // 
    // ds1 : | *  *  *  a1 |      ds2 : | #  #  #  a1 |
    //       | #  #  #  a2 |            | *  *  *  a2 |
    //       | #  #  #  a3 |            | #  #  #  a3 |
    //       | #  #  #  a4 |            | #  #  #  a4 |
    // 
    // ds3 : | #  #  #  a1 |      ds4 : | #  #  #  a1 |
    //       | #  #  #  a2 |            | #  #  #  a2 |
    //       | *  *  *  a3 |            | #  #  #  a3 |
    //       | #  #  #  a4 |            | *  *  *  a4 |
    // 
    // La formule du calcul du determinant par les cofacteurs est donnee par :
    // 
    //     det = a4.ds4 - a3.ds3 + a2.ds2 - a1.ds1
    //
    //     (notons l'alternance des signes)
    // 
    // On doit calculer les "ds", donc creer les matrices correspondants, et calculer
    // par recurrence leurs determinants. On remarque facilement que les matrices se
    // ressemblent à (toujours) 2 lignes pres. Donc pour eviter de toujours les recalculer,
    // on les initilise de facon "incrementale".
    // 
    
    // on modifie la bonne ligne pour la matrice des cofacteurs
    if(j < n-1)
      {
      for(i=0;i<n-1;i++)
        {
        cof[j][i] = mat[j+1][i];
        }
      }
    // la recurrence en personne !
    det += sgn * mat[j][n-1] * StupidDet(cof,n-1);
    // changement de signe (l'alternance remarquee)
    sgn *= -1.0;
    }

  DeleteMat(cof,n-1);
  return det;
  }
} // StupidDet()