题目链接：
给出一张 \(n\) 个点的有向图 \(G(V,E)\)。对于任意两个点 \(u,v\)（\(u\) 可以等于 \(v\)），\(u\) 向 \(v\) 的连边数为
\[
\sum_{i=1}^kout_{u,i}\times in_{v,i}
\]
其中 \(k\) 和数组 \(out,in\) 均已知，现在给出 \(m\) 个询问，每次询问给出三个参数 \(u, v, d\)，你需要回答从节点 \(u\) 出发，经过不超过 \(d\) 条边到达节点 \(v\) 的路径有多少种。答案模 \(10^9 + 7\)。
\(n \leq 1000, k \leq 20, m \leq 50, d < 2^{31}\)
时空限制：\(1\texttt s/256 \texttt{MB}\)
因为我们要求的是不超过 \(d\) 条边的方案数，因此新建一个虚拟点 \(n+1\)，并且令 \(out_{n+1,k+1}=in_{n+1,k+1}=1\)。
对于每个询问 \((u,v,d)\)，就临时令 \(out_{v,k+1}=1\)，其含义是到达了点 \(v\) 后就多走一步走到 \(n+1\)，并且 \(n+1\) 有一条连向自己的自环。这样就将询问转化成，从 \(u\) 走到 \(n+1\)，恰好走 \(d+1\) 步的方案数。
不妨令 \(n \leftarrow n+1,k\leftarrow k+1\)，\(A,B\) 对应如上更改后的 \(in\) 和 \(out\) 的矩阵，答案就是
\[
R=(A\times B^T)^{d+1}
\]
这个矩阵的第 \(u\) 行第 \(n\) 列，这样直接求是 \(\mathcal O(mn^3\log d)\) 的，非常逊色。
注意到 \(k \leq 20\)，以及 \(A,B\) 分别是 \(n\times k\) 和 \(k\times n\) 的矩阵，直接把它们乘起来变成 \(n\times n\) 的矩阵，然后快速幂显然非常浪费。实际上利用矩阵乘法的结合律，我们可以将 \(R\) 写成
\[
R=A(B^T\times A)^dB^T
\]
这样中间的矩阵 \(B^TA\) 就是 \(k\times k\) 的了，拿去矩阵快速幂就非常合适。
最后因为我们只要求 \(R_{u,n}\)，所以不需要把 \(R\) 求出来。预处理出 \(S=B^TA\) 之后，每次询问求出 \(S^d\)，然后 \(\mathcal O(k^2)\) 枚举一下就可以求出 \(R_{u,n}\)。
每次询问临时令 \(out_{v,k+1}=1\) 的时候，会更改 \(S\) 的某一列，这个只要 \(\mathcal O(k^2)\) 直接处理就好了。
总时间复杂度就是 \(\mathcal O\left(nk^2+mk^3\log d\right)\)。
#include <bits/stdc++.h>
template <class T>
inline void read(T &x)
{
static char ch;
while (!isdigit(ch = getchar()));
x = ch - '0';
while (isdigit(ch = getchar()))
x = x * 10 + ch - '0';
}
template <class T>
inline void putint(T x)
{
static char buf[15], *tail = buf;
if (!x)
putchar('0');
else
{
for (; x; x /= 10) *++tail = x % 10 + '0';
for (; tail != buf; --tail) putchar(*tail);
}
}
typedef long long s64;
const int mod = 1e9 + 7;
const s64 MOD = 4611000029LL * mod;
const int MaxK = 22;
const int MaxN = 1e3 + 5;
struct mat
{
int r, c;
int a[22][22];
mat(){}
mat(int _r, int _c):
r(_r), c(_c) {memset(a, 0, sizeof(a));}
inline void init()
{
for (int i = 1; i <= r; ++i)
a[i][i] = 1;
}
inline mat operator * (mat rhs) const
{
mat res(r, rhs.c);
for (int i = 1; i <= r; ++i)
for (int j = 1; j <= rhs.c; ++j)
{
s64 sum = 0;
for (int k = 1; k <= c; ++k)
{
sum += 1LL * a[i][k] * rhs.a[k][j];
if (sum >= MOD)
sum -= MOD;
}
res.a[i][j] = sum % mod;
}
return res;
}
inline mat operator ^ (int p) const
{
mat res(r, c), x = *this;
res.init();
for (; p; p >>= 1, x = x * x)
if (p & 1)
res = res * x;
return res;
}
}T;
int n, K, m;
int prod[MaxK][MaxK];
int A[MaxN][MaxK], B[MaxK][MaxN];
int main()
{
read(n), read(K);
for (int i = 1; i <= n; ++i)
{
for (int j = 1; j <= K; ++j)
{
read(A[i][j]);
A[i][j] %= mod;
}
for (int j = 1; j <= K; ++j)
{
read(B[j][i]);
B[j][i] %= mod;
}
}
++K, ++n;
A[n][K] = B[K][n] = 1;
for (int i = 1; i <= K; ++i)
for (int j = 1; j <= K; ++j)
{
s64 res = 0;
for (int k = 1; k <= n; ++k)
{
res += 1LL * B[i][k] * A[k][j];
if (res >= MOD)
res -= MOD;
}
prod[i][j] = res % mod;
}
T = mat(K, K);
read(m);
while (m--)
{
int u, v, d;
read(u), read(v), read(d);
A[v][K] = 1;
for (int i = 1; i <= K; ++i)
for (int j = 1; j <= K; ++j)
{
T.a[i][j] = prod[i][j];
if (j == K)
T.a[i][j] = (T.a[i][j] + 1LL * B[i][v] * A[v][j]) % mod;
}
T = T ^ d;
int res = 0;
for (int i = 1; i <= K; ++i)
{
int sum = 0;
for (int j = 1; j <= K; ++j)
sum = (sum + 1LL * T.a[i][j] * B[j][n]) % mod;
res = (res + 1LL * A[u][i] * sum) % mod;
}
A[v][K] = 0;
putint(res), putchar('\n');
}
return 0;
}原文：https://www.cnblogs.com/cyx0406/p/12322806.html