1 计算路径数
对于一个边权为 \(\bf1\) 的图（有向或无向）的邻接矩阵 \(G\)，考虑它的幂的意义是什么。
设 \(c_{i,j}\) 表示 \(i\) 和 \(j\) 之间是否有连边，则有
\[G=\begin{bmatrix} c_{1,1} & c_{1,2} & \cdots & c_{1,n} \\ c_{2,1} & c_{2,2} & \cdots & c_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ c_{n,1} & c_{n,2} & \cdots & c_{n,n} \end{bmatrix} \]那么
\[\left(G^2\right)_{i,j}=\sum_{k=1}^nc_{i,k}\times c_{k,j} \]把 \(c_{i,k}\) 和 \(c_{k,j}\) 都看作是布尔值，那么 \(c_{i,k}\times c_{k,j}\) 的意思就是 \(i\to k\to j\) 有连边。换句话说，\(i\) 到 \(j\) 存在长度为 \(2\) 的一条路径。
因此 \(\left(G^2\right)_{i,j}\) 就表示 \(i\) 到 \(j\) 长度为 \(2\) 的路径数。
如果再乘一次呢？
\[\left(G^3\right)_{i,j}=\sum_{k=1}^n\left(G^2\right)_{i,k}\times c_{k,j} \]可以发现，相当于枚举 \(k\) 为中转点，所以 \(\left(G^3\right)_{i,j}\) 表示 \(i\) 到 \(j\) 长度为 \(3\) 的路径数。
因此我们得到结论：
对于一个边权为 \(1\) 的图，邻接矩阵 \(G\) 的 \(k\) 次幂的一个元素 \(\left(G^k\right)_{i,j}\) 表示 \(i\) 到 \(j\) 长度为 \(k\) 的路径数。
P4159 [SCOI2009] 迷路
题意
该有向图有 \(n\) 个节点，节点从 \(1\) 至 \(n\) 编号，windy 从节点 \(1\) 出发，他必须恰好在 \(t\) 时刻到达节点 \(n\)。
通过某有向边的时间为给定的时间，范围为 \(1\) 到 \(9\) 的整数。
现在给出该有向图，你能告诉 windy 总共有多少种不同的路径吗？答案对 \(2009\) 取模。
对于 \(100\%\) 的数据，保证 \(2\le n\le10\)，\(1\le t\le10^9\)。
思路
如果边权为 \(1\)，就变成了上面讨论的问题。然后下面的拆边就需要智慧了。。
解决方法是：设一条边 \(u\to v\) 的权值是 \(w\)，那么我们建立 \(w-1\) 个虚点 \(v_{w-1},\dots,v_1\)，然后在新图上建立权值为 \(1\) 的边 \(u\to v_{w-1},v_{w-1}\to v_{w-2},\dots,v_2\to v_1,v_1\to v\)。这样一来，只有从 \(u\) 经过 \(w\) 条边到达 \(v\) 才算真正走完原来权值为 \(w\) 的有向边。
设新图的邻接矩阵为 \(G\)，答案即为 \(\left(G^t\right)_{1,n}\)。
代码
#include <cctype>
#include <cstring>
#include <iostream>
#define f(x, y, z) for (int x = (y); x <= (z); ++x)
using namespace std;
int constexpr N = 99;
int constexpr MOD = 2009;
int n, t;
inline int &AddEq(int &a, int const &b) { return (a += b) >= MOD ? (a -= MOD) : a; }
struct Matrix {
int a[N][N];
Matrix() { memset(a, 0, sizeof a); }
inline void clear() { memset(a, 0, sizeof a); f(i, 1, n) a[i][i] = 1; }
Matrix operator*(Matrix const &b) {
Matrix c = Matrix();
f(i, 1, n) f(j, 1, n) f(k, 1, n)
AddEq(c.a[i][j], a[i][k] * b.a[k][j] % MOD);
return c;
}
Matrix operator^(int x) {
Matrix res, p = *this;
res.clear();
while (x) {
if (x & 1) res = res * p;
p = p * p;
x >>= 1;
}
return res;
}
} mat;
inline int pos(int const &i, int const &j) { return i + n * j; }
int read() {
char ch;
while (!isdigit(ch = getchar()));
return ch - '0';
}
signed main() {
cin >> n >> t;
mat = Matrix();
f(i, 1, n) {
f(j, 1, 8) mat.a[pos(i, j)][pos(i, j - 1)] = 1;
f(j, 1, n) {
int w = read();
if (w) mat.a[i][pos(j, w - 1)] = 1;
}
}
n *= 9;
mat = mat ^ t;
cout << mat.a[1][n / 9] << '\n';
return 0;
}
P2151 [SDOI2009] HH去散步
题意
给定一张 \(N\) 个点 \(M\) 条边的无向图，边长度为 \(1\)，可能有 重边，且 不允许连续重复经过同一条边。
求长度为 \(t\) 的从 \(A\) 到 \(B\) 的路径共有多少种。答案模 \(45989\)。
对于 \(100\%\) 的数据，\(N\le50\)，\(M\le60\)，\(t\le2^{30}\)，\(0\le A,B\)。
思路
如果没有重边和不允许走回头路的限制，这道题仍然是一道裸题。考虑如何避免重边和回头路的影响。
接下来又是人类智慧环节：
我们把无向边拆成两条有向边，然后为每条有向边在新图中建立一个节点。
例如，存在 \(i-j-k\) 两条相邻的无向边，则在新图中建立 \((i,j)\to(j,k)\) 和 \((k,j)\to(j,i)\) 四个节点两条边。
统计答案时对于所有 \(i\) 和 \(j\)，把 \((A,i)\) 到 \((j,B)\) 的所有路径累加进答案即可。
代码
#include <cstring>
#include <iostream>
#define f(x, y, z) for (int x = (y); x <= (z); ++x)
using namespace std;
int constexpr N = 60, M = 130;
int constexpr MOD = 45989;
int n, m, t, a, b;
inline int &AddEq(int &a, int const &b) { return (a += b) >= MOD ? (a -= MOD) : a; }
int head[N], cnt = 1;
struct Edge {
int to, nxt;
} e[M];
inline void add(int from, int to) {
e[++cnt].to = to, e[cnt].nxt = head[from], head[from] = cnt;
return;
}
struct Matrix {
int a[M][M];
} I, g, ans, E;
Matrix Mul(Matrix a, Matrix b) {
Matrix c = E;
f(i, 1, cnt) f(j, 1, cnt) f(k, 1, cnt)
AddEq(c.a[i][j], a.a[i][k] * b.a[k][j] % MOD);
return c;
}
Matrix ksm(Matrix a, int x) {
Matrix res = I;
while (x) {
if (x & 1) res = Mul(res, a);
a = Mul(a, a);
x >>= 1;
}
return res;
}
void Init() {
for (int i = head[a]; i; i = e[i].nxt)
ans.a[1][i] = 1; //初始矩阵
f(i, 1, cnt)
for (int j = head[e[i].to]; j; j = e[j].nxt)
if (j ^ i ^ 1) ++g.a[i][j]; //转移矩阵
return;
}
signed main() {
ios::sync_with_stdio(false);
cin.tie(0), cout.tie(0);
cin >> n >> m >> t >> a >> b;
++a, ++b;
int u, v;
f(i, 1, m) {
cin >> u >> v;
++u, ++v;
add(u, v), add(v, u);
}
f(i, 1, cnt) I.a[i][i] = 1;
memset(E.a, 0, sizeof E.a);
Init();
ans = Mul(ans, ksm(g, t - 1));
int tot = 0;
f(i, 1, cnt) if (e[i].to == b) AddEq(tot, ans.a[1][i]);
cout << tot << '\n';
return 0;
}
2 广义矩阵乘法优化 DP
对于 \(p\times n\) 的矩阵 \(A\) 和 \(n\times q\) 的矩阵 \(B\)，定义广义矩阵乘法 \(\circ\)，得到的乘积 \(A\circ B\) 为 \(p\times q\) 矩阵，满足
\[(A\circ B)_{i,j}=\bigoplus_{k=1}^nA_{i,k}\otimes B_{k,j} \]其中 \(\oplus\) 和 \(\otimes\) 为两种二元运算，且满足：
- \(\otimes\) 有交换律和结合律，即 \(a\otimes b=b\otimes a\)，\((a\otimes b)\otimes c=a\otimes(b\otimes c)\)。
- \(\otimes\) 对 \(\oplus\) 有分配律，即 \((a\oplus b)\otimes c=(a\otimes c)\oplus(b\otimes c)\)。
这样得到的广义矩阵乘法和普通的矩阵乘法一样，不满足交换律，但满足结合律。
证明 下面把 \(M_1\circ M_2\) 省略为 \(M_1M_2\)。设 \(A\) 为 \(n\times m\) 矩阵，\(B\) 为 \(m\times p\) 矩阵， \(C\) 为 \(p\times q\) 矩阵。
\[\begin{aligned} &\begin{aligned} ((AB)C)_{i,j}&=\bigoplus_{k=1}^p(AB)_{i,k}\otimes C_{k,j} \\ &=\bigoplus_{k=1}^p\left(\bigoplus_{l=1}^mA_{i,l}\otimes B_{l,k}\right)\otimes C_{k,j}\qquad(1) \end{aligned}\\ &\begin{aligned} (A(BC))_{i,j}&=\bigoplus_{l=1}^mA_{i,l}\otimes(BC)_{l,j} \\ &=\bigoplus_{l=1}^mA_{i,l}\otimes\left(\bigoplus_{k=1}^pB_{l,k}\otimes C_{k,j}\right)\qquad(2) \end{aligned} \end{aligned} \]如果 \(\oplus\) 和 \(\otimes\) 满足条件，那么 \((1)\) 和 \((2)\) 都等于
\[\bigoplus_{k=1}^p\bigoplus_{l=1}^mA_{i,l}\otimes B_{l,k}\otimes C_{k,j} \]因此广义矩阵乘法有结合律。\(\kern{10pt}\square\)
于是我们可以利用快速幂来加速矩阵乘法，从而加速整个算法。
P2886 [USACO07NOV]Cow Relays G
题意
给定一张 \(T\) 条边的无向连通图，求从 \(S\) 到 \(E\) 经过 \(N\) 条边的最短路长度。
\(1\le N\le10^6\)，\(2\le T\le100\)，\(1\le u,v,w\le1000\)。
思路
首先将点离散化处理。因为图保证连通，所以点数 \(n\le T+1\)。
考虑 DP，设 \(f_x(i,j)\) 表示经过 \(x\) 条边从 \(i\) 到 \(j\) 的最短路长度。那么有
\[f_x(i,j)=\min_{(k,j)\in E}\left\{f_{x-1}(i,k)+w(k,j)\right\} \]其中 \((k,j)\in E\) 表示 \(k\) 和 \(j\) 有连边，\(w(k,j)\) 表示边 \((k,j)\) 的长度。由于我们枚举了所有的 \(k\)，因此这里并不需要枚举所有 \(f_{x-y}(i,k)+f_y(k,j)\)（\(y<x\)）。
我们发现，这与矩阵乘法的形式很相似。我们定义广义矩阵乘法
\[(AB)_{i,j}=\min_{k=1}^n\left\{A_{i,k}+B_{k,j}\right\} \]那么对于图的邻接矩阵 \(G\)，有
\[\left(G^k\right)_{i,j}=f_{k}(i,j) \]于是求出邻接矩阵的 \(N\) 次幂，输出 \(\left(G^N\right)_{S,E}\) 即可。
代码
注意初始化邻接矩阵（因为求的是最小值）。代码中的 k 即为题中的 \(N\)。
#include <cstring>
#include <iostream>
#include <algorithm>
#define f(x, y, z) for (int x = (y); x <= (z); ++x)
using namespace std;
int constexpr T = 110, N = 1010;
int k, t, s, e, n, raw[N], idx[N];
struct Edge {
int u, v, w;
} edge[T];
struct Matrix {
int x[T][T];
Matrix() { memset(x, 0x3f, sizeof x); }
Matrix operator*(Matrix const &o) {
Matrix res = Matrix();
f(i, 1, n) f(j, 1, n) f(k, 1, n)
res.x[i][j] = min(res.x[i][j], x[i][k] + o.x[k][j]);
return res;
}
} a, ans;
void ksm(int x) {
while (x) {
if (x & 1) ans = ans * a;
a = a * a;
x >>= 1;
}
return;
}
signed main() {
ios::sync_with_stdio(false);
cin.tie(0), cout.tie(0);
cin >> k >> t >> s >> e;
f(i, 1, t) {
cin >> edge[i].w >> edge[i].u >> edge[i].v;
raw[++n] = edge[i].u, raw[++n] = edge[i].v;
}
sort(raw + 1, raw + n + 1);
n = unique(raw + 1, raw + n + 1) - raw - 1;
ans = Matrix();
f(i, 1, n) idx[raw[i]] = i, ans.x[i][i] = 0;
s = idx[s], e = idx[e];
f(i, 1, t) {
int x = idx[edge[i].u], y = idx[edge[i].v], z = edge[i].w;
a.x[x][y] = a.x[y][x] = min(a.x[x][y], z);
}
ksm(k);
cout << ans.x[s][e] << '\n';
return 0;
}
P6569 [NOI Online #3 提高组] 魔法值
题意
H 国的交通由 \(n\) 座城市与 \(m\) 条道路构成，城市与道路都从 \(1\) 开始编号，其中 \(1\) 号城市是 H 国的首都。H 国中一条道路将把两个不同城市直接相连，且任意两个城市间至多有一条道路。
H 国是一个信奉魔法的国家，在第 \(j\) 天，\(i\) 号城市的魔法值为 \(f_{i,j}\)。H 国的魔法师已观测到第 0 天时所有城市的魔法值 \(f_{i,0}\)，且他们还发现，之后的每一天每个城市的魔法值，都将会变为所有与该城市直接相连的城市的前一天魔法值的异或值，即
\[f_{x,j}=f_{v_1,j-1}\oplus f_{v_2,j-1}\oplus\cdots\oplus f_{v_k,j-1} \]其中 \(j\ge1\)，\(v_1,v_2,\dots,v_k\) 是所有与 \(x\) 号城市直接相连的城市，\(\oplus\) 为异或运算。
现在 H 国的国王问了你 \(q\) 个问题，对于第 \(i\)（\(1\le i\le q\)）个问题你需要回答：第 \(a_i\) 天时首都的魔法值是多少。
\(1\le n,q\le100\)，\(1\le m\le\frac{n(n-1)}2\)，\(1\le a_i<2^{32}\)，\(0\le f_i<2^{32}\)。
思路
下面用 \(\oplus\) 表示异或。
在题目给的式子中，我们把 \(x\) 与某一点 \(y\) 是否连通用布尔值 \(c_{x,y}\) 表示，那么有
\[f_{x,j}=\bigoplus_{y=1}^nc_{x,y}\times f_{y,j-1} \]注意到，这个式子与矩阵乘法的式子比较相似。我们可以尝试定义一种广义矩阵乘法。
设
\[A_j=\begin{bmatrix} f_{1,j}\\f_{2,j}\\\vdots\\f_{n,j} \end{bmatrix}\qquad C=\begin{bmatrix} c_{1,1}&c_{1,2}&\cdots&c_{1,n}\\ c_{2,1}&c_{2,2}&\cdots&c_{2,n}\\ \vdots&\vdots&\ddots&\vdots\\ c_{n,1}&c_{n,2}&\cdots&c_{n,n}\\ \end{bmatrix} \]其中 \(A_j\) 是 \(\mathbb R^n\) 上的向量，\(C\) 是 \(n\) 阶方阵。定义广义矩阵乘法
\[(C\times A_{j})_i=\bigoplus_{k=1}^nc_{i,k}\times(A_j)_k=\bigoplus_{k=1}^nc_{i,k}\times f_{k,j}=f_{i,j+1} \]于是
\[C\times A_j=A_{j+1} \]但是，想要用快速幂解决这道题，我们还需要证明这个广义矩阵乘法具有结合律，即
\[A_{j+2}=C\times A_{j+1}=C\times\left(C\times A_{j}\right)=C^2\times A_j \]即
\[(C(CA))_{i,j}=\bigoplus_{k=1}^nc_{i,k}\times\left(\bigoplus_{t=1}^nc_{k,t}\times A_{t,j}\right)\\ ((CC)A)_{i,j}=\bigoplus_{t=1}^n\left(\bigoplus_{k=1}^nc_{i,k}\times c_{k,t}\right)\times A_{t,j} \]乍一看这个好像不太对，但是再多看一眼题解就会发现，\(C\) 是一个 \(0/1\) 矩阵。于是，不管两个式子中的哪一个，\(A_{t,j}\) 有贡献当且仅当 \(c_{i,k}\) 和 \(c_{k,t}\) 都是 \(1\)。那么这两个式子就是相等的了，于是就可以用快速幂解决 \(A_n=C^{n-1}\times A_1\)。
然而对于题目中的多次询问怎么办呢？我们可以一开始预处理 \(C\) 的 \(2^i\) 次幂，然后询问时把 \(a_i\) 二进制拆分，拼起来即可。
预处理复杂度 \(O(n^3\log a_i)\)，询问总复杂度 \(O(n^2q\log a_i)\)。
#include <cassert>
#include <cstring>
#include <iostream>
#define f(x, y, z) for (int x = (y); x <= (z); ++x)
#define int long long
using namespace std;
int constexpr N = 110;
int n, m, q;
int a[N];
struct Matrix {
int r, c;
int a[N][N];
Matrix(int const _r = 0, int const _c = 0): r(_r), c(_c) {}
Matrix operator*(Matrix const &B) const {
assert(c == B.r);
Matrix C(r, B.c);
f(i, 1, r) f(j, 1, B.c) {
C.a[i][j] = 0; //!!!!!!!!!!!
f(k, 1, c) C.a[i][j] ^= (a[i][k] * B.a[k][j]);
}
return C;
}
} c[33], f;
void Prework() {
f(i, 1, 31) c[i] = c[i - 1] * c[i - 1];
return;
}
void Solve() {
int x;
cin >> x;
f.r = 1, f.c = n;
memcpy(f.a[1] + 1, a + 1, sizeof(int) * n);
for (int w = 1, i = 0; w <= x; w <<= 1, ++i)
if (x & w) f = f * c[i];
cout << f.a[1][1] << '\n';
return;
}
signed main() {
ios::sync_with_stdio(false);
cin.tie(0), cout.tie(0);
cin >> n >> m >> q;
f(i, 1, n) cin >> a[i];
int u, v;
c[0].r = c[0].c = n;
f(i, 1, m) {
cin >> u >> v;
c[0].a[u][v] = c[0].a[v][u] = 1;
}
Prework();
while (q--) Solve();
return 0;
}