@description@

一个 n 点 m 边的有向轨道系统，每条边有一个费用，Kyoya Ootori 从点 1 出发去点 n。
由于这个轨道系统不太完美，所以每条边消耗的时间是[1, t]中的某个整数。现在已知第 i 条边消耗时间 j 的概率。如果 Kyoya Ootori 在严格大于 t 时刻到达点 n，则他要付罚款 x。
Kyoya Ootori 可以根据到达某个点的时间调整他的策略，求 Kyoya Ootori 在最优策略下，到达点 n 最小期望费用。

input
第一行 4 个整数 n, m, t, x（2 <= n <= 50， 1 <= m <= 100，1 <= t <= 20000，0 <= x <= 10^6）。
接下来 2m 行描述了边的信息：
第 2i 行三个整数 ai bi ci，表示有一条费用为 ci 的有向边由 ai 到 bi。保证一定存在从 1 到 n 的合法路径。
第 2i + 1 行包含 t 个整数，pi1,pi2,...,pit，表示第 i 条边消耗时间 j 的概率为 pij/10^5。保证 pi1 + pi2 + ... + pit = 10^5。
保证不存在单方向的重边。

output
输出最小期望费用。

sample input
4 4 5 1
1 2 0
50000 0 50000 0 0
2 3 0
10000 0 0 0 90000
3 4 0
100000 0 0 0 0
2 4 0
0 0 0 50000 50000
sample output
0.7000000000

@solution@

@part - 1@

有一个很 naive 的 dp：
令 \(dp[u][i]\) 表示在 i 时刻到达点 u，从点 u 出发到达点 n 的最小期望费用。
当 \(i > t\) 时，反正都会罚款了，直接破罐子破摔走最短路。即：
\[dp[u][i]=dist(u, n)+x\]
当 \(i <= t\) 时，则有转移：
\[dp[u][i] = \min\{\sum_{j=1}^{j<=t}(dp[v][i+j]\times p[e][j])+dist[e]\}\]
转移要求图中存在 (u, v) 这条有向边，其中 e 表示 (u, v) 这条边。
时间复杂度 \(O(m*t*t)\)，显然没这么简单。

@part - 2@

【注意下文中的变量 i, j, k 的意义一直在变化】

观察转移式，貌似是个卷积。
令 \(f[e][i] = \sum_{j=1}^{j<=t}(dp[v][i+j]\times p[e][j])\)，则 \(dp[u][i] = \min\{f[e][i]+dist[e]\}\)。

现在再来考虑怎么转换 \(f\) 的求解。
再对 \(f\) 进行变形，得到：\(f[e][i] = \sum_{j-k=i}(dp[v][j]\times p[e][k])\)。
这……难不成是一个减法卷积？

不要慌我们换个元。
令 \(dp'[v][j] = dp[v][2t-j]\)，上式变为 \(f[e][i] = \sum_{2t-j-k=i}(dp'[v][j]\times p[e][k])\)。
好像还是不对劲？

