矩阵乘法总结

矩阵乘法

矩阵乘法是用于优化一些递推式的方法。

定义

两个矩阵的乘法当且仅当第一个矩阵A的列数和第二个矩阵B的行数相等下才有定义。
假设A是$n\times m$ 的矩阵，B是$m\times p$ 的矩阵，那么他们的乘积C一定是一个$n\times p$ 的矩阵。其中C任意的一个元素值为：

$$c_{i,j}=a_{i,1}b_{1,j}+a_{i,2}b_{2,j}+a_{i,3}b_{3,j}+...+a_{i,m}b_{m,j}=\sum _{r=1}^{m}a_{i,r}b_{r,j}$$
我们将这个乘积记为: $C=A\times B$。
显然乘法法则满足结合律，各种分配律，但是不满足交换律（显然叉出来的矩阵都不一样）。即：
$$(AB)C=A(BC)$$ $$(A+B)C=AC+BC$$ $$C(A+B)=CA+CB$$

应用

比如运用矩阵乘法我们可以优化我们所熟知的斐波那契数列。
斐波那契数列的递推式为：

$$f(i)=f(i-1)+f(i-2)$$
可以看成矩阵 $$\begin{bmatrix}f(i-2)&f(i-1)\end{bmatrix}$$
我们只需要构造一个矩阵B，使得 $$\begin{bmatrix}f(i-2)&f(i-1)\end{bmatrix}\times B=\begin{bmatrix}f(i-1)&f(i)\end{bmatrix}$$
不难推出 $B=\begin{bmatrix}0&1\\1&1\end{bmatrix}$
正确性读者可以自行验证。

这样我们只需要从初始状态$\begin{bmatrix}0&1\end{bmatrix}\times B^{n-1}$就可以推出$n$项了。
$B^{n-1}$显然可以用快速幂求出。

看到这里肯定有一个疑问，初值”1”怎么给，就是使得$A\times B=A$。
这里引入单位矩阵的概念，其实就是上述问题的答案，显然单位矩阵是一个$m\times m$的矩阵。

用$\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}$表示单位矩阵。

例题

poj 3070,求斐波那契数列的第n项（n的范围在以前的我眼中是不正常的QAQ）

#include<cstdio>
#include<cstring>
using namespace std;
const int tt=10000;
int n;
struct jz{
    int r,c,s[2][2];
    void clean(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 work(){clean(2,2);for (int i=0;i<c;i++) s[i][i]=1;}
};
jz operator*(const jz &a,const jz &b){
        jz c;c.clean(2,2);
        for (int i=0;i<a.r;i++)
        for (int j=0;j<b.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])%tt;
        return c;
    }
jz a,ans;
jz qsm(jz w,int b){
    jz num;num.work();
    while(b){
        if (b%2==1) num=num*w;
        b=b>>1;if (b) w=w*w;
    }
    return num;
}
int main(){
    freopen("exam.in","r",stdin);
    freopen("exam.out","w",stdout);
    a.r=a.c=2;a.s[0][0]=0;a.s[0][1]=a.s[1][0]=a.s[1][1]=1;
    ans.r=1,ans.c=2;ans.s[0][0]=0;ans.s[0][1]=1;
    while(1){
        scanf("%d",&n);
        if (n==-1) return 0;
        if (!n){printf("0\n");continue;}
        printf("%d\n",(ans*qsm(a,n-1)).s[0][1]);
    }
}