定义：

设A为 的矩阵，B为 的矩阵，那么称 的矩阵C为矩阵A与B的乘积，记作 ，其中矩阵C中的第行第列元素可以表示为：

如下所示：
性质：
满足结合律、分配率，不满足交换律
应用：
可以用来加速形如dp[i]=a*dp[i-1]+b*dp[i-2]+...+c这样的dp（我只知道这个用途。。）
举个例子，f[i]=p*f[i-1]+q*f[i-2]，且第一项是a1，第二项是a2，求f[n]。那么f[i]和f[i-1]可以由f[i-1]和f[i-2]得到，就可以把它们分别放到一个长度为2的列向量B、C中，那么用一个2*2大小的存储系数的矩阵A乘上B就可以等于C（注意是n*n的矩阵乘上n*1的列向量得到一个n*1的列向量，由于不满足交换律所以不能写反）。那A怎么求呢？想到C[i]=A[i][1]*B[1]+A[i][2]*B[2]+...+A[i][n]*B[n]，就可以根据这个公式，直接在矩阵中填上就可以了。那么例子中就是：
(p,q)    (f[i-1])    (f[i])
(1,0) * (f[i-2]) = (f[i-1])
由于已经知道了前两项，只要做n-2次矩阵乘法就可以了，设状态数是k，那么时间复杂度是O(n*k^3)???
然而由于矩阵乘法满足交换律，即A*(A*(A*(...*B)))=(A*A*...*A)*B，那么显然可以用快速幂优化，于是时间复杂度就可以降为O(log n*k^3)，只要状态数不太大就好了。
附矩阵乘法+快速幂代码：
//定义一个矩阵,表示一个n*m大小的矩阵，根据题目需要定义maxk
struct mat
{
    ll a[maxk][maxk],n,m;
}x,y;
//矩阵乘法
mat mul(mat x,mat y)
{
    int i,j,k;
    mat ans;
    for(i=0;i<x.n;i++)
        for(j=0;j<y.m;j++)
            ans.a[i][j]=0;//初始化矩阵
    for(k=0;k<x.m;k++)
        for(i=0;i<x.n;i++)
            for(j=0;j<y.m;j++)
                ans.a[i][j]=(ans.a[i][j]+x.a[i][k]*y.a[k][j])%m;
    ans.n=x.n;ans.m=y.m;
    //注意m*p的矩阵乘上p*n的矩阵等于m*n的矩阵，做矩阵乘法时不要弄错
    return ans;
}
//快速幂
mat pow(mat x,ll y)
{
    if(y==1)return x;
    mat ans=pow(x,y/2);
    ans=mul(ans,ans);
    if(y&1)ans=mul(ans,x);
    return ans;
}