#include <iostream.h>

/* Implements an arbitrary sized matrix class. */

class matrix
{
  public:
  		matrix(const int row, const int col);
      ~matrix();
      long double det();
      long double set(int row, int col, long double val);
      void print();
      void equals(matrix *source);
      void inverse(void);
      long double value(int row, int col);
      void times_row(long double *row_mat, long double *results);

  private:
  		long double *values;
      int row, col;

};

/* Constructor: Initializes number of rows and cols in the matrix and sets
   all elements to zero. */

matrix::matrix(const int num_rows, const int num_cols)
{
	values = new long double[num_rows*num_cols];
   row = num_rows;
   col = num_cols;
   for (int i=0; i<row;i++)
   	{
      for (int j=0; j<col; j++)
      	set(i,j,0.0);
      }
}

/* Destructor: deletes the matrix. */
matrix::~matrix(void)
{
	delete [] values;
}

/* Sets the value of the matrix to val at the row, r, and column, c.  Returns
	val. */
long double matrix::set(int r, int c, long double val)
{
   *(values + row*c + r) = val;
   return val;
}

/* Takes the determinant of a SQUARE matrix recursively, using the cofactor
   method. */
long double matrix::det()
{
if (row == 2 && col == 2) // Base case for 2x2 matrix.
   return (*(values) * *(values+3) - *(values+1) * *(values+2));

else  // Recursive case
	{
   int t_col;                       // tracks column of temp that is being set
   long double the_sign = 1.0;      // tracks sign of cofactor
   long double result = 0.0;
   matrix temp(row-1, col-1);       // a temporary cofactor of the matrix

   for (int i=0; i<col; i++)        // loop over all necessary cofactors
   	{
      t_col = 0;
      for (int j=0; j<col; j++)     // double loop to set values of cofactors
      	{									// excluding the proper row and column
         if (i!=j)
            {
         	for (int k=1; k<row; k++)
	         	temp.set(k-1,t_col,*(values + row*j + k) );
            t_col++;
            }
         }

      // Determinant is sum over columns of:
      //   element in row=0,col=i * det(cofactor) with proper sign
      result += *(values + row*i) * temp.det()*the_sign;
      the_sign *= -1.0; // the sign simply flips as one moves along row=0
      }
   return (result);
   }

}

/* Prints matrix to screen, tab delimited. */
void matrix::print()
{
 for (int i=0; i<row; i++)
 	{
   for (int j=0; j<col; j++)
   	cout << *(values + row*j + i) << '\t';
   cout << endl;
   }
 cout << endl;
}

/* Set the calling matrix to the values of source.  Matrices must be same size
   or strange things will occur. */
void matrix::equals(matrix *source)
{
	for (int i=0;i<row*col;i++)
   	*(values + i) = *(source->values+i);
}

/* Takes the inverse of a matrix using Gauss-Jordan elimination. */
void matrix::inverse()
{
matrix old(row,col),temp(row, col),temp_old(row,col);
long double subval,divval;

// Create temp = unity matrix
for (int i=0;i<row;i++) temp.set(i,i,1.0);

// The inverse of the matrix is the unity matrix divided by the matrix.
// This method converts the orginal matrix to the unity matrix using
// eliminations.  It performs the same operations of the unity matrix.
// Therefore, at the end of the procedure the original matrix is the unity
// matrix, and the unity matrix is the inverse of the original matrix.

for (int i=0; i<col; i++) // eliminate all elements in col i, except row=i,col=i
	{
   old.equals(this);       // save old values temporarily
   temp_old.equals(&temp);

	for (int k=0; k<col; k++)
		{
 		for (int j=0; j<row; j++)
			{
     		if (j!=i)
        		{
	         subval = old.value(j,i) / old.value(i,i);
  		      set(j,k, old.value(j,k) - old.value(i,k)*subval );
     	      if (k<=i)
  				temp.set(j, k, temp_old.value(j,k) - temp_old.value(i,k)*subval );
         	}
      	}
	   }
   for (int n=0;n<col;n++)  // normalize remaining value to 1
   	{
      divval = 1.0/old.value(i,i);
      set(i,n,value(i,n)*divval );
      temp.set(i,n, temp.value(i,n) * divval );
      }

   }

equals(&temp);

}

/* returns the value of the matrix at row=r and column=c */
long double matrix::value(int r, int c)
{
  return ( *(values + row*c + r));
}

/* Multiplies a row matrix by the matrix calling the procedure.
   It assumes the row matrix has "row" number of elements and the results
   row matrix has "row" number of elements. */
void matrix::times_row(long double *row_mat, long double *results)
	{
   for (int i = 0; i<col; i++)
   	{
      results[i] = 0.0;
      for (int j=0; j<row; j++)
      	{
         results[i] += row_mat[j] * value(j,i);
         }
      }
   }