和快速幂一模一样，只是把数乘换成矩阵乘，只需要定义结构体矩阵然后重载一下乘法*即可。
注意：
$1$乘以任何数都等于这个数本身
单位矩阵乘以任何矩阵就等于这个矩阵本身

#define IO ios::sync_with_stdio(false);cin.tie();cout.tie(0)
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const ll mod=1e9+7;
const int N=110;
int n;
ll k;
struct node
{ll m[N][N];node(){memset(m,0,sizeof m);};node operator *(const node &b) const{node res;for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)for(int k=1;k<=n;k++)   res.m[i][j]=(res.m[i][j]+m[i][k]*b.m[k][j])%mod;return res;}
};
node qmi(node a,ll b)
{node res;for(int i=1;i<=n;i++)// 单位矩阵res.m[i][i]=1;while(b){if(b&1) res=res*a;a=a*a;b>>=1;}return res;
}
int main()
{IO;int T=1;//cin>>T;while(T--){node a;cin>>n>>k;for(int i=1;i<=n;i++)for(int j=1;j<=n;j++) cin>>a.m[i][j];node res=qmi(a,k);for(int i=1;i<=n;i++){for(int j=1;j<=n;j++) cout<<res.m[i][j]<<' ';cout<<'\n';}}return 0;
}

$$\left[\begin{array}{cc}0& 0\\ 1& 1\end{array}\right]\left[\begin{array}{c}{f}_{n-2}\\ f_{n-1}\end{array}\right]=\left[\begin{array}{c}{f}_{n-1}\\ f_n\end{array}\right]\to {\left[\begin{array}{cc}0& 0\\ 1& 1\end{array}\right]}^{n-2}\left[\begin{array}{c}{f}_1\\ f_2\end{array}\right]=\left[\begin{array}{c}{f}_{n-1}\\ f_n\end{array}\right]$$
矩阵乘法满足结合律，由此可以根据上述式子进行矩阵快速乘

#define IO ios::sync_with_stdio(false);cin.tie();cout.tie(0)
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const ll mod=1e9+7;
const int N=110;
int sz;//矩阵大小
ll n,k;
struct node
{ll m[N][N];node(){memset(m,0,sizeof m);};node operator *(const node &b) const{node res;for(int i=1;i<=sz;i++)for(int j=1;j<=sz;j++)for(int k=1;k<=sz;k++)   res.m[i][j]=(res.m[i][j]+m[i][k]*b.m[k][j])%mod;return res;}
};
node qmi(node a,ll b)
{node res;for(int i=1;i<=sz;i++)// 单位矩阵res.m[i][i]=1;while(b){if(b&1) res=res*a;a=a*a;b>>=1;}return res;
}
int main()
{IO;int T=1;//cin>>T;while(T--){cin>>n;if(n<=2) {cout<<1<<'\n';continue;}sz=2;node a;a.m[1][1]=0,a.m[1][2]=1,a.m[2][1]=1,a.m[2][2]=1;node res=qmi(a,n-2);ll ans=0;ans=(ans+res.m[1][2]+res.m[2][2])%mod;cout<<ans<<'\n';}return 0;
}

P1939 【模板】矩阵加速（数列）和上面这个题基本一样。

$$\left[\begin{array}{cc}a& 1\\ 0& 1\end{array}\right]\left[\begin{array}{c}{x}_{n-1}\\ c\end{array}\right]=\left[\begin{array}{c}{x}_n\\ c\end{array}\right]\to {\left[\begin{array}{cc}a& 1\\ 0& 1\end{array}\right]}^n\left[\begin{array}{c}{x}_0\\ c\end{array}\right]=\left[\begin{array}{c}{x}_n\\ c\end{array}\right]$$
这题比较dt的地方是两个数相乘会爆long long，于是上龟速乘防止乘积爆

#define IO ios::sync_with_stdio(false);cin.tie();cout.tie(0)
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=110;
int sz;//矩阵大小
ll mod;
ll mul(ll a,ll b)//龟速乘
{ll res=0;a%=mod;while(b){if(b&1) res=(res+a)%mod;a=(a+a)%mod;b>>=1;}return res;
}
struct node
{ll m[N][N];node(){memset(m,0,sizeof m);};node operator *(const node &b) const{node res;for(int i=1;i<=sz;i++)for(int j=1;j<=sz;j++)for(int k=1;k<=sz;k++)   res.m[i][j]=(res.m[i][j]+mul(m[i][k],b.m[k][j]))%mod;return res;}
};
node qmi(node a,ll b)
{node res;for(int i=1;i<=sz;i++)// 单位矩阵res.m[i][i]=1;while(b){if(b&1) res=res*a;a=a*a;b>>=1;}return res;
}
int main()
{IO;int T=1;//cin>>T;sz=2;ll a,c,x0,n,g;while(T--){cin>>mod>>a>>c>>x0>>n>>g;node now;now.m[1][2]=now.m[2][2]=1,now.m[1][1]=a;node res=qmi(now,n);ll ans=(mul(res.m[1][1],x0)+mul(res.m[1][2],c))%mod%g;cout<<ans<<'\n';}return 0;
}

