#include "matrix.h"

typedef double* DP;

#define NP(i) ((((i)%2)==0)? 1 : -1)


Matrix::Matrix(int row, int col) {
    mp=(double**)new DP [row];
    assert(mp!=0);
    for (int i=0;i<row;i++) {
        mp[i]=new double [col];
        assert(mp[i]!=0);
    }
    rows=row; cols=col;
    return;
};

Matrix::Matrix( Matrix & other) {
    int i, j, row=other.rows, col=other.cols;
    mp=(double**)new DP  [row];
    assert(mp!=0);
    for (i=0;i<row;i++) {
        mp[i]=new double[col];
        assert(mp[i]!=0);
    }
    for (i=0; i<row; i++)
        for (j=0; j<col; j++)
            mp[i][j]=other.mp[i][j];
    rows=row; cols=col;
    return;
};

Matrix::~Matrix() {
    for(int i=0; i<rows; i++) delete mp[i];
    delete mp;
    mp=0;
};

Matrix& Matrix::operator=( Matrix& other) {
    int i, j, row=other.rows, col=other.cols;
    mp=(double**)new DP  [row];
    assert(mp!=0);
    for (i=0;i<row;i++) {
        mp[i]=new double[col];
        assert(mp[i]!=0);
    }
    for (i=0; i<row; i++)
        for (j=0; j<col; j++)
            mp[i][j]=other.mp[i][j];
    return *this;
};

double& Matrix::operator()(int i, int j) {
    assert((0<=i)&&(i<=rows)&&(0<=j)&&(j<=cols));
    return mp[i-1][j-1];
};

Matrix operator+(Matrix& A, Matrix& B) {
    int i,j;
    Matrix res(A.rows, A.cols);
    assert((A.rows==B.rows)&&(A.cols==B.cols));
    for (i=1; i<=A.rows; i++)
        for (j=1; j<=A.cols; j++)
            res(i,j)=A(i,j)+B(i,j);
    return res;
};

Matrix operator*( Matrix& A,  Matrix& B) {
    int i,j,k;
    Matrix res(A.rows, B.cols);
    assert(A.cols==B.rows);
    for (i=1; i<=A.rows; i++)
        for (j=1; j<=B.cols; j++) {
            res(i,j)=0;
            for (k=1;k<=A.cols;k++)
                res(i,j)+=A(i,k)*B(k,j);
        }
    return res;
};

Matrix operator*(Matrix& A, double c) { return c*A;};

Matrix operator*(double c,  Matrix& A) {
    Matrix res(A);
    int i, j;
    for (i=1; i<=A.rows; i++)
        for (j=1; j<=A.cols; j++)
            res(i,j)*=c;
    return res;
};

istream& operator>>(istream& ist,  Matrix& A) {
    int i, j;
    for (i=1; i<=A.rows; i++)
        for (j=1; j<=A.cols; j++) {
            cout << "Enter (" << i << "," << j << ")'th entry: ";
            ist >> A(i,j);
        }
    return ist;
};

ostream& operator<<(ostream& ost,  Matrix& A) {
    int i, j;
    for (i=1; i<=A.rows; i++) {
        ost << "[ ";
        for (j=1; j<=A.cols; j++)
            cout << A(i,j) << "\t";
        ost << " ] " << endl;
    }
    return ost;
};

double& Vector::operator()(int row) {
    return *(mp[row-1]);
};

double operator|( Vector& first,  Vector& second      ) {
    assert(first.rows==second.rows);
    double result=0;
    for(int i=1; i<=first.rows; i++) result+=first(i)*second(i);
    return result;
};

double SquareMatrix::tr() {
    double result=0;
    for (int i=1; i<=rows; i++) result+=(*this)(i,i);
    return result;
};

Matrix Matrix::restrict(int i, int j) {
    Matrix res(rows-1, cols-1);
    for (int r=1; r<i; r++) {
        for (int c=1; c<j; c++)
            res(r,c)=(*this)(r,c);
        for (int c=j+1; c<=cols; c++)
            res(r,c-1)=(*this)(r,c);
    };
    for (int r=i+1; r<=rows; r++) {
        for (int c=1; c<j; c++)
            res(r-1,c)=(*this)(r,c);
        for (int c=j+1; c<=cols; c++)
            res(r-1,c-1)=(*this)(r,c);
    };
    return res;
};

double SquareMatrix::det() {
    double res=0;
    if (cols==1) return (*this)(1,1);
    else for (int i=1; i<=cols; i++)
            res+=NP(i)*(*this)(1,i)*((SquareMatrix)(this->restrict(1,i))).det();
    return res;
};

SquareMatrix SquareMatrix::adj() {
    SquareMatrix result(cols);
    for (int i=1; i<=rows; i++)
        for (int j=1; j<=cols; j++)
            result(j, i)=NP(i+j+1)*((SquareMatrix)restrict(i, j)).det();
    return result;
};

SquareMatrix SquareMatrix::inv() {
    SquareMatrix result(rows);
    double d=det();
    assert(d!=0);
    result=(1/d)*adj();
    return result;
};