@bzoj - 4386@ [POI2015] Wycieczki
@description@
给定一张 n 个点 m 条边的带权有向图，每条边的边权只可能是1，2，3中的一种。
将所有可能的路径按路径长度排序，请输出第 k 小的路径的长度，注意路径不一定是简单路径，即可以重复走同一个点。
input
第一行包含三个整数 n, m, k (1<=n<=40，1<=m<=1000，1<=k<=10^18)。
接下来 m 行，每行三个整数 u, v, c(1<=u,v<=n，u不等于v，1<=c<=3)，表示从 u 出发有一条到 v 的单向边，边长为 c 。
可能有重边。
output
包含一行一个正整数，即第 k 短的路径的长度，如果不存在，输出-1。
sample input
6 6 11
1 2 1
2 3 2
3 4 2
4 5 1
5 3 1
4 6 3
sample output
4
sample explain
长度为1的路径有1->2，5->3，4->5。
长度为2的路径有2->3，3->4，4->5->3。
长度为3的路径有4->6，1->2->3，3->4->5，5->3->4。
长度为4的路径有5->3->4->5。
@solution@
一种比较直观的想法是二分答案。二分出来 mid 后我们相当于是要计算路径长度 <= mid 的路径总数。
考虑怎么求路径长度 = x 的路径数量。讨论路径中最后一条边的边权是 1, 2, 3 中的哪一个，可以由路径长度为 x - 1, x - 2, x - 3 的状态转移过来。
我们令 G1 表示这样一个矩阵，其中 G1[i][j] 为 i->j 边权为 1 的边数量。G2, G3 同理。
再定义 F[k] 表示这样一系列矩阵，其中 F[k][i][j] 为 i 到 j 路径长度为 k 的路径总数。
则有：
这样，我们计算的路径长度 <= mid 的路径总数，就等于 (F[1] + F[2] + ... + F[mid]) 这样一个矩阵中所有元素之和。
观察上面那个转移式，它是一个线性递推式。我们不妨令 F[0] = I（单位矩阵），F[-1] = F[-2] = O（零矩阵）。我们可以把这个转移式写成矩阵嵌套矩阵的形式，，可以得到：
其中 S[k] 是一个统计元素和的，1*n 的向量。矩阵中的 i 是 1*n 的全 1 向量，o1 是 1*n 的全 0 向量，o2 是 n*1 的全 0 向量。
初始向量\(\begin{bmatrix}
S[0]\\
F[0]\\
F[-1]\\
F[-2]\\
\end{bmatrix}\)，最后将 S[mid] 的所有元素之和求出来与 k 比较大小。
这样我们就得到了一个 O(n3log2 k) 的算法，但并不够优秀。
考虑类比树上求某个点的第 k 个祖先的做法，即使用倍增。
记矩阵 A0 = \(\begin{bmatrix}
1 & i & o1 & o1\\
o2 & G1 & G2 & G3\\
o2 & I & O & O\\
o2 & O & I & O\\
\end{bmatrix}\)，Ai = A0(2i)。
对 A 使用倍增即可。时间复杂度 O(n^3log k)。
最后一个问题：路径数是指数级增长的，它有可能溢出 long long。
没事儿，long long 不行我们用 double，double 不行我们用 long double。
@accepted code@
#include<cstdio>
#include<iostream>
using namespace std;
const int MAXN = 40 + 5;
struct matrix{
long double m[3*MAXN][3*MAXN];
int r, c;
}G[64 + 5];
int n;
long double fun(matrix A) {
long double ret = 0;
for(int i=1;i<=n;i++)
ret += A.m[0][i];
return ret;
}
matrix operator * (matrix A, matrix B) {
matrix C; C.r = A.r, C.c = B.c;
for(int i=0;i<C.r;i++)
for(int j=0;j<C.c;j++) {
C.m[i][j] = 0;
for(int k=0;k<A.c;k++)
C.m[i][j] += A.m[i][k]*B.m[k][j];
}
return C;
}
int main() {
int m; long double k;
scanf("%d%d", &n, &m);
cin >> k; k += n;
G[0].r = G[0].c = 3*n + 1;
for(int i=0;i<=n;i++) G[0].m[0][i]++;
for(int i=1;i<=n;i++) G[0].m[i+n][i]++, G[0].m[i+2*n][i+n]++;
for(int i=1;i<=m;i++) {
int u, v, c; scanf("%d%d%d", &u, &v, &c);
G[0].m[u][v+c*n-n]++;
}
int i;
for(i=1;i<=65;i++) {
G[i] = G[i-1]*G[i-1];
if( fun(G[i]) >= k ) break;
if( i == 65 ) {
puts("-1");
return 0;
}
}
i--;
matrix M = G[i];
long long ret = 1LL<<i;
for(;i>=0;i--) {
if( fun(M*G[i]) < k ) {
ret += (1LL<<i);
M = M*G[i];
}
}
printf("%lld\n", ret);
}
@details@
网上一搜，全是什么拆点之类的。
心里一惊。难道我想错了。
最后发现其实通过拆点构造出来的矩阵和用递推构造出来的矩阵是一样的。
处理 long long 溢出这么复杂，还不如直接用 long double（虽然常数要大一些）
