矩阵乘法

首先推荐matrix67的一篇好文章，《十个利用矩阵乘法解决的经典题目》。讲的非常清晰。

 为什么他讲矩阵讲的很清晰呢？因为他叫matrix67。。（大雾

---

 矩阵乘法：P_i,j=ΣA_i,k*B_k,j

显而易见 ，矩阵乘法满足结合律但不满足交换律，程序里不要写反了。

---

 矩阵快速幂：因为矩阵乘法满足结合律，所以可以用快速幂加速。复杂度O(n³logK)。

注意通常把矩阵的快速幂写成迭代而不是递归，因为矩阵通常比较占地方，写成迭代比较好。

Quickmatrix(martix m,int k){
    //注意ans的初始化，或者在第一次使用ans时令ans=m;
    while (k){
        if (k&1)    ans*=m;
        m=m*m;
        k>>=1;
    }
}

---

 常见应用1：

利用矩阵乘法加速线性递推。因为矩阵本就是描述线性变换的，所以用来表示线性递推就很自然。

例如计算斐波那契数列的第n项，就可以根据斐波那契数列的定义构造一个矩阵，使其能将向量（F_n-2,F_n-1）变换到向量（F_n-1,F_n-2+F_n-1）

显然，我们可以得到，那么第n项就是这个2 x 2的矩阵自乘n次，再乘以(0,1)的结果。

---

 常见应用2：

利用矩阵乘法统计图中的路径方案数。

根据定义，某个图的邻接01矩阵自乘n次后，A_i,j就是从i到j经过n条边的路径方案数。

这个问题可以理解为一般意义上的计数问题。很多统计方案的问题，实质上就是其状态图中的路径方案数问题。

例如文章中的例题9

---

【例题】SCOI2009 迷路

给出一个有向有权图，问从1到n，长度为T的路径一共有多少条。

注意到题目中每条路径的长度不超过10，那么就可以拆点转化为应用2.

 附上代码：

#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <iostream>
#include <cstring>
#include <cmath>
#include <stack>
#include <ctime>
#define INF (9999999)
#define MAXN (100+10)
#define BreakTest //cout<<"Here!"<<endl;
using namespace std;
/*Gloable*/
int map[MAXN][MAXN];
int newpoint[MAXN][11],np=1;
int inside[MAXN];
/*ADT & reload*/
struct matrix{
    int p[MAXN][MAXN];
    int x,y;
    matrix(){memset(p,0,sizeof(p));}
}m;
typedef struct matrix matrix;
matrix operator * (matrix a ,matrix b){//a.x==b.y
    matrix ans;
    for(int i=1;i<=a.x;i++)
        for (int j=1;j<=b.y;j++)
            for (int k=1;k<=a.y;k++)
                ans.p[i][j]=(((a.p[i][k]%2009)*(b.p[k][j]%2009))%2009+ans.p[i][j]%2009)%2009;
    ans.x=a.x;
    ans.y=b.y;
    return ans;
}
/*Funciton*/
matrix quickmatrix(matrix m,long long k){
    matrix ans;
    bool exits=false;
    while (k){
        if (k&1==1){
            if (!exits) {ans=m;exits=true;}
            else
            ans=ans*m;
        }
        k>>=1;
        m=m*m;
    }
    return ans;
}
int main()
{
    int n;
    long long t;
    char str[15];
    int i,j,k;
    scanf("%d%lld",&n,&t);
    for (i=1;i<=n;i++){
        scanf("%s",str);
        for (j=0;str[j];j++){
            if (str[j]!='0')
                map[i][j+1]=str[j]-'0';
            else
                map[i][j+1]=-1;
        }
    }
    #ifdef DEBUG
    for (i=1;i<=n;i++){
        for(j=1;j<=n;j++)
            cout<<map[i][j]<<"  ";
        cout<<endl;
    }
    #endif // DEBUG



    for (i=1;i<=n;i++){
        int maxlen=-INF;
        for (j=1;j<=n;j++){
            if (map[j][i]>maxlen)   maxlen=map[j][i];
        }

        inside[i]=maxlen-1;
        for (j=0;j<=inside[i];j++){
            newpoint[i][j]=np++;
        }
    }

    #ifdef DEBUG
    for (i=1;i<=n;i++)
        cout<<inside[i]<<" ";
    cout<<endl;
    #endif // DEBUG



    for (i=1;i<=n;i++){
        for (j=1;j<=inside[i];j++){
            m.p[newpoint[i][j]][newpoint[i][j-1]]=1;
        }
    }
    for (i=1;i<=n;i++){
        for (j=1;j<=n;j++){
            if (map[i][j]!=-1){
                m.p[newpoint[i][0]][newpoint[j][map[i][j]-1]]=1;
            }
        }
    }

    #ifdef DEBUG
    cout<<np<<endl;

    for (i=1;i<np;i++){
        for(j=1;j<np;j++)
            cout<<m.p[i][j]<<" ";
        cout<<endl;
    }
    #endif // DEBUG



    m.x=m.y=np-1;
    matrix ans=quickmatrix(m,t);
    printf("%d\n",ans.p[newpoint[1][0]][newpoint[n][0]]);
    return 0;

}