[Codeforces 553E]Kyoya and Train(期望DP+Floyd+分治FFT)

[Codeforces 553E]Kyoya and Train(期望DP+Floyd+分治FFT)

题面

给出一个(n)个点(m)条边的有向图(可能有环),走每条边需要支付一个价格(c_i),需要的时间为([1,T])中随机的整数,时间为(j)的概率为(p_{i,j})。从(1)出发走到(n),如果到(n)的时间超过(T),就需要再支付(X)。找出一条路径,使得支付钱数的期望值最小。输出最小期望。

(n leq 50,m leq 100,T leq 20000)

分析

设(dp[x][t])表示当前走到了节点为(x),已经经过的时间为(t)时，从(x)到终点的最小花费。那么我们可以枚举从(x)出发的每条边(i:x o y)，显然可以写出下面的转移方程

[dp[x][t] = min( c[i]+sum_{k=1}^Tdp[y][t+k] imes p[i][k]) ]

注意到把(dp[y])反转成(dp'[y])后,乘积可以化为(sum_{k=1}^Tdp'[y][2T-(t+k)] imes p[i][k])),这是一个卷积的形式,可以用分治FFT优化

由于图上有环，按照点的顺序转移是有后效性的。但是如果我们按照时间的顺序分治,先计算出所有点中(t)的范围在([mid+1,r])中的(dp)值,用FFT计算出它们的代价(c[i]+sum_{k=1}^Tdp[y][t+k] imes p[i][k]),用这个代价来更新([l,mid])的dp值,就可以避免后效性。因为时间小的状态一定由时间大的状态转移过来。

初始值:

[dp[i][j]= egin{cases} 0,i=n,jin[1,T] \ + infty,i in[1,n-1],j in[1,T] \ X+mincost[i][n] ,j in[T+1,2T]end{cases} ]

最后一种情况中(mincost[i][n])表示(i)到(n)的最小代价,可以用Floyd预处理。这里实际上用了贪心的思想,既然现在已经迟到了,那之后无论走多快都会被罚(X)。不如选价格最小的路慢慢走。由于所有([T+1,2T])的dp值已经求出,可以做一次FFT计算代价,然后直接对([0,T])做分治。

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#define maxn 50
#define maxm 100
#define maxt 20000
#define INF 1e20
const double pi=acos(-1.0);
using namespace std;
struct com {
	double real;
	double imag;
	com() {

	}
	com(double _real,double _imag) {
		real=_real;
		imag=_imag;
	}
	com(double x) {
		real=x;
		imag=0;
	}
	void operator = (const com x) {
		this->real=x.real;
		this->imag=x.imag;
	}
	void operator = (const double x) {
		this->real=x;
		this->imag=0;
	}
	friend com operator + (com p,com q) {
		return com(p.real+q.real,p.imag+q.imag);
	}
	friend com operator + (com p,double q) {
		return com(p.real+q,p.imag);
	}
	void operator += (com q) {
		*this=*this+q;
	}
	void operator += (double q) {
		*this=*this+q;
	}
	friend com operator - (com p,com q) {
		return com(p.real-q.real,p.imag-q.imag);
	}
	friend com operator - (com p,double q) {
		return com(p.real-q,p.imag);
	}
	void operator -= (com q) {
		*this=*this-q;
	}
	void operator -= (double q) {
		*this=*this-q;
	}
	friend com operator * (com p,com q) {
		return com(p.real*q.real-p.imag*q.imag,p.real*q.imag+p.imag*q.real);
	}
	friend com operator * (com p,double q) {
		return com(p.real*q,p.imag*q);
	}
	void operator *= (com q) {
		*this=(*this)*q;
	}
	void operator *= (double q) {
		*this=(*this)*q;
	}
	friend com operator / (com p,double q) {
		return com(p.real/q,p.imag/q);
	}
	void operator /= (double q) {
		*this=(*this)/q;
	}
	void print() {
		printf("%lf + %lf i ",real,imag);
	}
};


