【Codeforces553E_CF553E】Kyoya and Train（概率_CDQ分治_FFT）
题目
我为什么要写这道题？因为说到 553 ，你有没有想到 ……
翻译
这个 Kyoya Ootori 怎么看都像是日语名字但是我是真查不出来对应的汉字是什么（好像是什么京屋鳳之类的），方便起见直接认为主人公叫张三。
题目名称：张三和火车
描述
张三想坐火车去学校。有 \(n\) 个火车站和在不同车站间开行的 \(m\) 条单向火车线路。张三现在在 \(1\) 号火车站，学校在 \(n\) 号火车站。他必须买票才能坐火车，并且坐火车也需要花费时间。然而，由于火车存在缺陷，所以到达目的地所需的时间是随机的。如果张三到达学校的时间严格大于 \(t\) ，他需要支付 \(x\) 元的罚款。
对每条火车线路会给出票价和花费时间的概率分布。更加形式化地，第 \(i\) 条线路的票价是 \(c_i\) ；对于 \(1\leq k\leq t\) ，概率分布 \(p_{i,k}\) 表示这条线路需要花费 \(k\) 单位时间的概率。张三搭乘火车所需的时间是互相独立的随机变量（此外，如果张三搭乘一条火车线路超过一次，这条线路可能花费不同的时间，这些值也互相独立）。
张三想让去学校的花费的期望最小（包括票价和迟到的罚款）。当然，张三会采取最优方案去学校，并且每当他到达一个火车站，他可以根据目前已经花费的时间重新计算他的最优方案。如果张三采取最优方案，他去学校的花费的期望是多少？
输入
第一行包含四个整数 \(n,m,t,x(2\leq n\leq 50,1\leq m\leq 100,1\leq t\leq 20000,0\leq x\leq 10^6)\) 。
接下来 \(2m\) 行描述 \(m\) 条火车线路。
第 \(2i\) 行包含三个整数 \(a_i,b_i,c_i\) ，表示第 \(i\) 条单向线路从 \(a_i\) 到 \(b_i\) ，票价为 \(c_i(1\leq a_i,b_i\leq n,a_i\neq b_i,0\leq c_i\leq 10^6)\) 。
第 \((2i+1)\) 行包含 \(t\) 个整数 \(p_{i,1},p_{i,2},\dots,p_{i,t}\) ，\(\frac{p_{i,k}}{100000}\) 是这条线路花费 \(k\) 单位时间的概率。 \((\) 对于 \(1\leq k\leq t,0\leq p_{i,k}\leq 100000,\sum_{k=1}^t p_{i,k}=100000)\) 。
保证任意一对车站之间每个方向只有不超过一条线路。
输出
输出一个实数，表示到学校花费的最小期望。和标准答案的绝对或相对误差不超过 \(10^{-6}\) 即为正确。
分析
其实思路挺显然的。
很明显，如果已经迟到了，因为迟到一分钟和迟到一天没区别，那么就不需要管时间了，直接沿着花钱最少的路线走到学校就可以了。因此只需要考虑还没有迟到时，即时间不超过 \(t\) 的决策。因为时间 \(a\) 可以转移到时间 \(a+t\) ，所以为了方便，实际上要考虑的范围是 \([0,2t]\) 。
看到 \(n\cdot t\) 这么小，考虑用 \(f_{i,j}\) 表示当前在 \(i\) ，已经过去了 \(j\) 单位时间。然后就有一个很显然的转移方程：
（其中边 \(e\) 起点为 \(u\) ，终点为 \(v\) ）
虽然图不是 DAG ，但只会从只会从时间晚的状态转移到时间早的状态，所以转移没有环。直接暴力做的复杂度是 \(O(mt^2)\) 。
仔细观察一下， \(\sum_{k=1}^tf_{v,i+k}\cdot p_{e,k}\) 是不是长得一脸多项式卷积的样子（什么？你说 \(i+k\) 和 \(k\) 是差为定值不是和为定值？你把其中一个翻转一下不就是和为定值了吗）。但很不幸的一点是，这个式子里用到了 \(f\) ，也就是卷积的结果。这是一个经典问题，详见 【洛谷4721】【模板】分治FFT（CDQ分治_NTT） 。
具体到这道题来说。因为 \(f\) 不是和式，而是对一堆和式取 min ，所以要设 \(g_{e,i}\) 表示 \(e\) 这条边对应的 \(\sum_{k=1}^tf_{v,i+k}\cdot p_{e,k}\) 。
对于每个分治的区间，先向下分治 \((mid,r]\) ，再枚举每一条边，用 FFT 计算 \((mid,r]\) 里的 \(f\) 对 \([l,mid]\) 里的 \(g\) 的贡献，最后向下分治 \([l,mid]\) 。因此，当分治到 \([l,l]\) 时，对 \(l\) 有贡献的区间（即 \((l,+\infin)\) ）已经全部算完了，所以直接用 \(g_{e,l}\) 更新 \(f_{u,l}\) 。
用 FFT 计算 \((mid,r]\) 对 \([l,mid]\) 的贡献时，因为先分治了右区间，所以已经知道了 \((mid,r]\) 的 \(f\) 。把整个区间翻转，这样就变成了传统的已知左边计算右边的模型，同时也把上面式子里的 \(i+k\) 变成了 \(i-k\) 。
我感觉我说了一团浆糊。其实思路是比较好理解的，上面两段没说清楚的话直接看代码吧。
分治的每一层的多项式长度之和是 \(O(T)\) 的，因此时间复杂度是 \(O(mT\log^2T)\) 。
代码
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cctype>
#include <vector>
#include <functional>
#include <queue>
#include <cmath>
using namespace std;
namespace zyt
{
template<typename T>
inline bool read(T &x)
{
char c;
bool f = false;
x = 0;
do
c = getchar();
while (c != EOF && c != '-' && !isdigit(c));
if (c == EOF)
return false;
if (c == '-')
f = true, c = getchar();
do
x = x * 10 + c - '0', c = getchar();
while (isdigit(c));
if (f)
x = -x;
return true;
}
inline bool read(double &x)
{
return ~scanf("%lf", &x);
}
template<typename T>
inline void write(T x)
{
char buf[20];
char *pos = buf;
if (x < 0)
putchar('-'), x = -x;
do
*pos++ = x % 10 + '0';
while (x /= 10);
while (pos > buf)
putchar(*--pos);
}
inline void write(const double x, const int fixed = 9)
{
printf("%.*f", fixed, x);
}
typedef pair<int, int> pii;
const int N = 60, M = 110, T = 2e4 + 10, B = 20, INF = 0x3f3f3f3f;
int n, m, t, x, dis[N], head[N], ecnt;
struct edge
{
int from, to, w, next;
}e[M];
double p[M][T], f[N][T << 1], g[M][T << 1];
void add(const int a, const int b, const int c)
{
e[ecnt] = (edge){a, b, c, head[a]}, head[a] = ecnt++;
}
void Dijkstra()
{
static bool vis[N];
priority_queue<pii, vector<pii>, greater<pii> > q;
memset(dis, INF, sizeof(int[n + 1]));
memset(vis, 0, sizeof(bool[n + 1]));
dis[n] = 0;
q.push(pii(0, n));
while (!q.empty())
{
int u = q.top().second;
q.pop();
if (vis[u])
continue;
vis[u] = true;
for (int i = head[u]; ~i; i = e[i].next)
{
int v = e[i].to;
if (dis[v] > dis[u] + e[i].w)
dis[v] = dis[u] + e[i].w, q.push(pii(dis[v], v));
}
}
}
namespace Polynomial
{
const int LEN = T * 8;
const double PI = acos(-1);
struct cpx
{
double x, y;
cpx(const double _x = 0, const double _y = 0)
: x(_x), y(_y) {}
cpx operator + (const cpx &b) const
{
return cpx(x + b.x, y + b.y);
}
cpx operator - (const cpx &b) const
{
return cpx(x - b.x, y - b.y);
}
cpx operator * (const cpx &b) const
{
return cpx(x * b.x - y * b.y, x * b.y + y * b.x);
}
};
cpx omega[B][LEN], winv[B][LEN];
int rev[B][LEN];
bool flag[B];
void init(const int n, const int lg2)
{
if (flag[lg2])
return;
flag[lg2] = true;
cpx w = cpx(cos(2 * PI / n), sin(2 * PI / n)), wi = cpx(cos(2 * PI / n), -sin(2 * PI / n));
omega[lg2][0] = winv[lg2][0] = cpx(1, 0);
for (int i = 1; i < n; i++)
{
omega[lg2][i] = omega[lg2][i - 1] * w;
winv[lg2][i] = winv[lg2][i - 1] * wi;
}
for (int i = 0; i < n; i++)
rev[lg2][i] = ((rev[lg2][i >> 1] >> 1) | ((i & 1) << (lg2 - 1)));
}
void fft(cpx *const a, const cpx *const w, const int *const rev, const int n)
{
for (int i = 0; i < n; i++)
if (i < rev[i])
swap(a[i], a[rev[i]]);
for (int l = 1; l < n; l <<= 1)
for (int i = 0; i < n; i += (l << 1))
for (int k = 0; k < l; k++)
{
cpx x = a[i + k], y = a[i + l + k] * w[n / (l << 1) * k];
a[i + k] = x + y;
a[i + l + k] = x - y;
}
}
void mul(const double *const a, const double *const b, double *const c, const int n)
{
static cpx x[LEN], y[LEN];
int m = 1, lg2 = 0;
while (m < (n + n - 1))
m <<= 1, ++lg2;
init(m, lg2);
for (int i = 0; i < n; i++)
x[i] = a[i], y[i] = b[i];
memset(x + n, 0, sizeof(cpx[m - n]));
memset(y + n, 0, sizeof(cpx[m - n]));
fft(x, omega[lg2], rev[lg2], m), fft(y, omega[lg2], rev[lg2], m);
for (int i = 0; i < m; i++)
x[i] = x[i] * y[i];
fft(x, winv[lg2], rev[lg2], m);
for (int i = 0; i < n; i++)
c[i] = x[i].x / m;
}
}
void solve(const int l, const int r)
{
if (l > t)
return;
if (l == r)
{
for (int i = 0; i < ecnt; i++)
f[e[i].to][l] = min(f[e[i].to][l], g[i][l] + e[i].w);
return;
}
int mid = (l + r) >> 1, len = r - l + 1;
solve(mid + 1, r);
for (int i = 0; i < ecnt; i++)
{
static double a[T << 1], b[T << 1];
a[0] = 0;
for (int j = 0; j < len; j++)
{
a[j] = (l + len - j - 1 > mid ? f[e[i].from][l + len - j - 1] : 0);
b[j] = (j <= t ? p[i][j] : 0);
}
Polynomial::mul(a, b, a, len);
for (int j = 0; j <= mid - l; j++)
g[i][j + l] += a[len - j - 1];
}
solve(l, mid);
}
int work()
{
read(n), read(m), read(t), read(x);
memset(head, -1, sizeof(int[n + 1]));
for (int i = 0; i < m; i++)
{
int a, b, c;
read(a), read(b), read(c);
add(b, a, c);
for (int j = 1; j <= t; j++)
read(p[i][j]), p[i][j] /= 1e5;
}
Dijkstra();
for (int i = 1; i <= n; i++)
{
for (int j = 0; j <= t; j++)
f[i][j] = (i == n ? 0 : INFINITY);
for (int j = t + 1; j <= (t << 1); j++)
f[i][j] = dis[i] + x;
}
solve(0, t << 1);
write(f[1][0]);
return 0;
}
}
int main()
{
return zyt::work();
}
