矩阵

template<class T>
struct Matrix {
    int N, M;
    vector<vector<T>> data;
    Matrix(int n, int m) : N(n), M(m) {
        data.resize(n);
        for (int i = 0 ; i < n ; i++) {
            data[i].resize(m);
        }
    }
    Matrix E() {
        Matrix<T> mat(N, M);
        for (int i = 0 ; i < N ; i++) {
            for (int j = 0 ; j < M ; j++) {
                if (i == j) mat[i][j] = T(1);
                 else mat[i][j] = T(0);
            }
        }
        return mat;
    }
    void clear() {
        for (int i = 0 ; i < N ; i++) {
            data[i].resize(0);
            data[i].resize(M);
        }
    }
    static Matrix qpow(Matrix base, long long b) {
        assert(base.N == base.M);
        Matrix res(base.N, base.M);
        while (b) {
            if (b & 1) res *= base;
            base *= base;
            b >>= 1;
        }
        return res;
    }
    Matrix& operator*=(const Matrix &x) {
        assert(M == x.N);
        Matrix res(N, x.M);
        for (int i = 0 ; i < N ; i++) {
            for (int k = 0 ; k < M ; k++) {
                T r = data[i][k];
                for (int j = 0 ; j < x.M ; j++) {
                    res[i][j] += r * x.data[k][j];
                }
            }
        }
        return *this = res;
    }
    Matrix operator*(const Matrix &x) {
        Matrix res = *this;
        return res *= x;
    }
    vector<T>& operator[](int x) {
        return data[x];
    }
    friend istream& operator>>(istream& is, Matrix &mat) {
        for (int i = 0 ; i < mat.N ; i++) {
            for (int j = 0 ; j < mat.M ; j++) {
                is >> mat.data[i][j];
            }
        }
        return is;
    }
    friend ostream& operator<<(ostream& os,  Matrix &mat) {
        for (int i = 0 ; i < mat.N ; i++) {
            for (int j = 0 ; j < mat.M ; j++) {
                os << mat.data[i][j] << " ";
            }
            os << "\n";
        }
        return os;
    }
}; // Matrix
// Matrix 矩阵类
// 行列下标从 0 开始
// 实现了矩阵的乘法运算和快速幂
// 实现了矩阵的输入和输出