/*
  naive matrix determinant calculation using 
  determinant expansion by minors implemented
  as a recursive function ;)
-----------------------------------------------
  mateusz@malczak.info
  http://devblog.malczak.info                    
*/

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

static int mC=0;
static int mD=0;
static long memA = 0;
static long memD = 0;

typedef struct __matrix 
{
 unsigned int N;
 double **data;
} matrix;

matrix* new_matrix(unsigned int n)
{
 matrix* m = (matrix*)malloc( sizeof(matrix) );
 memA += sizeof(matrix);
 m->N = n;
 m->data = (double**)malloc( n*sizeof(double*) );
 memA += n*n*sizeof(double);
  while ( n )  
  {
   m->data[--n] = (double*)malloc( m->N*sizeof(double) ); 
   memset(m->data[n],m->N*sizeof(double),0);
  };
 mC+=1;
 return (matrix *)m;
};

void free_matrix(matrix **m)
{
 memD -= (*m)->N*(*m)->N*sizeof(double);
 memD -= sizeof(matrix);
 while ( (*m)->N )
  free( (*m)->data[--(*m)->N] );
 free( (*m)->data );
 free( (*m) );
 mD+=1;
 *m = NULL; 
}

matrix* M( matrix *mtx, unsigned int M, unsigned int N )
{ 
 matrix* m = new_matrix(mtx->N-1);
 unsigned int i=0,j=0,ii=0,jj=0;
  for ( ; i<mtx->N; i++ )
  {
   if ( i!=M ) 
   {
   	for ( jj=0,j=0; j<mtx->N; j++ ) 
    		if ( j!=N ) 
	     		m->data[ii][jj++] = mtx->data[i][j];     
     ii++;
   };
  }; 
 return (matrix*)m;
}

double D( matrix* mtx )
{
 double det = 0.0;
 unsigned int m=0, n;
  if ( mtx->N>1 )
  {
  	for ( n=0; n<mtx->N; n++ )
	{
	   	det += pow(-1,n)*mtx->data[m][n]*D( M( (matrix*)mtx, m, n) );
	};
  } else det = mtx->data[0][0];
 free_matrix( &mtx );
 return det;
};

int main()
{
 int i,j;
 time_t t1=0, t2=0;
 double detv=0.0;
 printf("Enter matrix dimension : ");
 scanf("%d",&i);
 matrix *m = new_matrix(i);
 srand( time(NULL) );
  
  printf("Macierz wejsciowa\n");
  for (i=0; i<m->N; i++ )
  {
  	for (j=0; j<m->N; j++ )
	{
	  m->data[i][j] = rand()/(double)RAND_MAX;
	  	printf("%.3f ",m->data[i][j]);
	};
			printf("\n");
  };

 t1 = time(NULL);
  detv = D(m);  
 t2 = time(NULL);
 printf("det(M)=%f\n\n",detv); 

 printf("total time = %lf\n",difftime(t2,t1));
 printf("created matrices = %d, destroyed matrices = %d\n",mC,mD);
 printf("memory allocated   = %db ~ %.2fkb ~ %.2fMb\n",memA,memA/1024.0,memA/1024.0/1024.0);
 printf("memory deallocated = %db ~ %.2fkb ~ %.2fMb\n",memD,memD/1024.0,memD/1024.0/1024.0);
 
 return 1;
}