矩阵快速幂

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define md  1000000007
#define maxn 10
ll ans,b,n;
void mulp(ll x[9][9],ll y[9][9])
{
    ll z[maxn][maxn];
    memset(z,0,sizeof(z));
    for(int k=0; k<9; k++)
        for(int i=0; i<9; i++)
            for(int j=0; j<9; j++)
                z[i][j]=(z[i][j]+x[i][k]*y[k][j]%md)%md;
    for(int i=0; i<9; i++)
        for(int j=0; j<9; j++)
            x[i][j]=z[i][j];
}
int main()
{
    cin>>n;
    if(!n)
    {
        cout<<1<<endl;
        return 0;
    }
    b=n/5;
    ll a[9][9]= {0,16,0,9,0,4,0,1,0,
                 1,0,0,0,0,0,0,0,0,
                 0,1,0,0,0,0,0,0,0,
                 0,0,1,0,0,0,0,0,0,
                 0,0,0,1,0,0,0,0,0,
                 0,0,0,0,1,0,0,0,0,
                 0,0,0,0,0,1,0,0,0,
                 0,0,0,0,0,0,1,0,0,
                 0,0,0,0,0,0,0,1,0
                };
    ll e[9][9]= {1,0,0,0,0,0,0,0,0,
                 0,1,0,0,0,0,0,0,0,
                 0,0,1,0,0,0,0,0,0,
                 0,0,0,1,0,0,0,0,0,
                 0,0,0,0,1,0,0,0,0,
                 0,0,0,0,0,1,0,0,0,
                 0,0,0,0,0,0,1,0,0,
                 0,0,0,0,0,0,0,1,0,
                 0,0,0,0,0,0,0,0,1
                };
    while(b)
    {
        if(b%2)
            mulp(e,a);
        mulp(a,a);
        b/=2;
    }
    for(int i=0; i<9; i++)
        ans=(ans+e[i][0])%md;
    cout<<ans<<endl;
    return 0;
}

#include<iostream>
#include<iomanip>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#include<string>
#include<queue>
#include<stack>
#include<vector>
#include<map>
#include<set>
#include<cstdlib>
#define rep(i,a,n) for(int i=a;i<=n;i++)
#define ll long long
#define ull unsigned long long
#define mem(a,x) memset((a),(x),sizeof ((a)))//x只能是0或-1或false或true
#define debug(x) cout<<"X: "<<(x)<<endl
#define de cout<<"************"<<endl
#define lowbit(x) ((x)&(-x))
#define lson rt<<1
#define rson rt<<1|1
#define gcd(a,b) __gcd(a,b)
#define lcm(a,b) a*b/(__gcd(a,b))
#define inf 0x3f3f3f3f//1e9+6e7
#define eps 1e-8
#define mod 1e9+7
#define N 1000010
const double pi=acos(-1.0);
using namespace std;
int a[110][110],b[110][110];
int z[110][110];
void mulp(int x[][110],int y[][110],int n,int m,int p,int md)
{
    //n*m的矩阵，m*p的矩阵，两个矩阵相乘，模md
    //结果存在全局变量z中，z为一个n*p的矩阵
    memset(z,0,sizeof(z));
    for(int k=0; k<m; k++)
        for(int i=0; i<n; i++)
            for(int j=0; j<p; j++)
                z[i][j]=(z[i][j]+x[i][k]*y[k][j]%md)%md;
}
int main()
{
    int n,m,p;
    cin>>n>>m>>p;
    for(int i=0; i<n; i++)
        for(int j=0; j<m; j++)
            cin>>a[i][j];
    for(int i=0; i<m; i++)
        for(int j=0; j<p; j++)
            cin>>b[i][j];
    mulp(a,b,n,m,p,mod);
    for(int i=0; i<n; i++)
    {
        for(int j=0; j<p; j++)
            cout<<z[i][j]<<" ";
        cout<<endl;
    }
    return 0;
}
/*
2 3 2
1 2 3
4 6 7
5 2
1 4
3 1
output：
16 13
47 39
*/