矩阵运算 快速幂

1. 矩阵乘法必须满足第一个矩阵的列等于第二个矩阵的行。
2. 结果矩阵的行是第一个矩阵的行，列是第二个矩阵的列。
3. 矩阵乘法行列各元素相乘。
4. 矩阵快速幂利用幂的二分优化。
5. 利用运算符重载进行一系列更方便的操作。

// 整数快速幂取模
template <typename Type>
long long ipow(Type x, Type n, Type mod){
    long long result = 1LL;
    while(n) {
        if(n & 1)
            result = (result*x) % mod;
        x = ((long long)x*x) % mod;
        n >>= 1;
    }

    return result;
}

//矩阵运算类
struct mat{
    int r = 0;
    int c = 0;
    vector<vector<long long> >rect;
    //construct
    mat(int r, int c) : r(r),c(c){
        rect.resize(r);
        for(int i = 0; i < r; ++i)
            rect[i].resize(c);
    }
    mat() = default;
    //input and output
    void in(){
        for(int i = 0; i < r; ++i)
            for(int j = 0; j < c; ++j)
                scanf("%I64d",&rect[i][j]);
    }
    void out(){
        for(int i = 0; i < r; ++i){
            for(int j = 0; j < c; ++j){
                if(j == 0)
                    printf("%I64d",rect[i][0]);
                else
                    printf(" %I64d",rect[i][j]);
            }
            printf("\n");
        }

    }
    //init
    void init(){
        for(int i = 0; i < r; ++i)
            for(int j = 0; j < c; ++j)
                rect[i][j] = (i == j);
    }
    //overload
    mat operator * (mat& rhs){
        mat ans(r,rhs.c);
        long long res;
        if(c == rhs.r){
            for(int i = 0; i < r; ++i){
                for(int j = 0; j < rhs.c; ++j){
                    res = 0LL;
                    for(int k = 0; k < c; ++k){//k means c or rhs.r
                        res += rect[i][k] * rhs.rect[k][j];
                    }
                    ans.rect[i][j] = res;
                }
            }
        }
        return ans;
    }
    void operator *= (mat& rhs){
        mat ans(r,rhs.c);
        long long res;
        if(c == rhs.r){
            for(int i = 0; i < r; ++i){
                for(int j = 0; j < rhs.c; ++j){
                    res = 0LL;
                    for(int k = 0; k < c; ++k){//k means c or rhs.r
                        res += rect[i][k] * rhs.rect[k][j];
                    }
                    ans.rect[i][j] = res;
                }
            }
        }
        rect = ans.rect;
    }
    mat operator ^ (int k){//r == c
        mat ans(r,c);
        ans.init();
        for(;k;k >>= 1){
            if(k & 1)
                ans *= (*this);
            (*this) *= (*this);
        }
        return ans;
    }
    void operator ^= (int k){//r == c
        mat ans(r,c);
        ans.init();
        for(;k;k >>= 1){
            if(k & 1)
                ans *= (*this);
            (*this) *= (*this);
        }
        (*this) = ans;
    }
};