Matrix Computation

#include<bits/stdc++.h>
#define ll long long

// Support:
//   Matrix multiplication
//   Solving system of equations
template <typename T> struct Matrix 
{ 
    vector< vector<T> > data; 
    int nrow, ncol;
    T MMOD;

    Matrix (int nrow_, int ncol_, T MMOD_=-1) : nrow(nrow_), ncol(ncol_), MMOD(MMOD_) {

        data.clear();
        data.resize(nrow, vector<T>(0, 0));
        for (int i = 0; i < nrow; ++i) {
            data[i].resize(ncol, 0);
        }
    }
    
    Matrix (int nrow_, int ncol_, vector< vector<T> > data_, T MMOD_=-1) : nrow(nrow_), ncol(ncol_), data(data_), MMOD(MMOD_) {}
    
    Matrix<T> (const Matrix<T> &mt) {
        nrow = mt.nrow;
        ncol = mt.ncol;
        data = mt.data;
        MMOD = mt.MMOD;
    }
    
    static Matrix identityMatrix(int nrow, int ncol, T MMOD_) {
        Matrix<T> m (nrow, ncol, MMOD_);
        for (int i = 0; i < nrow; ++i)
            m.data[i][i] = 1;
        return m;
    }
    
    Matrix operator * (Matrix b) const {
        
        if (ncol != b.nrow) {
            cout << "multiplication error: matrices' size not valid \n";
            return Matrix(0, 0, MMOD);
        }
      
        Matrix<T> res(nrow, b.ncol, MMOD);
        for (int i = 0; i < nrow; ++i) 
            for (int j = 0; j < b.ncol; ++j)
                for (int k = 0; k < ncol; ++k) {
                    
                    res.data[i][j] += data[i][k]*b.data[k][j];
                
                    if (MMOD > 0 && (res.data[i][j] > MMOD || res.data[i][j] < 0))
                        res.data[i][j] = ((res.data[i][j] % MMOD) + MMOD) % MMOD;
                }
        
        return res;
    }
    
    Matrix pow(ll exp) {
        
        if (ncol != nrow) {
            cout << "Power error: matrices' size not valid \n";
            return Matrix(0, 0, MMOD);
        }
        
        if (exp == 0) 
            return identityMatrix(nrow, ncol, MMOD);
        
        if (exp == 1)
//          return (*this);
            return Matrix<T> (nrow, ncol, data, MMOD);
            
        Matrix res = pow(exp >> 1);
        res = res * res;
        if (exp & 1)
            res = res * (*this);
//      cout << "exp: " << exp << endl;
//      res.print();
        return res;
    }
    
    vector<double> solveSystemEquations() {
        // Use Gauss-Elimination
        // Intend to solve system of linear equations of form:
        //
        //      a[0][0]*x[0]   + a[0][1]*x[1]   + a[0][2]*x[2]   + ... + a[0][n-1]*x[n-1]   = a[0][n]
        //      a[1][0]*x[0]   + a[1][1]*x[1]   + a[1][2]*x[2]   + ... + a[1][n-1]*x[n-1]   = a[1][n]
        //                                              ...
        //      a[n-1][0]*x[0] + a[n-1][1]*x[1] + a[n-1][2]*x[2] + ... + a[n-1][n-1]*x[n-1] = a[n-1][n]
        //
        // (n rows and n+1 columns)
        
        // pre-check
        T tmp = 1.1;
        if (tmp == 1) {
            cout << "Matrix is of type int (or long long) which is not suitable for gauss-elimination \n";
        }
        
        auto gdata = data;
                
        // do gauss-elimination
        for (int i = 0; i < nrow; ++i) {
            int pivot = i;
            for (int j = i+1; j < nrow; ++j)
                if (fabs(gdata[j][i]) > fabs(gdata[pivot][i]))
                    pivot = j;
            
            if (pivot != i)
                swap(gdata[i], gdata[pivot]);
            
            if (fabs(gdata[i][i]) < EPS)
                continue;
            
            for (int k = ncol-1; k >= i; k--)
                gdata[i][k] /= gdata[i][i];
            
            for (int j = i+1; j < nrow; ++j) {
                for (int k = ncol-1; k >= i; --k)
                    gdata[j][k] -= gdata[j][i] * gdata[i][k];
            }
        }
        
        // find x[0..n-1]
        for (int i = nrow-2; i >= 0; --i) {
            if (fabs(gdata[i][i]) > EPS)  // means: if gdata[i][i] == 1
                for (int k = i+1; k < ncol-1; ++k) 
                    if (isfinite(gdata[i][k]))
                        gdata[i][ncol-1] -= gdata[i][k] * gdata[k][ncol-1];
            else {
                gdata[i][ncol-1] = NAN;
            }
        }
        
        vector<double> x;
        for (int i = 0; i < nrow; ++i) {
            x.push_back(gdata[i][ncol-1]);
        }

        return x;
        
    }
    
    void print() {
        cout << "print matrix " << nrow << " x " << ncol << endl;
        for (int i = 0; i < nrow; ++i) {
            for (int j = 0; j < ncol; ++j)
                cout << data[i][j] << " ";
            cout << "\n";
        }
    }
  
};