【SDOI2009】HH去散步（dp，矩阵快速幂）

看到题目第一眼，最暴力的想法就是设 \(dp(i,j)\) 表示以 \(a\) 为起点、\(i\) 为终点，路径长度为 \(j\) 的路径数。

容易得到状态转移方程：\(dp(u,j)=\sum_{(v,u)}dp(v,j-1)\)。

但是我们发现一个问题：题目中要求“人物不能立刻沿着刚刚走来的路走回”，但显然这种 dp 方式是不符合要求的，所以想怎么解决。

这里用到一个常用的小trick：我们先把双向边拆成两条单向边，共 \(2m\) 条，并编号。然后设 \(dp(i,j)\) 表示以 \(a\) 为起点、刚好走完编号为 \(i\) 的这条单向边，路径长度为 \(j\) 的路径数。你也可以理解成把边当成点。

状态转移方程同样也容易得到：

如图，设单向边 \((u,v)\) 的编号为 \(i\)，所有终点为 \(u\) 的单向边的边集为 \(S_i\)（也就是红圈圈出来的部分）。

显然有：\(dp(i,j)=\sum_{k\in S_i}dp(k,j-1)\)。

把 dp 的第二维滚动一下，即可 \(O(mt)\) 获得 30pts。

发现问题的瓶颈主要在 \(t\) 上，思考如何优化。

仔细观察 dp 式，发现和 【NOI Online #3 提高组】魔法值 的做法很相似。

于是联想到用矩阵快速幂来加速 dp 过程。

设计矩阵：（第二维被滚动掉了）

\(\begin{bmatrix} dp(1)\\ dp(2)\\ \cdots\\ dp(2m) \end{bmatrix}\)

然后根据上面推出的式子：\(dp(i,j)=\sum_{k\in S}dp(k,j-1)\)，可以设计出转移矩阵：（这里的 \(c_{i,j}\) 表示 \(j\) 是否属于集合 \(S_i\)，属于则 \(c_{i,j}=1\)，否则 \(c_{i,j}=0\)）

\(B=\begin{bmatrix} c_{1,1},&c_{1,2},&\dots,&c_{1,2m}\\ c_{2,1},&c_{2,2},&\dots,&c_{2,2m}\\ \cdots&\cdots&\cdots&\cdots\\ c_{2m,1},&c_{2m,2},&\dots,&c_{2m,2m}\\ \end{bmatrix}\)

设初始矩阵是 \(A\)（表示所有 \(dp(i,1)\) 的值），那么答案矩阵即为：\(Ans=A\times B^{t-1}\)，表示的是所有 \(dp(i,t)\) 的值。

那么答案就是终点为 \(b\) 的每一条单向边的 \(dp\) 值的总和。

总时间复杂度为 \(O(m^3\log t)\)，可以跑过。

---

顺带提一下我一开始的错误想法：设 \(dp(i,j)\) 表示以 \(a\) 为起点、\(i\) 为终点，路径长度为 \(j\) 的路径数，然后找到与 \(i\) 距离为 \(2\) 的点 \(v\)，用 \(dp(v,j-2)\) 来更新 \(dp(i,j)\)，这样就能保证不会出现 \(i\) 走到一个点又走回头路的情况了。但是下面这种图就能卡掉这个算法：

询问 \(t=3\)，\(a=2\)，\(b=3\)。显然答案为 \(0\)，但是上面那种做法的答案是 \(1\)。

至于为什么是错的我不详细阐述了，如果也是像我一样相同想法的可以结合这个例子自己回去模拟一下并想一想。

正解的代码及注释如下：

#include<bits/stdc++.h>
 
#define N 25
#define M 65
#define mod 45989
 
using namespace std;
 
int n,m,t,a,b;
int cnt=1,head[N],to[M<<1],nxt[M<<1];
 
struct Matrix
{
    int a[M<<1][M<<1];
    Matrix(){memset(a,0,sizeof(a));}
    void init(){for(int i=1;i<=cnt;i++) a[i][i]=1;}
}st,ch;
 
Matrix operator * (Matrix a,Matrix b)
{
    Matrix c;
    for(int i=1;i<=cnt;i++)
        for(int j=1;j<=cnt;j++)
            for(int k=1;k<=cnt;k++)
                c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%mod;
    return c;
}
  
Matrix poww(Matrix a,int b)
{
    Matrix ans;
    ans.init();
    while(b)
    {
        if(b&1) ans=ans*a;
        a=a*a;
        b>>=1;
    }
    return ans;
}
//以上是矩阵快速幂板子
 
void adde(int u,int v)
{
    to[++cnt]=v;
    nxt[cnt]=head[u];
    head[u]=cnt;
}
 
void init()
{
    for(int i=head[a];i;i=nxt[i])
        st.a[1][i]=1;//设计初始矩阵
    for(int i=2;i<=cnt;i++)//枚举每条边i，这里其实是先枚举Si再去更新i
    {
        int u=to[i];
        for(int j=head[u];j;j=nxt[j])
            if(j!=(i^1))
                ch.a[i][j]++;//设计转移矩阵
    }
}
 
int main()
{
    scanf("%d%d%d%d%d",&n,&m,&t,&a,&b);
    a++,b++;
    for(int i=1;i<=m;i++)
    {
        int u,v;
        scanf("%d%d",&u,&v);
        u++,v++;
        adde(u,v),adde(v,u);//两条单向边
    }
    init();
    Matrix tmp=st*poww(ch,t-1);
    int ans=0;
    for(int i=2;i<=cnt;i++)
        if(to[i]==b)//枚举每一条终点是b的单向边
            ans=(ans+tmp.a[1][i])%mod;
    printf("%d\n",ans);
    return 0;
}