【Codeforces553E_CF553E】Kyoya and Train(概率_CDQ分治_FFT)

题目

Codeforces 553E

我为什么要写这道题?因为说到 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\) 单位时间。然后就有一个很显然的转移方程:

\[f_{u,i}=\begin{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(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();
}
posted @ 2020-03-20 12:04  Inspector_Javert  阅读(155)  评论(0编辑  收藏  举报