题解 P4102 【[HEOI2014]林中路径 】

· · 题解

咕了几天终于来了~

这个题首先看数据范围发现是 n \le 100 的点对之间路径计数，其中带有路径长度的限制，那很明显就跟邻接矩阵的幂脱不了关系，那么主要的难点就在于如何处理路径长度的平方这件事。

大体的思路是递归，先求出 m = \lfloor \frac{k}{2} \rfloor 时的答案，再利用 m 时的答案推出 k 的答案。考虑维护三个矩阵 A,B,C 分别代表当前 所有路径长度的 0,1,2 次方之和。

简单起见，这里只讨论 k 为偶数，即 k=2m 时的情况，奇数情况是基本同理的。考虑对于某一个有序点对 (u,v)，令 a_i 为从 u 到 v 长度恰为 i 的路径数量，则有：

A_{u,v} = \sum_{i=1}^m a_i,\ B_{u,v} = \sum_{i=1}^m ia_i,\ C_{u,v} = \sum_{i=1}^m i^2a_i

而对于最长路径限制为 k=2m 时的 A,B,C（不妨记为 A',B',C'）有：

A'_{u,v} = \sum_{i=1}^{2m} a_i = A_{u,v}+\sum_{i=1}^m a_{i+m} B'_{u,v} = \sum_{i=1}^{2m} ia_i = B_{u,v}+\sum_{i=1}^m ia_{i+m}+m\sum_{i=1}^m a_{i+m} C'_{u,v} = \sum_{i=1}^{2m} i^2a_i = C_{u,v}+\sum_{i=1}^m i^2a_{i+m}+2m\sum_{i=1}^m ia_{i+m}+m^2\sum_{i=1}^m a_{i+m}

不妨令原图的邻接矩阵为 G。我们所熟知的是，对于任意的 i，全体 a_i 到全体 a_{i+m} 的变化是可以用全体 a_i 对应的矩阵乘上 G^m 来刻画的，那么对于 a_i 的前缀和（即 A 到 A' 的变化）其实也是完全相同的。于是对照上式即得：

A' = A+A \cdot G^m

同理更有

B' = B+B \cdot G^m+m \cdot A \cdot G^m C' = C+C \cdot G^m+2m \cdot B \cdot G^m+m^2 \cdot A \cdot G^m

那么只要在递归过程中记录 G^m 即可。k 为奇数的情况，只要再加上一次 m=1 即可，是完全类似的。复杂度 O(n^3 \log k + q)。可能略卡常。

代码：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#define mod 1000000007

using namespace std;
typedef long long ll;
const int MAXN = 105; 
inline int readint()
{
    int res = 0, f = 1;
    char c = 0;
    while(!isdigit(c))
    {
        c = getchar();
        if(c=='-')
            f = -1;
    }
    while(isdigit(c))
        res = res*10+c-'0', c = getchar();
    return res*f;
}
int n,m,k,q;
inline int modd(int x)
    { return x>=mod?x-mod:x; } 
struct Mat
{
    int a[MAXN][MAXN];
    Mat()
        { memset(a,0,sizeof(a)); } 
    Mat operator+(const Mat t) const
    {
        Mat res;
        for(int i = 1; i<=n; i++)
            for(int j = 1; j<=n; j++)
                res.a[i][j] = modd(a[i][j]+t.a[i][j]);
        return res;
    }
    Mat operator*(const Mat t) const
    {
        Mat res;
        for(int i = 1; i<=n; i++)
            for(int j = 1; j<=n; j++)
                if(t.a[i][j])
                    for(int k = 1; k<=n; k++)
                        res.a[i][k] = (1ll*res.a[i][k]+1ll*a[j][k]*t.a[i][j])%mod;
        return res; 
    }
    Mat operator*(const int t) const
    {
        Mat res;
        for(int i = 1; i<=n; i++)
            for(int j = 1; j<=n; j++)
                res.a[i][j] = 1ll*a[i][j]*t%mod;
        return res; 
    }
}G,c,ans0,ans1,ans2;
inline Mat unit()
{
    Mat res;
    for(int i = 1; i<=n; i++)
        res.a[i][i] = 1;
    return res;
} 
inline void solve(int t)
{
    if(t==1)
    {
        c = G, ans0 = unit();
        return;
    }
    solve(t/2);
    int x = t/2;
    ans2 = ans2+ans2*c+ans1*c*(2*x)+ans0*c*(1ll*x*x%mod);
    ans1 = ans1+ans1*c+ans0*c*x;
    ans0 = ans0+ans0*c;
    c = c*c;
    if(t&1)
    {
        ans0 = ans0+c;
        ans1 = ans1+c*(t-1);
        ans2 = ans2+c*(1ll*(t-1)*(t-1)%mod);
        c = c*G;
    }
}

int main()
{
    n = readint(), m = readint(), k = readint(), q = readint();
    for(int i = 1; i<=m; i++)
    {
        int u = readint(), v = readint();
        G.a[u][v]++;
    }
    solve(k+1);
    while(q--)
    {
        int x = readint(), y = readint();
        printf("%d\n",ans2.a[x][y]); 
    }
    return 0;
}