[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 \to y\)，显然可以写出下面的转移方程
\[dp[x][t] = \min( c[i]+\sum_{k=1}^Tdp[y][t+k] \times p[i][k])
\]
\]
注意到把\(dp[y]\)反转成\(dp'[y]\)后,乘积可以化为\(\sum_{k=1}^Tdp'[y][2T-(t+k)] \times p[i][k])\),这是一个卷积的形式,可以用分治FFT优化
由于图上有环，按照点的顺序转移是有后效性的。但是如果我们按照时间的顺序分治,先计算出所有点中\(t\)的范围在\([mid+1,r]\)中的\(dp\)值,用FFT计算出它们的代价\(c[i]+\sum_{k=1}^Tdp[y][t+k] \times p[i][k]\),用这个代价来更新\([l,mid]\)的dp值,就可以避免后效性。因为时间小的状态一定由时间大的状态转移过来。
初始值:
\[dp[i][j]= \begin{cases} 0,i=n,j\in[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\n",dp[1][0]);
}