【Matrix】矩阵乘法

矩阵乘法: 

Mat Mult(Mat a, Mat b){  // Mat 为自定义类型    Mat c;    memset(c.m, 0, sizeof(c.m));    c.r= a.r;  c.c=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.m[i][j]= (c.m[i][j]+a.m[i][k]*b.m[k][j]);                    }            }    return c;}

快速幂:

Mat Calc(int x){     if(x==1) {              return yuan;  // 原始矩阵              }     Mat ret= Calc(x/2);     ret= Mult(ret, ret);     if(x%2==1) ret=Mult(ret, yuan);     return ret;     }

 

snack：

给定一个有向图， 问在规定时间从起点出发沿边行走并收集所有物品(边上有一定的物品，收集花费时间2，否则花费时间1)的总方案数。

时间  10^9,   物品4,  顶点数为25

考虑DP：  f[i,  j, k](i   所在顶点,   j 经过时间,     k收集物品的状态)  =  { f[v,      j-1,    k]  }+  {  f[v,  j-2,  x] } 分别为收集该边和不收集改边，  由于时间一维巨大， 考虑使用矩阵乘法。

定义一维矩阵G= {sum,    f1,   f2,   f3 …, fn ,   g1,   g2,   g3, … gn} 为到顶点1的总方案数，  时间i的每个顶点的每种状态, i-1时间的顶点与状态。

于是可以构造一个二位矩阵M { 符合上面的运算的 }  

分析时间复杂度  25* 16 * 2+ 1= 801    O(N^3)的时间复杂度接受不了。

 

所以可以考虑容斥原理：  把ABCD含都加进方案里，再把ABC, ABD,  ACD , BCD减去 ……

这样构造矩阵的时候比较和谐， 时间复杂度也变成了 O(16 *  (25*2+1)^3 *  log2(10^9) )

代码：

#include<cstdio>#include<cstring>using namespace std;#define N 25#define M 51#define Edge 500#define Mod 5557struct Mat{ int m[M][M], r , c ;};const int data[16]={1, -1, -1 , 1, -1, 1, 1,-1, -1, 1, 1, -1, 1, -1, -1, 1};Mat yuan, I;int n,m,T;int state[Edge], x[Edge], y[Edge];int Change(char c){    if (c=='B') return 1;     if (c=='J') return 2;    if (c=='M') return 4;     if (c=='P') return 8;    }void Make_State(int &s){     char c;       while(c=getchar(), c<'A' ||  c>'Z');     s=Change(c);  while(c=getchar(), c>='A'&& c<='Z') s+=Change(c);     }void Init(){     scanf("%d%d",&n,&m);     for(int i=0;i<m;i++){             scanf("%d%d",&x[i],&y[i]);             Make_State(state[i]);             }     scanf("%d",&T);     }     Mat Mult(Mat a, Mat b){    Mat c;    memset(c.m, 0, sizeof(c.m));    c.r= a.r;  c.c=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.m[i][j]= (c.m[i][j]+a.m[i][k]*b.m[k][j]) % Mod;                    }            }    return c;}Mat Calc(int x){     if(x==1) {              return yuan;              }     Mat ret= Calc(x/2);     ret= Mult(ret, ret);     if(x%2==1) ret=Mult(ret, yuan);     return ret;     }     void mat(int S){     memset(yuan.m, 0, sizeof(yuan.m));     yuan.r=yuan.c=2*n;     yuan.m[0][0]= yuan.m[0][1]=1;     for(int i=1; i<=n; i++){yuan.m[n+i][i]=1;}     for(int i=0; i<m; i++){             int u=x[i], v=y[i];             yuan.m[v][u] =1;               if((state[i] | S)==S) yuan.m[v][u+n] =1;             }     }     int Solve(){    Mat amat=Calc(T+1);    memset(I.m, 0, sizeof(I.m));    I.r=n*2, I.c=0, I.m[0][0]=0; I.m[1][0]=1;    Mat ans= Mult(amat, I);    return ans.m[0][0];    }         void Work(){     int Ans=0;     for(int i=0;i<=15;i++){             mat(15-i);             int tt=Solve();             Ans= (Ans+ data[i]*tt ) % Mod;             }     Ans= (Ans+Mod)%Mod;     printf("%d\n",Ans);     }int main(){    Init();    Work();    return 0;}