不要慌我们再换个元。令 \(f'[e][i] = f[e][2t-i]\)，则：
\[f'[e][i] = f[e][2t-i] = \sum_{2t-j-k=2t-i}(dp'[v][j]\times p[e][k]) = \sum_{j+k=i}(dp'[v][j]\times p[e][k])\]
此时 \(dp'[u][i] = dp[u][2t-i]=\min\{f[e][2t-i]+dist[e]\}=\min\{f'[e][i]+dist[e]\}\)。

感性理解一下：我们上面做的这些相当于把 \(f\) 和 \(dp\) 翻转了一下变为 \(f'\) 与 \(dp'\)。

这里为什么不把 \(p\) 变成 \(p'\) 呢？等会儿就知道了 qwq。

至少它终于变成了我们熟悉的卷积形式。

@part - 3@

然而问题又来了。
我们并不能一次性 get 到某个点所有的 dp 值。因为假如原图中有环，这个点的 dp 值可以更新环上其他点的 dp 值，反过来又回来更新这个点的其他 dp 值。
所以一般的 FFT 并不能搞。

但是，这个问题和 dp 本质的无后效性又有区别。
感性解释一下：这个 dp 有两维，第一维是点的编号，是不具有无后效性的，所以不能直接 用卷积优化。但是第二维是时间，是具有无后效性的，所以可以用 dp 。

怎么办呢？我们用 cdq 分治套 FFT 来解决这一问题（好像学名叫作 “分治 FFT”）。
考虑抽象出这样一个模型：

\(C\) 是 \(A,B\) 的卷积，其中 \(B\) 是一开始就已知的。依次给出 \(a_0, a_1,\dots\) 的值，当给出 \(a_i\) 的值时，需要立即算出 \(c_i\) 的值。

在这个题目中，\(f'[e]\) 是 \(dp'[v]\) 与 \(p[e]\) 的卷积，其中 \(p[e]\) 是一开始已知的，\(dp'[v]\) 是一个个知道的。
再来看我们上面所提到的那个问题，为什么不把 \(p\) 换成 \(p'\)。现在答案就很明显了，因为 \(dp\) 是从后往前依次知道的，所以我们必须要把 \(dp\) 换成 \(dp'\)，才能运用接下来的方法。

假如我们已知 \([le, mid]\) 内所有的 \(A[i]\)，则 \([le, mid]\) 对 \([mid+1, ri]\) 的贡献为：
\[C[i] += \sum_{le \le j\le mid}^{j+k=i}A[j]*B[k]（mid+1 \le i \le ri）\]
可以求出 k 的范围为 1 <= k <= ri-le。这是一个长度为 ri - le 的卷积，可以用 FFT 来优化。

cdq 分治共有 \(\log t\) 层，每一层的 FFT 时间复杂度为 \(O(t\log t)\)，一共有 m 条边，故要进行 \(O(m)\) 次 FFT。所以总时间复杂度 \(O(mt\log^2t)\)

@accepted code@

#include<cmath>
#include<vector>
#include<cstdio>
#include<algorithm>
using namespace std;
const int MAXN = 50;
const int MAXM = 100;
const int MAXT = 10*20000;
const int INF = int(1E9);
const double PI = acos(-1);
struct complex{double r, i;complex(double _r=0, double _i=0):r(_r), i(_i){}
}tmp1[MAXT + 5], tmp2[MAXT + 5], tmp3[MAXT + 5];
typedef complex cplx;
cplx operator + (cplx a, cplx b) {return cplx(a.r+b.r, a.i+b.i);}
cplx operator - (cplx a, cplx b) {return cplx(a.r-b.r, a.i-b.i);}
cplx operator * (cplx a, cplx b) {return cplx(a.r*b.r-a.i*b.i, a.i*b.r+a.r*b.i);}
cplx w[2][MAXT + 5];
void init() {int s; for(s = 2;s <= MAXT;s <<= 1); s >>= 1;w[0][s] = cplx(cos(2*PI/s), sin(2*PI/s));w[1][s] = cplx(w[0][s].r, -w[0][s].i);for(s >>= 1;s >= 2;s >>= 1) {w[0][s] = w[0][s<<1] * w[0][s<<1];w[1][s] = w[1][s<<1] * w[1][s<<1];}
}
void fft(cplx *A, int n, int type) {for(int i=0,j=0;i<n;i++) {if( i > j ) swap(A[i], A[j]);for(int l=(n>>1);(j^=l)<l;l>>=1);}for(int s=2;s<=n;s<<=1) {int t = (s>>1); cplx u = w[type][s];for(int i=0;i<n;i+=s) {cplx p = cplx(1, 0);for(int j=0;j<t;j++, p=p*u) {cplx x = A[i+j], y = A[i+j+t];A[i+j] = x + p*y; A[i+j+t] = x - p*y;}}}if( type ) {for(int i=0;i<n;i++)A[i].r /= n;}
}
struct edge{int u, v, w;double p[MAXT + 5], f[MAXT + 5];
}e[MAXM + 5];
int G[MAXN + 5][MAXN + 5];
int n, m, t, x;
void floyd() {for(int k=1;k<=n;k++)for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)G[i][j] = min(G[i][j], G[i][k] + G[k][j]);
}
double dp[MAXN + 5][MAXT + 5];
void divide_conquer(int le, int ri) {if( le == ri ) {if( le >= t )for(int i=1;i<=m;i++)dp[e[i].u][le] = min(dp[e[i].u][le], e[i].f[le] + e[i].w);return ;}int mid = (le + ri) >> 1, len;divide_conquer(le, mid);for(len = 1;len < ri-le;len <<= 1);for(int i=1;i<=m;i++) {for(int j=le;j<=mid;j++)tmp1[j-le] = cplx(dp[e[i].v][j], 0);for(int j=mid-le+1;j<len;j++)tmp1[j] = cplx(0, 0);for(int j=1;j<=ri-le;j++)tmp2[j-1] = cplx(e[i].p[j], 0);for(int j=ri-le;j<len;j++)tmp2[j] = cplx(0, 0);fft(tmp1, len, 0); fft(tmp2, len, 0);for(int j=0;j<len;j++)tmp3[j] = tmp1[j]*tmp2[j];fft(tmp3, len, 1);for(int j=mid+1;j<=ri;j++)e[i].f[j] += tmp3[j-le-1].r;}divide_conquer(mid+1, ri);
}
int main() {init();scanf("%d%d%d%d", &n, &m, &t, &x);for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)if( i != j ) G[i][j] = INF;for(int i=1;i<=m;i++) {scanf("%d%d%d", &e[i].u, &e[i].v, &e[i].w);G[e[i].u][e[i].v] = min(G[e[i].u][e[i].v], e[i].w);for(int j=1;j<=t;j++) {int p; scanf("%d", &p);e[i].p[j] = 1.0*p/100000;}}floyd();for(int i=1;i<=n;i++) {for(int j=0;j<t;j++)dp[i][j] = G[i][n] + x;for(int j=t;j<=2*t;j++)if( i == n ) dp[i][j] = 0;else dp[i][j] = INF;}divide_conquer(0, 2*t);printf("%lf\n", dp[1][2*t]);
}

@details@

边写博客边理清思路，写完博客再写代码的我 www。
注意下标的转化，变来变去的很容易脑袋一晕就什么也看不出来了。