/************************************************************
 !*  This program calculates the determinant of a complex    *
 !* square matrix using Function CFindDet (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, 0.000000)                     *
 !*                                                          *
 !* -------------------------------------------------------- *
 !*                                                          *
 !*                      C++ Version By J-P Moreau, Paris.   *
 !***********************************************************/ 
#include <stdio.h>
#include <math.h>

#define NMAX  25
#define FALSE  0
#define TRUE   1

  //complex number
  typedef struct {
     double R,I;  //algebraic form
  } COMPLEX;
 
  typedef COMPLEX Matc[NMAX][NMAX];

  int N;
  Matc A;
  COMPLEX D;


  double CABS(COMPLEX Z) {
    return sqrt(Z.R*Z.R + Z.I*Z.I);
  }

  void CADD(COMPLEX Z1, COMPLEX Z2, COMPLEX *Z) {
    Z->R = Z1.R + Z2.R;
    Z->I = Z1.I + Z2.I;
  }

  void CDIF(COMPLEX Z1, COMPLEX Z2, COMPLEX *Z) {
    Z->R = Z1.R - Z2.R;
    Z->I = Z1.I - Z2.I;
  }

  void CMUL(COMPLEX Z1, COMPLEX Z2, COMPLEX *Z) {
    Z->R = Z1.R*Z2.R - Z1.I*Z2.I;
    Z->I = Z1.R*Z2.I + Z1.I*Z2.R;
  }

  void CDIV(COMPLEX Z1, COMPLEX Z2, COMPLEX *Z) {
    double d; COMPLEX C;
    d = Z2.R*Z2.R+Z2.I*Z2.I;
    if (d<1e-12)
      printf(" Complex Divide by zero!\n");
    else {
      C.R=Z2.R; C.I=-Z2.I;
      CMUL(Z1,C,Z);
      Z->R=Z->R/d; Z->I=Z->I/d;
    }
  }
 
 /*---------------------------------------------------------------------------------------------------
 !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
 !----------------------------------------------------------------------------------------------------*/
 void CFindDet(Matc matrix, int n, COMPLEX *cdet) {
     COMPLEX m, temp, CONE, CZERO;
     int i, j, k, l;
     bool DetExists = TRUE;
     l = 1;
	 CZERO.R=0,0; CZERO.I=0.0;
	 CONE.R=1,0; CONE.I=0.0;
     // Convert to upper triangular form
     for (k = 1; k<n; k++) {
         if (matrix[k][k].R == CZERO.R && matrix[k][k].I == CZERO.I) {
             DetExists = FALSE;
             for (i = k+1; i<=n; i++) {
                 if (matrix[i][k].R != CZERO.R || matrix[i][k].I != CZERO.I) {
                     for (j = 1; j<=n; j++) {
                         temp.R = matrix[i][j].R;
                         temp.I = matrix[i][j].I;
                         matrix[i][j].R = matrix[k][j].R;
                         matrix[i][j].I = matrix[k][j].I;
                         matrix[k][j].R = temp.R;
                         matrix[k][j].I = temp.I;
                     }
                     DetExists = TRUE;
                     l=-l;
                     break;
                 }
             }
             if (DetExists == FALSE) {
                 cdet->R = CZERO.R;
                 cdet->I = CZERO.I;
                 return;
             }
         }
         for (j = k+1; j<=n; j++) {
             // m = matrix(j,k)/matrix(k,k)
			 CDIV(matrix[j][k],matrix[k][k],&m);
             for (i = k+1; i<=n; i++) {
                 // matrix(j,i) = matrix(j,i) - m*matrix(k,i)
                 CMUL(m,matrix[k][i],&temp);
		 matrix[j][i].R = matrix[j][i].R - temp.R;
                 matrix[j][i].I = matrix[j][i].I - temp.I;
             }
         }
     }
     
     // Calculate determinant by finding product of diagonal elements
     cdet->R = 1.0; cdet->I = 0.0;
     for (i = 1; i<=n; i++) {
         // CDet = CDet * matrix(i,i)
	 CMUL(*cdet,matrix[i][i],&temp);
	 cdet->R = temp.R;
         cdet->I = temp.I;
     }
     
}

void main() {
 int N=4;

 A[1][1].R = 0.0; A[1][1].I = 0.0;
 A[2][1].R = 1.0; A[2][1].I = 0.0;
 A[3][1].R = 1.0; A[3][1].I = 0.0;
 A[4][1].R = 1.0; A[4][1].I = 0.0;

 A[1][2].R = -1.0; A[1][2].I = 0.0;
 A[2][2].R =  0.0; A[2][2].I = 0.0;
 A[3][2].R =  1.0; A[3][2].I = 0.0;
 A[4][2].R =  1.0; A[4][2].I = 0.0;

 A[1][3].R = -1.0; A[1][3].I = 0.0;
 A[2][3].R = -1.0; A[2][3].I = 0.0;
 A[3][3].R =  0.0; A[3][3].I = 0.0;
 A[4][3].R =  1.0; A[4][3].I = 0.0;

 A[1][4].R = -1.0; A[1][4].I = 0.0;
 A[2][4].R = -1.0; A[2][4].I = 0.0;
 A[3][4].R = -1.0; A[3][4].I = 0.0;
 A[4][4].R =  0.0; A[4][4].I = 0.0;

 
 CFindDet(A, N, &D);


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

}

// end of file cfinddet.cpp