cf553 E Kyoya and Train （分治FFT）

题意

给一张有向图，每条边的通过时间是一个1~t的给定离散概率分布。每条边还有一个通过费用。你现在要从1号点到n号点，假如你到达n号点的时间比t大，那么你需要额外付出X的费用。请求出最优策略下的最小期望费用。
$n\le 50,m\le 100,t\le 20000$

分析

• 先搞清楚最优策略是什么意思：
• 设$f[t][x]$表示时刻t在x时，到达n的最小期望费用。
• 根据定义，则有$f[t][x]=min(c_e+{\sum }_{a}f[t+a][y]\times P_a)$
• $f[t^{\prime }][x],t^{\prime }>t$的值：此时必定付出额外费用，因此是走花费最短路，可以通过预处理得出。
• 这个dp可以通过分治fft来优化。 对于t’>t的情况我是先通过预处理算好贡献的。
• 具体来说，对时间分治：现在要计算$t\in [mid+1,R]$对$t\in [L,mid]$的状态的贡献。枚举每一条边，可以发现这个转移就是一个减法卷积。做一个FFT就可以了。
• 实现上，要注意将次数界控制在32768，否则有被卡常的风险。

#include <bits/stdc++.h>
using namespace std;
typedef double db;
const int N = 3.5e4;
const db E = 1e5, inf = 1e9, pi = acos(-1);
const int SZ = sizeof E;
int n,m,t,W;
struct edge{
	int from, to, c;
	db p[N];
} e[110];
int d[55][55];
db upd[20010][110];
db f[20010][55];
db A[N], B[N], C[N];
complex<db> w[N];
int M, h[N], MM = 32768;

void fft(complex<db> *a) {
	for(int i = 0; i < M; i++) 
		if (h[i] > i) swap(a[i], a[h[i]]);
	int st = MM;
	for(int m = 1; m < M; m <<= 1) {
		st >>= 1;
		for(int i = 0; i < M; i += (m << 1)) {
			for(int j = 0; j < m; j++) {
				complex<db> z = a[i + j + m] * w[j * st];
				a[i + j + m] = a[i + j] - z;
				a[i + j] = a[i + j] + z;
			}
		}
	}
}

void ifft(complex<db> *a) {
	reverse(a + 1, a + M);
	fft(a);
	db z = 1.0 / M;
	for(int i = 0; i < M; i++) a[i].real(a[i].real() * z);
}

int mx;
void calc(int sz,int cc) {
	int sa = sz - cc;
	memset(C, 0, (sz + 5) * SZ);
	for(M = 1; M <= sa + sz; M <<= 1);
	static complex<db> a[N], b[N];
	for(int i = 1; i < M; i++) 
		h[i] = (h[i >> 1] >> 1) + (i & 1) * (M >> 1);
	memset(a, 0, M * sizeof a[0]);
	memset(b, 0, M * sizeof b[0]);
	for(int i = 0; i <= sz; i++) {
		a[i].real(A[i]);
		b[i].real(B[sz - i]);
	}
	fft(a);
	fft(b);
	for(int i = 0; i < M; i++) a[i] = a[i] * b[i];
	ifft(a);
	for(int i = 0; i < sz; i++) {
		C[i] = a[i + sz - cc].real();
	}
}

void divide(int L, int R) {
	if (L > R) return;
	if (L == R) {
		for(int i = 1; i <= n; i++) f[L][i] = inf;
		for(int i = 1; i <= m; i++) {
			f[L][e[i].from] = min(f[L][e[i].from], upd[L][i] + e[i].c);
		}
		f[L][n] = 0;
		return;
	}

	int mid = L + R >> 1, sz = R - L + 1;
	divide(mid + 1, R);

	for(int i = 1; i <= m; i++) {
		int x = e[i].from, y = e[i].to;
		memset(A, 0, (sz + 5) * SZ);
		memset(B, 0, (sz + 5) * SZ);
		int cc = mid - L + 1;
		for(int j = mid + 1; j <= R; j++) {
			A[j - L - cc] = f[j][y];
		}
		for(int j = 1; j <= R - L + 1; j++) {
			B[j] = e[i].p[j];
		}
		calc(sz, cc);
		for(int j = L; j <= mid; j++) {
			upd[j][i] += C[j - L];
		}
	}
	divide(L, mid);
}

void read(int &a) {
	char c; while((c=getchar()) > '9' || c < '0');
	a = c - '0'; while ((c=getchar()) >= '0' && c <= '9') a = a * 10 + c - '0';
}

int main() {
	int g = clock();
	freopen("e.in","r",stdin);
	cin >> n >> m >> t >> W;
	memset(d,127/2,sizeof d);
	db WW = 1.0 / E;
	for(int i = 1; i <= m; i++) {
		scanf("%d %d %d", &e[i].from, &e[i].to, &e[i].c);
		d[e[i].from][e[i].to] = min(d[e[i].from][e[i].to], e[i].c);
		for(int j = 1; j <= t; j++) {
			int z; read(z);
			e[i].p[j] = z * WW;
		}
	}
	for(int i = 1; i <= n; i++) d[i][i] = 0;
	for(int k = 1; k <= n; k++) {
		for(int i = 1; i <= n; i++) if (i != k) {
			for(int j = 1; j <= n; j++) if (j != k && j != i) {
				d[i][j] = min(d[i][j], d[i][k] + d[k][j]);
			}
		}
	}
	for(int i = 1; i <= m; i++) {
		db po = 0;
		for(int j = 1; j <= t; j++) {
			po += e[i].p[t - j + 1];
			upd[j][i] = (d[e[i].to][n] + W) * po;
			if (upd[j][i] > inf) upd[j][i] = inf;
		}
	}
	w[0].real(1);
	for(int i = 1; i < MM; i++) {
		w[i] = complex<db>(cos(i*2*pi/MM),sin(i*2*pi/MM));
	} 
	divide(0, t);
	printf("%.10lf\n",f[0][1]);
	cerr << clock() - g << endl;
}