题链:
http://codeforces.com/problemset/problem/623/E
题解:
FFT,DP
题意:
一个有向图,给出每条边的起点u,终点v,费用c,以及花费每种时间的概率P[e][j](表示走第e条边花费时间为j的概率)
现在需要从1号点走到n号点,如果不能在T个单位时间内到达,则达到后还要另外支付X的费用。
求出所需支付的最小期望费用。
先是一个暴力的DP方法:
(考虑到每条边的耗时至少为1,可以把状态设为类似分层图的形式)
定义$F[i][t]$为$t$时刻在$i$点时,到达终点n所需要的最小期望费用。
不难得到DP转移,枚举每条i点的出边e:
$$F[i][t]=min(c[e]+sum_{j=1}^{T}P[e][j] imes F[v[e]][t+j])$$
然后令$G[e][t]$表示$t$时刻从${u[e]}$出发到达n点所需要的最小期望费用。
即$$G[e][t]=c[e]+sum_{j=1}^{T}P[e][j] imes F[v[e]][t+j]$$
所以$$F[i][t]=min(G[e][t])quad (u[e]=i)$$
这个复杂度为$O(mT^2)$
考虑优化,注意到求G的式子有点像卷积的形式。
的确,我们只需要把$P[e]$数组翻转,即可得到:
$$G[e][t]=G'[e][T+t]=c[e]+sum_{j=1}^{T}P[e][T-j] imes F[v[e]][t+j]$$
这个就可以用FFT做了。
但是我们并不知道P数组所有的值,又怎么办呢?
这时采用分治的方法,(类似CDQ那种)
基于这样一个事实:
$$D_r=sum_{i=0}^{n-1}f_ig_{r-i}=sum_{i=0}^{k}f_ig_{r-i}+sum_{i=k+1}^{n-1}f_ig_{r-i}$$
用文字描述就是:计算卷积时,可以把贡献分成几段来分别计算,
或者说用FFT计算D的值时,可以先计算$sum_{i=0}^{k}f_ig_{r-i}这个卷积,再加上sum_{i=k+1}^{n-1}f_ig_{r-i}这个卷积$
所以我们对时间进行分治:
1.对于当前的分治层$t∈[l,r]$,我们先递归处理右边[mid,r],得到右边所有的F[i][mid~r],然后用FFT求值去贡献到G[e][l~mid-1]。
2.继续递归处理左边[l,mid-1]。
3.当到达递归最底层即$l=r$时,我们用G去求出所有的$F[i][l]$,因为此时的转移来源$G[e][t](t>l)$都计算完毕。
算法的大体流程就这样,复杂度$O(mTlog^2T)$
具体的实现还需特别注意:
1.FFT时的下标转化感觉有点绕,(我是膜的别人的做法)。
2.当t>给出的T时,由于已经超时,所以就直接按费用最短路前往n点。
3.由于每次求G[e][t]时,需要的来源状态的时间t'<=t+T,所以数组的时间那一维开为2*T即可。
代码:
#include<bits/stdc++.h> #define MAXN 55 #define MAXM 105 #define MAXT 40005 using namespace std; const double Pi=acos(-1); struct Complex{ double real,image; Complex(double _real=0,double _image=0):real(_real),image(_image){} Complex operator - () const{return Complex(-real,-image);} friend Complex operator + (const Complex &A,const Complex &B){return Complex(A.real+B.real,A.image+B.image);} friend Complex operator - (const Complex &A,const Complex &B){return A+(-B);} friend Complex operator * (const Complex &A,const Complex &B){return Complex(A.real*B.real-A.image*B.image,A.image*B.real+A.real*B.image);} }A[65538],B[65538],nul(0,0); double F[MAXN][MAXT],G[MAXM][MAXT],P[MAXM][MAXT]; int U[MAXM],V[MAXM],C[MAXM],S[MAXN][MAXN]; int order[65538]; int N,M,T,X; void read(int &x){ static int sign; static char ch; x=0; sign=1; ch=getchar(); while(ch<'0'||'9'<ch){if(ch=='-')sign=-1;ch=getchar();} while('0'<=ch&&ch<='9'){x=x*10+ch-'0';ch=getchar();} x=x*sign; } void Floyd(){ memset(S,0x3f,sizeof(S)); for(int i=1;i<=N;i++) S[i][i]=0; for(int e=1;e<=M;e++) S[U[e]][V[e]]=C[e]; 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) S[i][j]=min(S[i][j],S[i][k]+S[k][j]); } void FFT(Complex *Y,int n,int sign){ for(int i=0;i<n;i++) if(i<order[i]) swap(Y[i],Y[order[i]]); for(int d=2;d<=n;d<<=1){ Complex dw(cos(2*Pi/d),sin(sign*2*Pi/d)),w,tmp; for(int i=0;w=Complex(1,0),i<n;i+=d) for(int k=i;k<i+d/2;w=w*dw,k++) tmp=w*Y[k+d/2],Y[k+d/2]=Y[k]-tmp,Y[k]=Y[k]+tmp; } if(sign==-1) for(int i=0;i<n;i++) Y[i].real/=1.0*n; } void contribution(int l,int mid,int r){ static int m,n,len; m=r-l+1; for(n=1,len=0;n<=m+r-mid;n<<=1) len++; for(int i=1;i<n;i++) order[i]=(order[i>>1]>>1)|((i&1)<<(len-1)); for(int e=1;e<=M;e++){ for(int i=0;i<n;i++) A[i]=B[i]=nul; for(int i=1;i<m;i++) B[i]=Complex(P[e][i],0); for(int i=mid;i<=r;i++) A[r-i]=Complex(F[V[e]][i],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++) G[e][i]+=A[r-i].real; } } void solve(int l,int r){ if(l==r){ for(int i=1;i<N;i++) F[i][l]=S[i][N]+X; for(int e=1;e<=M;e++) F[U[e]][l]=min(F[U[e]][l],C[e]+G[e][l]); return; } int mid=(l+r+1)/2;//向上取整 solve(mid,r); contribution(l,mid,r); solve(l,mid-1); } int main(){ read(N); read(M); read(T); read(X); for(int e=1,x;e<=M;e++){ read(U[e]); read(V[e]); read(C[e]); for(int i=1;i<=T;i++) read(x),P[e][i]=1.0*x/100000; } Floyd(); for(int i=1;i<=N;i++) for(int j=T+1;j<=2*T;j++) F[i][j]=S[i][N]+X; for(int j=1;j<=T;j++) F[N][j]=0; contribution(0,T+1,2*T); solve(0,T); printf("%.8lf ",F[1][0]); return 0; }