【Codeforces553E_CF553E】Kyoya and Train（概率_CDQ分治_FFT）

题目

Codeforces 553E

我为什么要写这道题？因为说到 553 ，你有没有想到 ……

翻译

这个 Kyoya Ootori 怎么看都像是日语名字但是我是真查不出来对应的汉字是什么（好像是什么京屋鳳之类的），方便起见直接认为主人公叫张三。

题目名称：张三和火车

描述

张三想坐火车去学校。有 (n) 个火车站和在不同车站间开行的 (m) 条单向火车线路。张三现在在 (1) 号火车站，学校在 (n) 号火车站。他必须买票才能坐火车，并且坐火车也需要花费时间。然而，由于火车存在缺陷，所以到达目的地所需的时间是随机的。如果张三到达学校的时间严格大于 (t) ，他需要支付 (x) 元的罚款。

对每条火车线路会给出票价和花费时间的概率分布。更加形式化地，第 (i) 条线路的票价是 (c_i) ；对于 (1leq kleq t) ，概率分布 (p_{i,k}) 表示这条线路需要花费 (k) 单位时间的概率。张三搭乘火车所需的时间是互相独立的随机变量（此外，如果张三搭乘一条火车线路超过一次，这条线路可能花费不同的时间，这些值也互相独立）。

张三想让去学校的花费的期望最小（包括票价和迟到的罚款）。当然，张三会采取最优方案去学校，并且每当他到达一个火车站，他可以根据目前已经花费的时间重新计算他的最优方案。如果张三采取最优方案，他去学校的花费的期望是多少？

输入

第一行包含四个整数 (n,m,t,x(2leq nleq 50,1leq mleq 100,1leq tleq 20000,0leq xleq 10^6)) 。

接下来 (2m) 行描述 (m) 条火车线路。

第 (2i) 行包含三个整数 (a_i,b_i,c_i) ，表示第 (i) 条单向线路从 (a_i) 到 (b_i) ，票价为 (c_i(1leq a_i,b_ileq n,a_i eq b_i,0leq c_ileq 10^6)) 。

第 ((2i+1)) 行包含 (t) 个整数 (p_{i,1},p_{i,2},dots,p_{i,t}) ，(frac{p_{i,k}}{100000}) 是这条线路花费 (k) 单位时间的概率。 (() 对于 (1leq kleq t,0leq p_{i,k}leq 100000,sum_{k=1}^t p_{i,k}=100000)) 。

保证任意一对车站之间每个方向只有不超过一条线路。

输出

输出一个实数，表示到学校花费的最小期望。和标准答案的绝对或相对误差不超过 (10^{-6}) 即为正确。

分析

其实思路挺显然的。

很明显，如果已经迟到了，因为迟到一分钟和迟到一天没区别，那么就不需要管时间了，直接沿着花钱最少的路线走到学校就可以了。因此只需要考虑还没有迟到时，即时间不超过 (t) 的决策。因为时间 (a) 可以转移到时间 (a+t) ，所以为了方便，实际上要考虑的范围是 ([0,2t]) 。

看到 (ncdot t) 这么小，考虑用 (f_{i,j}) 表示当前在 (i) ，已经过去了 (j) 单位时间。然后就有一个很显然的转移方程：

[f_{u,i}=egin{cases}d_u+x(i>t)\ min({a|a=c_e+sum_{k=1}^{t}f_{v,i+k}cdot p_{e,k},mathrm{from}_e=u,mathrm{to}_e=v}) mathrm{otherwise}end{cases} ]

（其中边 (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(mTlog^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();
}