矩阵乘法与斐波那契

矩阵乘法

A*B，只有当A的列数等于B的行数时才有意义，即A为$n*m$的矩阵，B为$m*L$的矩阵。矩阵的乘法定义如下：

$$\begin{bmatrix}a&b&c\\d&e&f\\\end{bmatrix}*\begin{bmatrix}g&h\\i&j\\k&l\\\end{bmatrix}=\begin{bmatrix}ag+bi+ck&ah+bj+cl\\dg+ei+fk&dh+ej+fl\\\end{bmatrix}$$
即 $$C_{i,j}=\sum_{k=1}^mA[i][k]+B[k][j]$$
由定义可知，矩阵乘法满足结合律，但不满足交换律。

应用

矩阵乘法可以用来优化很多线性递推式，下面以求斐波那契为例。
我们知道，$f[n]=f[n-1]+f[n-2]$，可以构造一个矩阵A，使得A*$\begin{bmatrix}f[n]\\f[n-1]\end{bmatrix}=\begin{bmatrix}f[n]+f[n-1]\\f[n]\end{bmatrix}$，这样，初始矩阵$\begin{bmatrix}1\\0\end{bmatrix}$每乘一次A数组，矩阵的第一位就往下推了一个，那么，要求第n项只要乘n-1次，而由于矩阵的乘法是满足结合律的，所以可以用快速幂来求矩阵A的n-1次幂，复杂度大大降低。
那么，如何设计矩阵A呢？
显然，A是一个2*2的矩阵，设$A=\begin{bmatrix}a&b\\c&d\end{bmatrix}$则$A*\begin{bmatrix}f[n]\\f[n-1]\end{bmatrix}=\begin{bmatrix}a*f[n]+b*f[n-1]\\c*f[n]+d*f[n-1]\end{bmatrix}$，显然，a=1,b=1,c=1,d=0，这样，我们就构造出

$$A=\begin{bmatrix}1&1\\1&0\end{bmatrix}$$
剩下的事就交给快速幂了。
下面附上poj3070的代码：

#include<cstdio>
#include<cstring>
#include<algorithm>
#define maxn 5
#define tt 10000
using namespace std;
inline char nc(){
    static char buf[100000],*i=buf,*j=buf;
    return i==j&&(j=(i=buf)+fread(buf,1,100000,stdin),i==j)?EOF:*i++;
}
inline int _read(){
    char ch=nc();int sum=0,p=1;
    while(ch!='-'&&!(ch>='0'&&ch<='9'))ch=nc();
    if(ch=='-')p=-1,ch=nc();
    while(ch>='0'&&ch<='9')sum=sum*10+ch-48,ch=nc();
    return sum*p;
}
struct matrix{
    int n,m,a[maxn][maxn];
    matrix operator *(matrix&b){
        matrix c;c.n=n;c.m=b.m;
        for(int i=1;i<=c.n;i++)
         for(int j=1;j<=c.m;j++){
            c.a[i][j]=0;
            for(int k=1;k<=m;k++) (c.a[i][j]+=a[i][k]*b.a[k][j])%=tt;
         }
        return c;
    }
};
int n;
matrix G,ans;
matrix power(matrix x,int y){
    if(y==1)return x;
    matrix c=power(x,y>>1);
    if(y&1)return c*c*x;
      else return c*c;
}
int main(){
    freopen("fib.in","r",stdin);
    freopen("fib.out","w",stdout);
    n=_read();
    while(n!=-1){
        G.n=G.m=2;G.a[1][1]=G.a[1][2]=G.a[2][1]=1;G.a[2][2]=0;
        ans.n=2;ans.m=1;ans.a[1][1]=1;ans.a[2][1]=0;
        if(n==0){
            printf("0\n");n=_read();
            continue;
        }
        if(n==1){
            printf("1\n");n=_read();
            continue;
        }
        G=power(G,n-1);
        ans=G*ans;
        printf("%d\n",ans.a[1][1]);
        n=_read();
    }
    return 0;
}