矩阵乘法

定义

矩阵乘法是矩阵的一种运算。设A是a*b的矩阵，B是c*d的矩阵，只有当b等于c时，矩阵乘法才能进行。
满足结合律，分配律，但不一定满足交换律（b等于c时a并不一定等于d，所以不能交换）。

乘积

设C=A*B，则：
$C_{i,j}=\sum_{k=1}^{b}A_{i,k}*B_{k,j}$
比如：
$\begin{bmatrix}1&0&2\\-1&3&1\\\end{bmatrix}* \begin{bmatrix}3&1\\2&1\\1&0\\\end{bmatrix}= \begin{bmatrix}1*3+0*2+2*1&1*1+0*1+2*0\\-1*3+3*2+1*1&-1*1+3*1+1*0\\\end{bmatrix}= \begin{bmatrix}5&1\\4&2\end{bmatrix}$

作用

矩阵乘法在信息学中的应用通常是加快递推，只需要将递推表示成矩阵的形式，就可以利用矩阵乘法来做了。下面以最经典的斐波那契数来举例。

斐波那契数：$f(i)=f(i-1)+f(i-2)$
构造这样的状态矩阵：
$\begin{bmatrix}f(n-1)&f(n)\end{bmatrix}$
再构造这样的变换矩阵：
$\begin{bmatrix}a&b\\c&d\end{bmatrix}$
根据
$\begin{bmatrix}f(n-2)&f(n-1)\end{bmatrix}* \begin{bmatrix}a&b\\c&d\end{bmatrix}= \begin{bmatrix}a*f(n-2)+c*f(n-1)&b*f(n-2)+d*f(n-1)\end{bmatrix}= \begin{bmatrix}f(n-1)&f(n)\end{bmatrix}$
可以算出变换矩阵的系数：
$\begin{bmatrix}0&1\\1&1\\\end{bmatrix}$
那么也就是说$\begin{bmatrix}f(n-2)&f(n-1)\end{bmatrix}* \begin{bmatrix}0&1\\1&1\end{bmatrix}= \begin{bmatrix}f(n-1)&f(n)\end{bmatrix}$。而我们已知了f(0)和f(1)，变换一次就可以变成f(1)和f(2)，再变换一次就可以变成f(2)和f(3)，以此类推，我们只要变换n-1次就可以得到f(n-1)和f(n)了。

由于矩阵乘法满足结合律，我们可以这么求f(n-1)和f(n)：
$\begin{bmatrix}f(0)&f(1)\end{bmatrix}* \begin{bmatrix}0&1\\1&1\end{bmatrix}^{n-1}= \begin{bmatrix}f(n-1)&f(n)\end{bmatrix}$
既然是乘法，就可以用快速幂。那么只需要$O(2^3*log_2(n))$的时间就可以求出斐波那契第n项了。

不过矩阵的“1”是什么？矩阵中的“1”被称为单位矩阵，是形如这样的矩阵：
$\begin{bmatrix} 1&0&\cdots&0&0\\ 0&1&\cdots&0&0\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&1&0\\ 0&0&\cdots&0&1\\ \end{bmatrix}$
其中i,i为1，其余为0。
一个a*b的矩阵乘上b*b的单位矩阵显然会变成自己。

模板

POJ3070，求斐波那契第n项。

#include<cstdio>
#include<cstring>
using namespace std;
const int MOD=10000;

int n;
struct Matrix
{
    int r,c,s[2][2];
    void zero(int R,int C) //零矩阵
    {
        r=R;c=C;
        for (int i=0;i<r;i++)
        for (int j=0;j<c;j++)
            s[i][j]=0;
    }
    void unit(int C) {zero(C,C);for (int i=0;i<c;i++) s[i][i]=1;}
    //单位矩阵
};
Matrix fst,T,ans,c;
Matrix operator * (const Matrix &a,const Matrix &b) //矩阵乘法
{
    c.zero(a.r,b.c);
    for (int i=0;i<c.r;i++)
    for (int j=0;j<c.c;j++)
    for (int k=0;k<a.c;k++)
        c.s[i][j]=(c.s[i][j]+a.s[i][k]*b.s[k][j])%MOD;
    return c;
}

Matrix power(Matrix w,int b) //快速幂
{
    Matrix s;s.unit(w.c);
    while (b)
    {
        if (b&1) s=s*w;
        b>>=1;
        if (b) w=w*w;
    }
    return s;
}
int main()
{
    freopen("Matrix_mul.in","r",stdin);
    freopen("Matrix_mul.out","w",stdout);
    fst.r=1;fst.c=2;fst.s[0][0]=0;fst.s[0][1]=1;
    T.r=T.c=2;T.s[0][0]=0;T.s[0][1]=T.s[1][0]=T.s[1][1]=1;
    for (scanf("%d",&n);~n;scanf("%d",&n))
    {
        if (!n) printf("0\n"); else
        printf("%d\n",(fst*power(T,n-1)).s[0][1]);
    }
    return 0;
}