void fft(com *x,int *rev,int n,int type) {//orz Fourier
	for(int i=0; i<n; i++) if(i<rev[i]) swap(x[i],x[rev[i]]);
	for(int len=1; len<n; len*=2) {
		int sz=len*2;
		com wn1=com(cos(2*pi/sz),type*sin(2*pi/sz));
		for(int l=0; l<n; l+=sz) {
			int r=l+len-1;
			com wnk=1;
			for(int i=l; i<=r; i++) {
				com tmp=x[i+len];
				x[i+len]=x[i]-wnk*tmp;
				x[i]=x[i]+wnk*tmp;
				wnk=wnk*wn1;
			}
		}
	}
	if(type==-1) for(int i=0; i<n; i++) x[i]/=n;
}

int n,m,T,X;
struct edge {
	int from;
	int to;
	int cost;
	double p[maxt+5];
} E[maxm+5];
int dist[maxn+5][maxn+5];
void floyd() {
	for(int k=1; k<=n; k++) {
		for(int i=1; i<=n; i++) {
			for(int j=1; j<=n; j++) {
				dist[i][j]=min(dist[i][j],dist[i][k]+dist[k][j]);
			}
		}
	}
}

double sum[maxm+5][maxt*2+5];
double dp[maxn+5][maxt*2+5];
void calc(int l,int mid,int r) { //计算时间在[mid+1,r]的dp,对时间在[l,mid]的dp的贡献
	static com a[maxt*10+5],b[maxt*10+5];
	static int rev[maxt*10+5];
	for(int e=1; e<=m; e++) {//枚举转移的边 
		int N=1,L=0;
		int tlim=min(T,r-l+1);
		while(N<r-mid+tlim) {
			N*=2;
			L++;
		}
		for(int i=0;i<N;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
		for(int i=0;i<N;i++){
			a[i]=0;
			b[i]=0;
		}
		//计算[mid+1,r]转移时的代价，是卷积形式 
		for(int i=0;i<r-mid;i++){
			a[i]=dp[E[e].to][r-i];//反转
		}
		for(int i=0;i<=tlim;i++){
			b[i]=E[e].p[i];
		}
		fft(a,rev,N,1);
		fft(b,rev,N,1);
		for(int i=0;i<N;i++) a[i]*=b[i];
		fft(a,rev,N,-1);
		//把答案累计到[l,mid]上 
		for(int i=l;i<=mid;i++) sum[e][i]+=a[r-i].real;
	}
}
void divide(int l,int r){//orz CDQ 
	if(l==r){
		//所有需要的数已经计算完毕,可以转移了
		for(int i=1;i<=m;i++){
			dp[E[i].from][l]=min(dp[E[i].from][l],sum[i][l]+E[i].cost);
		} 
		return;
	}
	int mid=(l+r)>>1;
	divide(mid+1,r);
	calc(l,mid,r);
	divide(l,mid);
}

int main(){
	scanf("%d %d %d %d",&n,&m,&T,&X);
	memset(dist,0x3f,sizeof(dist));
	for(int i=1;i<=m;i++){
		scanf("%d %d %d",&E[i].from,&E[i].to,&E[i].cost);
		dist[E[i].from][E[i].to]=min(dist[E[i].from][E[i].to],E[i].cost);
		for(int j=1;j<=T;j++){
			scanf("%lf",&E[i].p[j]); 
			E[i].p[j]/=1e5;
		} 
	}
	for(int i=1;i<=n;i++) dist[i][i]=0;
	floyd();
	for(int i=1;i<=n;i++){
		for(int j=0;j<=T;j++){
			//初始化不超时的花费 
			if(i!=n) dp[i][j]=INF;
			else dp[i][j]=0;
			//如果超时,就随便选距离最小的去 
			dp[i][j+T]=X+dist[i][n];
		}
	}
	calc(0,T,2*T+1);
	divide(0,T);
	printf("%.6f
",dp[1][0]);
}