矩阵快速幂的实现

//
#include <algorithm>
#include <iostream>
#include <utility>
#include <sstream>
#include <cstring>
#include <cstdio>
#include <vector>
#include <queue>
#include <stack>
#include <cmath>
#include <map>
#include <set>
using namespace std;

typedef long long ll;
const int MOD = 542;
//int MOD = 99990001;
const int INF = 0x3f3f3f3f;
const ll INFF = 0x3f3f3f3f3f3f3f3fLL;
const double EPS = 1e-9;
const double OO = 1e20;
const double PI = acos(-1.0); //M_PI;
const int fx[] = {-1, 1, 0, 0};
const int fy[] = {0, 0, -1, 1};
const int maxn=50 + 5;
struct matrix{
        int mat[maxn][maxn];
};
//这里只适用于n*n的矩阵，当然要转换成其他形式的矩阵也是相当轻松的
matrix matrix_multiply(matrix a,matrix b,int n){
        matrix c;
        memset(c.mat,0,sizeof(c.mat));
        //c[i][j]= Σ a[i][k]*b[k][j] (0<=k<=n)
//        for(int i=0;i<n;i++)
//                for(int j=0;j<n;j++)
//                        for(int k=0;k<n;k++)
//                                if(a.mat[i][k] && b.mat[k][j]) c.mat[i][j]=(c.mat[i][j] + (a.mat[i][k]*b.mat[k][j])%MOD)%MOD;
//上面的写法虽然是对的，但是可以进一步的优化，由于矩阵a中的某一个元素那么后面关于这个元素的所有的运算都没有必要进行了
        for(int i=0;i<n;i++)
        for(int k=0;k<n;k++){
                if(a.mat[i][k]!=0){
                        for(int j=0;j<n;j++){
                                if(b.mat[k][j]==0) continue;
                                c.mat[i][j]=(c.mat[i][j] + (a.mat[i][k]*b.mat[k][j])%MOD)%MOD;
                        }
                }
        }
        return c;
}
matrix matrix_index(matrix a,int mi,int n){
        matrix t;
        for(int i=0;i<n;i++)
                for(int j=0;j<n;j++)
                        t.mat[i][j]=(i==j);
        while(mi){
                if(mi&1) t=matrix_multiply(t,a,n);
                a=matrix_multiply(a,a,n);
                mi>>=1;
        }
        return t;
}
int main(){
        int n,N;cin >> n >> N;
        matrix a,c;
        memset(a.mat,0,sizeof(a.mat));memset(c.mat,0,sizeof(c.mat));
        for(int i=0;i<n;i++)
                for(int j=0;j<n;j++)
                        scanf("%d",&a.mat[i][j]);
        c=matrix_index(a,N,n);
        for(int i=0;i<n;i++){
                for(int j=0;j<n;j++)
                        cout << c.mat[i][j] << " ";
                cout << endl;
        }
        return 0;
}