/* Kelvin Cheung
   1.12
   PS5
*/

#include<iostream.h>
#include<stdlib.h>
#include<math.h>
#include<fstream.h>

class Matrix {
  friend class Solver;
protected:
  float **numbers;
  int rows,cols;
public:
  Matrix() {}
  Matrix(ifstream&,int,int);
  Matrix(int,int);
  Matrix(Matrix,int,int);
  Matrix& operator+(const Matrix&);
  Matrix& operator-(const Matrix&);
  Matrix& operator*(const Matrix&);
  float determinant();
  void row_multiply_and_add(int,float,int);
  void row_multiply(int,float);
  Matrix inverse();
  void print() {
    for(int i=0;i<rows;i++) {
      for(int j=0;j<cols;j++)
	cout << numbers[i][j] << "   ";
      cout << endl;
    }
  }
};

class Vector : public Matrix {
friend class Solver;
private:
  int size;
public:
  Vector () {}
  Vector(ifstream&,int);
  Vector(int);
};

class Solver : public Matrix {
public:
  Solver(Matrix,Vector);
  void gauss();
};

Matrix::Matrix(ifstream &infile, int r, int c)
     // This is the constructor to be used for reading from a text file
     // n is the dimension of the matrix
{
  rows=r;
  cols=c;

  numbers=new float*[r];
  for(int i=0;i<r;i++) {
    numbers[i]=new float[c];
    for(int j=0;j<c;j++)
      infile >> numbers[i][j];
  }
}

Matrix::Matrix(int r, int c)
     // This is a constructor to create a r x c matrix with no values
{
  rows=r;
  cols=c;

  numbers=new float*[rows];
  for(int i=0;i<r;i++)
    numbers[i]=new float[c];
}

Matrix::Matrix(Matrix m, int r, int c)
     // creates a matrix from the m matrix which takes out the r row and
     // c column (used for taking determinants)
{
  rows=m.rows-1;
  cols=m.cols-1;
  numbers=new float*[rows];
  for(int i=0;i<rows;i++)
    numbers[i]=new float[cols];
  for(int i=0;i<r;i++) {
    for(int j=0;j<c;j++)
      numbers[i][j]=m.numbers[i][j];
    for(int j=c+1;j<m.cols;j++)
      numbers[i][j-1]=m.numbers[i][j];
  }
  for(int i=r+1;i<m.rows;i++) {
    for(int j=0;j<c;j++)
      numbers[i-1][j]=m.numbers[i][j];
    for(int j=c+1;j<m.cols;j++)
      numbers[i-1][j-1]=m.numbers[i][j];
  }
}

Matrix& Matrix::operator+(const Matrix &m)
     // overloaded + operator
{
  Matrix sum=Matrix(rows,cols);
  for(int i=0;i<rows;i++)
    for(int j=0;j<cols;j++)
      sum.numbers[i][j]=numbers[i][j]+m.numbers[i][j];
  return sum;
}

Matrix& Matrix::operator-(const Matrix &m)
     // overloaded - operator
{
  Matrix diff=Matrix(rows,cols);
  for(int i=0;i<rows;i++)
    for(int j=0;j<cols;j++)
      diff.numbers[i][j]=numbers[i][j]-m.numbers[i][j];
  return diff;
}

Matrix& Matrix::operator*(const Matrix &m)
     // overloaded * operator
{
  Matrix product=Matrix(rows,m.cols);
  for(int i=0;i<rows;i++)
    for(int j=0;j<m.cols;j++) {
      product.numbers[i][j]=0;
      for(int k=0;k<cols;k++)
	product.numbers[i][j]+=numbers[i][k]*m.numbers[k][j];
    }
  return product;
}

float Matrix::determinant()
     // computes determinant.  It does this recursively.
{
  float det=0;
  Matrix *temp;

  if(rows==1)
    return numbers[0][0];
  for(int i=0;i<cols;i++) {
    temp=new Matrix(*this,0,i);
    if(i%2==0)
      det+=numbers[0][i]*temp->determinant();
    else
      det-=numbers[0][i]*temp->determinant();
  }
  return det;
}

void Matrix::row_multiply(int r, float factor)
     // multiplies row r by factor.
{
  for(int i=0;i<cols;i++)
    numbers[r][i]=numbers[r][i]*factor;
}

void Matrix::row_multiply_and_add(int r1, float factor, int r2)
     // changes row r1 by adding its current value to r2 multiplied by factor.
     // this is used for Gaussian elimination
{
  for(int i=0;i<cols;i++)
    numbers[r1][i]=numbers[r1][i]+numbers[r2][i]*factor;
}

Matrix Matrix::inverse()
     // computes inverse of a square matrix.  It does this with the Gauss-
     // Jordan method which operates on the fact that [A I] multiplied by
     // inv(A) gives you [I inv(A)]
{
  float factor;
  Matrix inv=Matrix(rows,cols);
  Matrix copy=Matrix(rows,cols);

  for(int i=0;i<rows;i++) {
    for(int j=0;j<cols;j++) {
      inv.numbers[i][j]=0;
      copy.numbers[i][j]=numbers[i][j];
    }
    inv.numbers[i][i]=1;
  }

  for(int i=0;i<cols;i++)
    for(int j=i+1;j<rows;j++) {
      factor=-copy.numbers[j][i]/copy.numbers[i][i];
      copy.row_multiply_and_add(j,factor,i);
      inv.row_multiply_and_add(j,factor,i);
    }
  
  for(int i=cols-1;i>=0;i--)
    for(int j=i-1;j>=0;j--) {
      factor=-copy.numbers[j][i]/copy.numbers[i][i];
      copy.row_multiply_and_add(j,factor,i);
      inv.row_multiply_and_add(j,factor,i);
    }

  for(int i=0;i<rows;i++)
    inv.row_multiply(i,1/copy.numbers[i][i]);

  return inv;
}

Vector::Vector(ifstream &infile, int n) : Matrix(infile,n,1)
     // constructor to create the vector from a text file.  n is the size
{
  size=n;
}

Vector::Vector(int n) : Matrix(n,1)
     // constructor that creates a size n vector with no values
{
  size=n;
}

Solver::Solver(Matrix m, Vector v)
     // creates an augmented matrix with the matrix m and the vector v
     // stapled on the end
{
  rows=m.rows;
  cols=m.cols+1;
  numbers=new float*[rows];
  for(int i=0;i<rows;i++) {
    numbers[i]=new float[cols];
    for(int j=0;j<cols;j++)
      numbers[i][j]=m.numbers[i][j];
    numbers[i][cols-1]=v.numbers[i][0];
  }
}

void Solver::gauss()
     // solves the set of equations in the augmented matrix using Gaussian
     // elimination.  it does this by using a few member functions from the
     // Matrix class which were used in the Gauss-Jordan inverse method
{
  Matrix copy(rows,cols);
  Vector soln(rows);

  for(int i=0;i<rows;i++)
    for(int j=0;j<cols;j++)
      copy.numbers[i][j]=numbers[i][j];

  for(int i=0;i<rows;i++)
    for(int j=i+1;j<rows;j++)
      copy.row_multiply_and_add(j,-copy.numbers[j][i]/copy.numbers[i][i],i);
  for(int i=rows-1;i>=0;i--)
    for(int j=i-1;j>=0;j--)
      copy.row_multiply_and_add(j,-copy.numbers[j][i]/copy.numbers[i][i],i);

  for(int i=0;i<copy.rows;i++) {
    copy.row_multiply(i,1/copy.numbers[i][i]);
    cout << copy.numbers[i][cols-1] << endl;
  }
}