/************************************************************
 !* This program calculates the determinant of a real square *
 !* matrix using Function FindDet (see ref. below).          *
 !* -------------------------------------------------------- *
 !* SAMPLE RUN:                                              *
 !*                                                          *
 !* Calculate the deteminant of square matrix:               *
 !*                                                          *
 !*           ! 0  -1  -1  -1 !                              *
 !*           ! 1   0  -1  -1 !                              *
 !*           ! 1   1   0  -1 !                              *
 !*           ! 1   1   1   0 !                              *
 !*                                                          *
 !*  Determinant =  1.000000                                 *
 !*                                                          *
 !* -------------------------------------------------------- *
 !*              C++ Version By Jean-Pierre Moreau, Paris.   *
 !***********************************************************/ 
#include <stdio.h>
#include <math.h>

#define N  4

#define FALSE  0
#define TRUE   1

typedef double MAT[N+1][N+1];  // index 0 not used here

 /*---------------------------------------------------------------------------------------------------
 !Function to find the determinant of a square matrix
 !Author : Louisda16th a.k.a Ashwith J. Rego
 !Description: The subroutine is based on two key points:
 !1] A determinant is unaltered when row operations are performed: Hence, using this principle,
 !row operations (column operations would work as well) are used
 !to convert the matrix into upper traingular form
 !2]The determinant of a triangular matrix is obtained by finding the product of the diagonal elements
 !--------------------------------------------------------------------------------------------------*/
 double FindDet(MAT matrix, int n) {
     double m, tmp;
     int i, j, k, l;
     bool DetExists = TRUE;
     l = 1;
     // Convert to upper triangular form
     for (k = 1; k<n; k++) {
         if (matrix[k][k] == 0) {
             DetExists = FALSE;
             for (i = k+1; i<=n; i++)
                 if (matrix[i][k] != 0) {
                     for (j = 1; j<=n; j++) {
                         tmp = matrix[i][j];
                         matrix[i][j]= matrix[k][j];
                         matrix[k][j] = tmp;
                     }
                     DetExists = TRUE;
                     l=-l;
                     goto e10;
                 }
e10:         if (DetExists == FALSE) {
                 return 0.0;
             }
         }
         for (j = k+1; j<=n; j++) {
             m = matrix[j][k]/matrix[k][k];
             for (i = k+1; i<=n; i++)
                 matrix[j][i] -= m*matrix[k][i];
         }
     } //end of k loop
     
     // Calculate determinant by finding product of diagonal elements
     tmp = l;
     for (i = 1; i<=n; i++)
         tmp *= matrix[i][i];

	 return tmp;
     
 } // FindDet()

 
void main() {

 MAT A;
 double D;

 A[1][1] =  0; A[1][2] =  1; A[1][3] =  1; A[1][4] = 1;
 A[2][1] = -1; A[2][2] =  0; A[2][3] =  1; A[2][4] = 1;
 A[3][1] = -1; A[3][2] = -1; A[3][3] =  0; A[3][4] = 1;
 A[4][1] = -1; A[4][2] = -1; A[4][3] = -1; A[4][4] = 0;

 D = FindDet(A, N);

 printf("\n Determinant = %f\n\n", D);

}

// end of file TFindDet.cpp