@codeforces

目录

• @description@
• @solution@
  • @part - 1@
  • @part - 2@
  • @part - 3@
• @accepted code@
• @details@

---

@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] imes 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] imes p[e][j]))，则 (dp[u][i] = min{f[e][i]+dist[e]})。

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

不要慌我们换个元。
令 (dp'[v][j] = dp[v][2t-j])，上式变为 (f[e][i] = sum_{2t-j-k=i}(dp'[v][j] imes 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] imes p[e][k]) = sum_{j+k=i}(dp'[v][j] imes 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 jle 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(tlog t))，一共有 m 条边，故要进行 (O(m)) 次 FFT。所以总时间复杂度 (O(mtlog^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
", dp[1][2*t]);
}

@details@

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