题面

给出一个\(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]);
}