[BZOJ1898][ZJOI2005]沼泽鳄鱼（矩阵乘法）

先把无向边拆成两条有向边。
如果不考虑食人鱼，那么此题就是一个简单的矩乘问题。建立矩阵$P$，如果存在边$i->j$则$P[i][j]=1$，否则$P[i][j]=0$。
最后$P^K[Start][End]$就是最终结果。
考虑食人鱼。注意到题目里食人鱼的运动周期长度只有$2$，$3$，$4$，可以得出，第$t$个时间单位时是否可以在这个点上（没有食人鱼经过），取决于$t\mod 12$的值。所以这里，先预处理出一个数组$vis[j][i]$，如果当前的时间取模$12$的值为$j$，$vis[j][i]$为假表示此时可以在点$i$上，为真表示此时不能在点$i$上。
构造出$12$个矩阵$P_0,P_1,P_2,...,P_{11}$，分别储存当前单位时间模$12$的值为 $0,1,2,...,11$时的邻接矩阵。
对于一条边$i->j$，如果$vis[k][j]$为假，那么$P_k[i][j]=1$，否则如果不存在边$i->j$或者$vis[k][j]$为真，那么$P_k[i][j]=0$。
此时令矩阵$Q=P_1*P_2*P_3*...*P_{11}*P_0$。
并求矩阵$F=Q^{\lfloor\frac{K}{12}\rfloor}$，
就得到了在恰好$\lfloor\frac{K}{12}\rfloor*12$的时间内任意两点互相到达的方案数。由于剩下的不会超过$12$个时间单位，所以对于剩余的时间单位，暴力转移即可。
代码：

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
inline int read() {
    int res = 0; bool bo = 0; char c;
    while (((c = getchar()) < '0' || c > '9') && c != '-');
    if (c == '-') bo = 1; else res = c - 48;
    while ((c = getchar()) >= '0' && c <= '9')
        res = (res << 3) + (res << 1) + (c - 48);
    return bo ? ~res + 1 : res;
}
const int N = 55, M = 3005, PYZ = 1e4;
struct cyx {
    int m, n, v[N][N];
    cyx() {}
    cyx(int _m, int _n) :
        m(_m), n(_n) {memset(v, 0, sizeof(v));}
    friend inline cyx operator * (cyx a, cyx b) {
        cyx res = cyx(a.m, b.n); int i, j, k;
        for (i = 1; i <= res.m; i++) for (j = 1; j <= res.n; j++)
        for (k = 1; k <= a.n; k++)
            (res.v[i][j] += a.v[i][k] * b.v[k][j]) %= PYZ;
        return res;
    }
    friend inline cyx operator ^ (cyx a, int b) {
        int i, d = b; cyx c = a, res = cyx(a.m, a.n);
        for (i = 1; i <= res.m; i++) res.v[i][i] = 1;
        while (d) {
            if (d & 1) res = res * c;
            c = c * c;
            d >>= 1;
        }
        return res;
    }
} P[15], Q, F;
int n, m, S, T, K, nf, ecnt, nxt[M], adj[N], go[M];
bool vis[15][N];
void add_edge(int u, int v) {
    nxt[++ecnt] = adj[u]; adj[u] = ecnt; go[ecnt] = v;
}
int main() {
    int i, j, k, x, y; n = read(); m = read(); S = read() + 1;
    T = read() + 1; K = read(); for (i = 1; i <= m; i++) {
        x = read() + 1; y = read() + 1;
        add_edge(x, y); add_edge(y, x);
    }
    nf = read(); for (i = 1; i <= nf; i++) {
        x = read(); for (j = 0; j < x; j++) {
            y = read() + 1;
            for (k = j; k < 12; k += x)
                vis[k][y] = 1;
        }
    }
    Q = cyx(n, n); for (i = 1; i <= n; i++) Q.v[i][i] = 1;
    for (j = 0; j < 12; j++) {
        P[j] = cyx(n, n);
        for (i = 1; i <= n; i++)
            for (int e = adj[i]; e; e = nxt[e])
                if (!vis[j][go[e]]) P[j].v[i][go[e]] = 1;
        if (j) Q = Q * P[j];
    }
    Q = Q * P[0]; F = Q ^ (K / 12);
    for (j = 1; j <= K % 12; j++) F = F * P[j];
    printf("%d\n", F.v[S][T]);
    return 0;
}