BZOJ 1875(DP+矩阵快速幂)

题面

传送门

分析

容易想到根据点来dp，设dp[i][j]表示到i点路径长度为j的方案数
状态转移方程为dp[i][k]=∑(i,j)∈Edp[j][k−1]$dp[i][k]=\sum _{\left(i,j\right)\in E}dp\left[j\right][k-1]$
但这样得出的结果是错误的，因为它没有考虑一个点经过多次的情况

因此，我们按边来dp，因为每条边只能经过一次，所以不会出现上面的问题
将每条无向边拆成两条有向边
设dp[i][j]表示当前走到到编号为i的边路径长度为j的方案数
dp[i][k]=∑from[i]=to[j],i与j为不为一对反向边dp[j][k−1]$$\begin{gathered} dp[i][k]=\sum _{from[i]=to[j], \\ i与j为不为一对反向边}dp[j][k-1] \end{gathered}$$
这样的时间复杂度为O(tm)$O(tm)$，显然是会超时的

注意到矩阵乘法优化dp的条件
前一个阶段到后一个阶段的映射是线性的，并且这个映射是不变的
此题中k为阶段，可以发现映射显然是不变的常量
我们用一个m×m$m\times m$矩阵
对于边i,j如果满足to[i] == from[j] &&i,j不为反向边这个条件，那么(i,j)是1，反之就是0
初始的矩阵为1×m$1\times m$的矩阵，对于从起点出发的每一条边i，我们将(1,i）设为1，反之为0.
对于答案矩阵，我们枚举指向终点的每一条边，将它在矩阵中对应位置的值相加即可

时间复杂度O(m3log2t)$O(m^{3}lo{g}_{2}t)$

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#define mod 45989
#define maxn 45
#define maxm 65
#define maxs maxm*2
using namespace std;
int n,m,k,s,t;
struct edge{
    int from;
    int to;
    int next;
}E[maxm<<1];
int head[maxn];
int sz=1;//从1开始存储，则第i条边的反向边编号为i^1
void add_edge(int u,int v){
    sz++;
    E[sz].from=u;
    E[sz].to=v;
    E[sz].next=head[u];
    head[u]=sz;
}

struct matrix{
    int n;
    int m;
    long long a[maxs][maxs];
    void print(){
        for(int i=1;i<=n;i++){
            for(int j=1;j<=m;j++){
                printf("%d ",a[i][j]);
            } 
            printf("
");
        }
        printf("
");
    }
    matrix(){
        memset(a,0,sizeof(a));
    }
    matrix(int x,int y){
        n=x;
        m=y;
        memset(a,0,sizeof(a));
    }
    friend matrix operator *(matrix u,matrix v){
        matrix ans=matrix(u.n,v.m);
        for(int i=1;i<=u.n;i++){
            for(int j=1;j<=v.m;j++){
                for(int k=1;k<=u.m;k++){
                    ans.a[i][j]+=u.a[i][k]*v.a[k][j]%mod;
                }
                ans.a[i][j]%=mod;
            }
        }
        return ans;
    }
}; 
matrix fast_pow(matrix x,int k){
    matrix ans=matrix(x.n,x.m);
    for(int i=1;i<=x.n;i++){
        ans.a[i][i]=1;
    }
    while(k){
        if(k&1) ans=ans*x;
        x=x*x;
        k>>=1;
    }
    return ans;
}
int main(){
    scanf("%d %d %d %d %d",&n,&m,&k,&s,&t);
    int u,v;
    for(int i=1;i<=m;i++){
        scanf("%d %d",&u,&v);
        add_edge(u,v);
        add_edge(v,u);
    }
    matrix f=matrix(sz,sz),p=matrix(sz,sz);
    for(int i=head[s];i;i=E[i].next){
        f.a[1][i]=1;
    }
    for(int i=2;i<=sz;i++){
        for(int j=2;j<=sz;j++){
            if(E[i].to==E[j].from&&i!=(j^1)){
                p.a[i][j]=1;
            }
        }
    }
//  f.print();
//  p.print();
    f=f*fast_pow(p,k-1);
//  f.print();
    long long ans=0;
    for(int i=head[t];i;i=E[i].next){
        ans=ans+f.a[1][i^1];
        ans%=mod;
    }
    printf("%lld
",ans%mod);
}