【CDQ分治】【FFT】【DP】【最短路】CodeForces553E Kyoya and Train

最新推荐文章于 2022-05-22 09:29:56 发布

616156

最新推荐文章于 2022-05-22 09:29:56 发布

阅读量274

点赞数

分类专栏： DP CDQ分治 FFT 多项式最短路

本文链接：https://blog.csdn.net/qq_34454069/article/details/85242683

版权

DP 同时被 3 个专栏收录

100 篇文章 1 订阅

订阅专栏

多项式

22 篇文章 0 订阅

订阅专栏

FFT

9 篇文章 0 订阅

订阅专栏

分析：

首先，很容易想到一个会T的DP：
定义 $f (i, j)$ 表示在i号点，用了j单位时间，到达目的地的最小期望代价。

转移很显然： $f(i,j)=min\{Cost_{i->v}+\sum f(v,j+k)*P_{i->v,k}\}$

这个式子显然可以看做是一个卷积的形式，但如果直接卷起来显然要出问题：
这个转移式如果用FFT优化，那么首先必须得到所有 $f (v, j + k)$ 的最终数值。并且当前的 $f (i, j)$ 还处于一个都没算出来的阶段。在不是DAG的图中，这样就有可能出现问题：即某个 $f (i, j)$ （未算出）的状态会转移到某个 $f (v, j + k)$ （已算出）的状态中，换言之，这个DP就失去了无后效性。

但这并不意味着就不能FFT优化了。

很显然，由于 $P$ 数组的最小下标为1，即走一条边必定会消耗至少1的时间。那么计算 $f (i, j)$ 时，只需要保证所有 $f (v, j^{'}) [j^{'} > j]$ 的状态，即所有时间维度比它小的都算出来就可以了。

那么，就可以借助CDQ分治的思想（其实和差不多就是CDQ）

假设当前要处理时间维度在 $[l, r]$ 区间内的状态。

那么可以先处理 $[m i d + 1, r]$ 的答案，求出其值后，就可以用它们来贡献 $[l, m i d]$ 的答案。然后再计算 $[l, m i d]$ 内部贡献的答案。

所以，这题麻烦之处就在于：每次FFT都是用于处理 $[l, m i d]$ 与 $[m i d + 1, r]$ 这两部分的卷积，由于每一层卷积的长度不同，所以下标转换起来比较恶心。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
#include<cmath>
#define SF scanf
#define PF printf
#define MAXN 105
#define MAXP 55
#define MAXM 200010
#define INF 0x3FFFFFFF
using namespace std;
const double Pi=acos(-1);
struct Edge{
	int u,v,c;
}edge[MAXN];
struct cpx{
	double r,i;	
	cpx() {}
	cpx(double r1,double i1):r(r1),i(i1) {}
	cpx operator + (const cpx &a) const {
		return cpx(r+a.r,i+a.i);	
	}
	cpx operator - (const cpx &a) const {
		return cpx(r-a.r,i-a.i);	
	}
	cpx operator * (const cpx &a) const {
		return cpx(r*a.r-i*a.i,r*a.i+i*a.r);
	}
}A[MAXM],B[MAXM];
double p[MAXN][MAXM],f[MAXP][MAXM],S[MAXN][MAXM];
int dist[MAXN][MAXN];
int n,m,t,x;
void FFT(cpx a[],int N,int flag){
	int i,j,k;
	for(i=1,j=0;i<N;i++){
		for(int d=N;j^=d>>=1,~j&d;);
		if(i<j)
			swap(a[i],a[j]);
	}
	for(i=1;i<N;i<<=1){
		cpx wn=cpx(cos(Pi/i),flag*sin(Pi/i));
		for(j=0;j<N;j+=(i<<1)){
			cpx w=cpx(1,0);
			for(k=0;k<i;k++,w=w*wn){
				cpx X=a[j+k],Y=w*a[i+j+k];
				a[j+k]=X+Y;
				a[i+j+k]=X-Y;	
			}
		}
	}
	if(flag==-1)
		for(int i=0;i<N;i++)
			a[i].r/=N;
}
void cal(int l,int r,int mid){
	int len=r-l;
	for(int j=1;j<=m;j++){
		int v=edge[j].v;
		int N=1;
		while(N<len+r-mid+1)
			N<<=1;
		for(int i=0;i<N;i++)
			A[i]=B[i]=cpx(0,0);
		for(int i=0;i<=r-mid;i++)
			A[i]=cpx(f[v][r-i],0);
		for(int i=0;i<len;i++)
			B[i]=cpx(p[j][i+1],0);
		FFT(A,N,1);
		FFT(B,N,1);
		for(int i=0;i<N;i++)
			A[i]=A[i]*B[i];
		FFT(A,N,-1);
		for(int i=l;i<mid;i++)
			S[j][i]+=A[r-i-1].r;
	}
}
void CDQ(int l,int r){
	if(l==r){
		for(int i=1;i<=m;i++){
			int u=edge[i].u;
			f[u][l]=min(f[u][l],S[i][l]+edge[i].c);
		}
		return ;
	}	
	int mid=(l+r)>>1;
	CDQ(mid+1,r);
	cal(l,r,mid+1);
	CDQ(l,mid);
}
int main(){
	SF("%d%d%d%d",&n,&m,&t,&x);
	int u,v,len;
	memset(dist,0x3f,sizeof dist);
	for(int i=1;i<=n;i++)
		dist[i][i]=0;
	for(int i=1;i<=m;i++){
		SF("%d%d%d",&u,&v,&len);	
		edge[i].u=u;
		edge[i].v=v;
		edge[i].c=len;
		dist[u][v]=min(dist[u][v],len);
		for(int j=1;j<=t;j++){
			SF("%lf",&p[i][j]);
			p[i][j]/=100000.0;
		}
	}
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			for(int k=1;k<=n;k++)
				dist[j][k]=min(dist[j][k],dist[j][i]+dist[i][k]);
//	for(int i=1;i<=n;i++)
//		PF("[%d %d]\n",i,dist[i][n]);
	for(int i=0;i<=t;i++)
		f[n][i]=0;
	for(int i=t+1;i<=2*t;i++)
		f[n][i]=x;
	for(int i=1;i<n;i++)
		for(int j=0;j<=2*t;j++){
			if(j>t)
				f[i][j]=dist[i][n]+x;
			else
				f[i][j]=INF;
		}
	cal(1,t*2,t+1);//update from the "illegal situation"
	CDQ(0,t);
	PF("%lf",f[1][0]);
}