矩阵乘法优化dp
考虑P1115 最大子段和如何做？
不难得知设计dp即可。
状态表示：$f_i$表示以第$i$个位置为结尾最大字段和，$g_i=max(f_1,f_2,\dots ,f_i)$
状态转移：$f_i=max(f_{i-1}+a_i,a_i)$，$g_i=max(g_{i-1},f_i)$
最终答案即是$g_n$

总所周知递推不满足结合律，换句话说就是你必须一步一步递推，不过矩阵乘法能够优化递推 斐波那契前 n 项和 ，而优化的方式就是使计算过程具有结合律，那么如果我们通过矩阵操作使一些不具有结合律的东西具有结合律那么我们就能用线段树维护这个东西，花费log代价显然使非常优秀的。

效仿矩阵乘法优化斐波那契的方法寻找矩阵
$$\left[\begin{array}{ccc}{a}_i& -\infty & a_i\\ a_i& 0& a_i\\ -\infty & -\infty & 0\end{array}\right]?\left[\begin{array}{c}{f}_{n-1}\\ g_{n-1}\\ 0\end{array}\right]=\left[\begin{array}{c}{f}_n\\ g_n\\ 0\end{array}\right]$$
$?$表示一个运算
$$\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]?\left[\begin{array}{c}e\\ f\end{array}\right]=\left[\begin{array}{c}max(a+e,b+f)\\ max(c+e,d+f)\end{array}\right]$$
不难发现上述定义的新运算是具有结合律的！！！
对比此运算和矩阵乘法与运算不难发现：
矩阵乘法中的“乘”相当于这里的“加”而矩阵乘法中的“加”相当于这里的“max”
矩阵乘法满足结合律实际上使“乘”对“加”满足分配了，而这里“加”对“max”同样满足分配率于是新运算具有结合律。

满足结合律并且区间查询单点修改无疑线段树，只需要每个节点维护一个矩阵即可。

考虑答案在哪?
定义：
$$A_i=\left[\begin{array}{ccc}{a}_i& -\infty & a_i\\ a_i& 0& a_i\\ -\infty & -\infty & 0\end{array}\right]$$
$$B_i=A_1?A_2?\dots ?A_i$$
那么再看此题 P1115 最大子段和，不难得出
$$B_n?\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]=\left[\begin{array}{c}{f}_n\\ g_n\\ 0\end{array}\right]$$
答案就是$g_n$，不难发现答案就蕴藏着$B_n$矩阵中，分析一下不难得知如果$B_n=\left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{12}\\ b_{21}& b_{22}& b_{23}\\ b_{31}& b_{32}& b_{33}\end{array}\right]$那么答案就是$max(b_{21},b_{23})$

而本题有了之前的工作只需要套个线段树就是基本的区间修改单点查询问题。

#define IO ios::sync_with_stdio(false);cin.tie();cout.tie(0)
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const int maxn=3,INF=0x3f3f3f3f;
const int N=50010;//投机取巧
int n,q,a;
int sz;//矩阵大小
struct node
{int l,r;int m[maxn][maxn];node(){l=r=0;memset(m,0,sizeof m);};node operator *(const node &b) const{node res;res.l=l,res.r=b.r;memset(res.m,-0x3f,sizeof res.m);for(int i=0;i<sz;i++)for(int j=0;j<sz;j++)for(int k=0;k<sz;k++)   res.m[i][j]=max(res.m[i][j],m[i][k]+b.m[k][j]);return res;}int ans(){return max(m[1][0],m[1][2]);}
}tree[N<<2];
void pushup(int u)
{tree[u]=tree[u<<1]*tree[u<<1|1];
}
void build(int u,int l,int r)
{if(l==r){cin>>a;//这样输入省空间tree[u].l=tree[u].r=r;tree[u].m[0][0]=tree[u].m[1][0]=tree[u].m[0][2]=tree[u].m[1][2]=a;tree[u].m[0][1]=tree[u].m[2][0]=tree[u].m[2][1]=-INF;return;}int mid=l+r>>1;build(u<<1,l,mid);build(u<<1|1,mid+1,r);pushup(u);
}
void modify(int u,int pos,int x)
{if(tree[u].l==tree[u].r){tree[u].m[0][0]=tree[u].m[1][0]=tree[u].m[0][2]=tree[u].m[1][2]=x;return;}int mid=tree[u].l+tree[u].r>>1;if(pos<=mid) modify(u<<1,pos,x);else modify(u<<1|1,pos,x);pushup(u);
}
node query(int u,int l,int r)
{if(tree[u].l>=l&&tree[u].r<=r) return tree[u];int mid=tree[u].l+tree[u].r>>1;if(r<=mid) return query(u<<1,l,r);else if(l>mid)return query(u<<1|1,l,r);else return query(u<<1,l,r)*query(u<<1|1,l,r);
}
int main()
{IO;cin>>n;sz=3;build(1,1,n);cin>>q;while(q--){int op,x,y;cin>>op>>x>>y;if(op==1){if(x>y) swap(x,y);cout<<query(1,x,y).ans()<<'\n';}elsemodify(1,x,y);}return 0;
}