自动搬运

查看原文

来自洛谷，原作者为

	Caro23333 Dream Theater 脑残粉

搬运于2025-08-24 21:56:31，当前版本为作者最后更新于2020-07-04 14:24:53，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

---

咕了几天终于来了~

这个题首先看数据范围发现是 $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;
}