很早之前写的 直接贴过来了

http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1113

矩阵快速幂：首先前置技能:  快速幂 + 矩阵乘法。

1  快速幂
1.1 快速乘法

题目：http://newoj.acmclub.cn/problems/2088 或者 acwing上的64位乘法

1.1.1

引用自2009年国家集训队论文，骆可强：《论程序底层优化的一些方法与技巧》 （膜膜膜）

可以根据需要换成uLL   (unsigned long long) 

// O(1)

#include <stdio.h>
typedef long long LL;
inline LL mult_mod(LL a, LL b, LL m){ // a * b % m
    a = (a % m + m) % m;
    b = (b % m + m) % m;
    LL ans = a * b - (LL)((long double) a * b / m + 0.5) * m;
    return ans < 0? ans + m : ans;
}

int main(){
    LL x, y, m;
    while(~scanf("%lld%lld%lld", &x, &y, &m)){
        printf("%lld\n", mult_mod(x, y, m));
    }
    return 0;
}

1.1.2  O(log b)

原理：

由于计算机底层设计的原因，做加法往往比乘法快的多，因此将乘法转换为加法计算将会大大提高（大数，比较小的数也没必要）乘法运算的速度，除此之外，当我们计算a*b%mod的时候，往往较大的数计算a＊b会超出long long int的范围，这个时候使用快速乘法方法也能解决上述问题． 
　　快速乘法的原理就是利用乘法分配率来将a＊b转化为多个式子相加的形式求解（注意这时使用乘法分配率的时候后面的一个乘数转化为二进制的形式计算）．举个栗子 
　　20＊14　＝　20＊(1110)2 = 20＊(2^3)＊1 + 20＊(2^2)＊1+20＊(2^1)＊１+20＊(2^0)＊0 = 160+80+40=280. 
　　这样计算了4次 而直接算的话 需要14次计算哦

typedef long long LL;
LL multi(LL a, LL b, LL m){
    // a * b % m
    LL ret = 0;
    while(b > 0){
        if(b & 1){
            ret = (ret + a) % m;
        }
        b >>= 1;
        a = (a << 1) % m;
    }
    return ret;
}

关于用哪个 emmm 大家自己决定吧， 一般都是用log b那个 

1.2  快速幂

其实快速幂的思想和log那个快速乘的思想差不多 也是将b转为二进制  不多赘述了

丢上百度链接 ： https://baike.baidu.com/item/快速幂/5500243?fr=aladdin

用到快速乘是因为 当a 和 b 太大的时候 用快速乘优化

#include <stdio.h>
typedef long long LL;
LL multi(LL a, LL b, LL m){
    // a * b % m
    LL ret = 0;
    while(b > 0){
        if(b & 1){
            ret = (ret + a) % m;
        }
        b >>= 1;
        a = (a << 1) % m;
    }
    return ret;
}

LL quick_mod(LL a, LL b, LL m){
    LL ans = 1;
    while(b){
        if(b & 1){
            ans = multi(ans,a,m);
            b--;
        }
        b /= 2;
        a = multi(a, a, m);
    }
    return ans;
}

2 矩阵乘法

2.1矩阵乘法

矩阵属于线性代数的东西。 受教于MR.ZZZ， 嘤嘤嘤 赵哥教了我好多东西。qaq

2.1.1 矩阵

由 m × n 个数aij排成的m行n列的数表称为m行n列的矩阵，简称m × n矩阵。记作：

2.1.2 矩阵乘法

前提： 矩阵A 乘 矩阵B  可以想乘的前提是 A的列数 == B的行数

两个矩阵的乘法仅当第一个矩阵A的列数和另一个矩阵B的行数相等时才能定义

矩阵A (m * a)  * 矩阵B(a * n) = 矩阵C(m * n) 

矩阵乘法如何用代码实现呢？

题目链接：http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1137

#include <stdio.h>
const int maxn = 1000 + 10;
int a[maxn][maxn], b[maxn][maxn];
long long c[maxn][maxn];
int main(){
    int n;
    scanf("%d", &n);
    for(int i = 1; i <= n; ++i){
        for(int j = 1; j <= n; ++j){
            scanf("%d", &a[i][j]);
        }
    }
    for(int i = 1; i <= n; ++i){
        for(int j = 1; j <= n; ++j){
            scanf("%d", &b[i][j]);
        }
    }
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++){
            for(int k = 1; k <= n; k++){
                c[i][j] += a[i][k] * b[k][j];
            }
            if(j == n){
                printf("%lld\n", c[i][j]);
            }
            else{
                printf("%lld ", c[i][j]);
            }
        }
    }
    return 0;
}

2.2 单位矩阵

在矩阵的乘法中，有一种矩阵起着特殊的作用，如同数的乘法中的1，这种矩阵被称为单位矩阵。它是个方阵，从左上角到右下角的对角线（称为主对角线）上的元素均为1。除此以外全都为0。

单位矩阵的乘法： 矩阵A * 单位矩阵 = 矩阵A

3 矩阵快速幂
3.1 有了前两个前置技能后， 那么矩阵快速幂就很好理解了。 我们可以把矩阵看作一个特殊的数字，然后去用快速幂的思想去写，最开始的ans为单位矩阵(单位矩阵的特性)。我们重写一下乘法就好了，具体看代码吧。

题目：http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1113

//
//  1113 - 矩阵快速幂.cpp
//  数论
//
//  Created by Terry on 2018/9/27.
//  Copyright © 2018年 Terry. All rights reserved.
//
//  1: 单位矩阵 * 矩阵A = 矩阵A
//  2: 矩阵快速幂 可以对应 快速幂去理解 就是重写一下乘法  a
//      把快速幂里面的数字乘法 重写成 矩阵乘法
#include <stdio.h>
int read(){
    int x = 0, f = 1;
    char ch=getchar();
    while(ch < '0' || ch > '9'){if(ch=='-') f = -1; ch = getchar();}
    while(ch >= '0' && ch <= '9'){x = x * 10 + ch - '0'; ch = getchar();}
    return x * f;
}
typedef long long LL;
const int mod = 1e9 + 7;
const int maxn = 100 + 10;
struct Matrix{
    int m[maxn][maxn];
}unit;
int n;
Matrix operator * (Matrix a, Matrix b){
    Matrix ret;
    LL x;
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            x = 0;
            for(int k = 0; k < n; ++k){
                x += ((LL)a.m[i][k] * b.m[k][j]) % mod;
            }
            ret.m[i][j] = x % mod;
        }
    }
    return ret;
}
void init_unit(){
    // 单位矩阵
    for(int i  = 0; i < maxn; i++){
        unit.m[i][i] = 1;
    }
}
Matrix pow_mat(Matrix a, LL n){
    Matrix ans = unit;
    while (n) {
        if(n & 1){
            ans = ans * a;
        }
        a = a * a;
        n >>= 1;
    }
    return ans;
}
int main(){
    init_unit();
    n = read();
    LL x = read();
    Matrix a;
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            a.m[i][j] = read();
        }
    }
    a = pow_mat(a, x);
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            if(j == n - 1){
                printf("%d\n", a.m[i][j]);
            }
            else{
                printf("%d ", a.m[i][j]);
            }
        }
    }
    return 0;
}