Description
题目链接：
- BZOJ #4362
- 印象里是 FJOI2018 day1T3 的原题，具体题意是否是这样，以及数据范围不太记得了，但记得做法差不多。
给出一张 (n) 个点的有向图 (G(V,E))。对于任意两个点 (u,v)（(u) 可以等于 (v)），(u) 向 (v) 的连边数为
其中 (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 exttt s/256 exttt{MB})
Solution
因为我们要求的是不超过 (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,kleftarrow k+1)，(A,B) 对应如上更改后的 (in) 和 (out) 的矩阵，答案就是
这个矩阵的第 (u) 行第 (n) 列，这样直接求是 (mathcal O(mn^3log d)) 的，非常逊色。
注意到 (k leq 20)，以及 (A,B) 分别是 (n imes k) 和 (k imes n) 的矩阵，直接把它们乘起来变成 (n imes n) 的矩阵，然后快速幂显然非常浪费。实际上利用矩阵乘法的结合律，我们可以将 (R) 写成
这样中间的矩阵 (B^TA) 就是 (k imes 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 Oleft(nk^2+mk^3log d ight))。
#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('
');
}
return 0;
}