[ZJOI2004]沼泽鳄鱼

[ZJOI2004]沼泽鳄鱼

题目大意：

给定一个无向图，以及起点S和终点T。每单位时间可以从一个点走向联通的另一个点，问从起点S出发，经过K单位时间到达终点T的方案总数。
最多50个点，K小于等于10^9。
另：某一周期（周期时间为2,3或4）内一些点不能走到。

问题分析：

到达某一点的方案总数，需要想到由加法原理进行递推。记$f(u,k)$ 为k秒后到达u点的方案总数，那么有 $f(u,k)=\sum_{(u,v)\in E} f(v,k-1)$

数据范围点数较小，时间极大，立刻触发使用快速幂倍增的想法。由于递推式是有限点的线性递推式，可以写成矩阵形式$(M,其中M_{u,v}=M_{v,u}=1)$，因而确实可以使用矩阵快速幂。

再来考虑某一时刻有些点无法到达，即$f(u,k)=0$ ，这时只要令$M_{u,i}=0$ 即可。由于周期只有2,3,4三种情况，取其最小公倍数12，计算12秒内的状态转移矩阵，并相乘作为一个转移周期。这样做K/12次整周期快速幂和K%12次部分快速幂即得K次后的方案数。

注意矩阵是左乘，复杂度$O(logN*N^3)$

小结：

利用矩阵刻画固定的、线性的状态转移方程，以求利用矩阵快速幂进行加速，是动态规划的一种常见加速技巧。曾在AC自动机DP问题中出现(POj DNA Sequence)。
局限：

1. 状态数要少
2. 转移要是线性的
3. 转移方程是固定的，或者在周期的意义下是固定的

其它方法，状态数很多，但转移的状态有限，相对当前状态固定，转移次数较少：
例如： $f(i)=\sum_{k=1}^{l} f(i-k)$，生成函数$(x+x^2+...+x^l)^i$
因为 $(f(0)+f(1)x+...+f(n-l)x^{n-l}+...+f(n-1)x^{n-1})*(x+x^2+...+x^l)$对于$x^n$ 的系数就是$f(n)=\sum_{k=1}^{l} f(n-k)$
但是这个貌似只有理论上的用途，实际求f(n)，最后的多项式有l*k次，如果n很大，大于10^5，FFT+快速幂仍然吃不消（$O(logN*(NL)log(NL))$）,除非L很大，N较小，这个时候就可以做了。对于特定的题目n很大还是有用的，比如hdu6050，要先求通项再将特征根代入生成函数，因为只用到了理论概念并不涉及机器计算。

#include <bits/stdc++.h>
using namespace std;
const int mo=10000;
inline void add(int &x,int y){
  x+=y; if (x>=mo) x-=mo;
}
struct matrix{
  int n,m;
  int a[60][60];
  matrix operator*(const matrix &u) const {
    int i,j,k;
    matrix ans;
    ans.n=n;ans.m=m;
    for (i=1;i<=n;i++)
    for (j=1;j<=m;j++){
        ans.a[i][j]=0;
        for (k=1;k<=m;k++) add(ans.a[i][j],a[i][k]*u.a[k][j]%mo);
    }
    return ans;
  }

  void unit(){
    int i,j;
    for (i=1;i<=n;i++)
    for (j=1;j<=m;j++){
        a[i][j]= i==j ? 1 :0;
    }
  }
  void clear(int x,int y){
    int i,j;
    n=x;m=y;
    for (i=1;i<=n;i++)
    for (j=1;j<=m;j++){
        a[i][j]=0;
    }
  }
  matrix operator^(const int& index) const {
    int p=index;
    matrix ans;ans.n=n;ans.m=m;
    ans.unit();
    if (index==0) return ans;
    matrix t=*this;
    while (p){
      if (p&1) ans=ans*t;
      t=t*t;p>>=1;
    }
    return ans;
  }

};


matrix g[20],pe,re;
int in[60],p[5],out[60];

int main()
{
    freopen("swamp.in","r",stdin);
    freopen("swamp.out","w",stdout);
    int n,m,s,e,k,i,j,x,y,nfish,t;
    cin>>n>>m>>s>>e>>k;
    s++;e++;
    for (i=0;i<12;i++) g[i].clear(n,n);
    for (i=1;i<=m;i++){
      scanf("%d %d",&x,&y);
      x++;y++;
      if (g[0].a[x][y]) continue;
      for (j=0;j<12;j++) g[j].a[x][y]=g[j].a[y][x]=1;
    }
    cin>>nfish;
    while (nfish--){
      cin>>t;
      for (i=0;i<t;i++) {cin>>p[i];p[i]++;}
      for (i=0;i<12;i++){
        for (j=1;j<=n;j++) g[i].a[p[i%t]][j]=0;
      }
    }
    int index=k/12;
    k%=12;
    re.n=re.m=n;
    re.unit();pe=re;
    for (i=1;i<=12;i++){
       pe=g[i%12]*pe;
       if (i==k) re=pe;
    }

    pe=pe^index;
    pe=re*pe;
    for (i=1;i<=12;i++) in[i]=out[i]=0;
    in[s]=1;
    for(i=1;i<=n;i++)
      for (j=1;j<=n;j++)
        add(out[i],pe.a[i][j]*in[j]%mo);
    cout<<out[e]<<endl;
    return 0;
}