Codeforces 553E：Kyoya and Train （最短路+概率DP+分治+FFT）

题目传送门：http://codeforces.com/contest/553/problem/E

---

题目大意：给出一幅n个点m条边的有向图，并给出参数T，你要从1号点走到n号点。经过每一条边都要花费时间和金钱，第i条边需要花费cost[i]的金钱，并且经过该边花费时间为t的概率是$p[i][t](1<=t<=T)$。如果最后你花费的总时间大于T，你就要付出额外X的金钱。请最小化期望花费的金钱。n<=50，m<=100，T<=20000。

---

题目分析：myy论文题。我们可以先想出一个很暴力的DP：f[i][t]表示在t时刻到达点i后，走完剩余路程的期望最小金钱花费。很明显边界条件为f[n][t]=0(0<=t<=T)，f[node][t]=dis[node]+X(T+1<=t)，其中dis[node]表示node到n的最短路长度。后半部分其实很容易理解，因为既然已经超过了时间限制，那么不如直接走最短路。

对于其它状态，我们假设第i条边由u到v，那么有：

$$f[u][t]=cost[i]+\sum_{j=1}^{T}p[i][j]*f[v][t+j]$$

这样时间复杂度是$O(mT^2)$的。我们观察到p[i][j]和f[v][t+j]第二维的下标之差相等，故可以将其中一个数组翻转后用FFT快速求出和式的值。又因为DP数组要按照第二维从大到小算，所以需要分治。时间复杂度为$O(mT\log^2(T))$。

说实话，这题是我一个多月前A掉的，今晚新年夜忽然想起来还没有写题解。具体细节我也忘得差不多了，大概只记得三点。一是可以另开一个g数组记录和式的值；二是FFT时下标的变换很烦；三是DP数组的边界要处理好。大概就这么多了吧，写完之后还是挺容易A的。

---

CODE：

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
using namespace std;

const int maxt=20010;
const int maxN=100000;
const int maxn=53;
const int maxm=105;
const double pi=acos(-1.0);

struct Complex
{
    double X,Y;
    Complex (double a=0.0,double b=0.0) : X(a),Y(b) {}
} ;

Complex operator+(Complex a,Complex b){return Complex(a.X+b.X,a.Y+b.Y);}
Complex operator-(Complex a,Complex b){return Complex(a.X-b.X,a.Y-b.Y);}
Complex operator*(Complex a,Complex b){return Complex(a.X*b.X-a.Y*b.Y,a.X*b.Y+a.Y*b.X);}

int Rev[maxN];
Complex A[maxN];
Complex B[maxN];
int N,Lg;

double g[maxm][maxt];
double f[maxn][maxt];

int dis[maxn];
bool vis[maxn];

double p[maxm][maxt];
double sum[maxm][maxt];

int u[maxm];
int v[maxm];
int w[maxm];
int n,m,t,cost;

void Dijkstra()
{
    for (int i=0; i<n; i++) dis[i]=1e9;
    for (int i=1; i<n; i++)
    {
        int x=0;
        for (int j=1; j<=n; j++)
            if ( !vis[j] && dis[j]<dis[x] ) x=j;
        vis[x]=true;
        for (int j=1; j<=m; j++)
            if (v[j]==x) dis[ u[j] ]=min(dis[ u[j] ],dis[x]+w[j]);
    }
}

void DFT(Complex *a,double f)
{
    for (int i=0; i<N; i++)
        if (i<Rev[i]) swap(a[i],a[ Rev[i] ]);
    for (int len=2; len<=N; len<<=1)
    {
        int mid=len>>1;
        double ang=2.0*pi/((double)len);
        Complex e( cos(ang) , f*sin(ang) );
        for (Complex *p=a; p!=a+N; p+=len)
        {
            Complex wn(1.0,0.0);
            for (int i=0; i<mid; i++)
            {
                Complex temp=wn*p[mid+i];
                p[mid+i]=p[i]-temp;
                p[i]=p[i]+temp;
                wn=wn*e;
            }
        }
    }
}

void FFT()
{
    DFT(A,1.0);
    DFT(B,1.0);
    for (int i=0; i<N; i++) A[i]=A[i]*B[i];
    DFT(A,-1.0);
    for (int i=0; i<N; i++) A[i].X/=((double)N);
    //for (int i=0; i<N; i++) if ( fabs(A[i].Y)>1e-6 ) exit(0);
}

void Solve(int L,int R)
{
    if (L==R)
    {
        for (int i=1; i<=m; i++)
            g[i][L]+=( (1.0-sum[i][t-L])*((double)(dis[ v[i] ]+cost)) ),
            f[ u[i] ][L]=min(f[ u[i] ][L],g[i][L]+(double)w[i]); //!!!!!
        return;
    }

    int mid=(L+R)>>1;
    Solve(mid+1,R);

    N=1,Lg=0;
    while (N<2*R-L-mid) N<<=1,Lg++;
    for (int i=0; i<N; i++) Rev[i]=0;
    for (int i=0; i<N; i++)
        for (int j=0; j<Lg; j++)
            if ( i&(1<<j) ) Rev[i]|=( 1<<(Lg-j-1) );

    for (int i=1; i<=m; i++)
    {
        for (int j=1; j<=R-L; j++) A[j-1]=Complex(p[i][j],0.0);
        for (int j=1; j<=R-mid; j++) B[j-1]=Complex(f[ v[i] ][mid+j],0.0);
        for (int j=0; R-mid-1-j>j; j++) swap(B[j],B[R-mid-1-j]);
        for (int j=R-L; j<N; j++) A[j]=Complex(0.0,0.0);
        for (int j=R-mid; j<N; j++) B[j]=Complex(0.0,0.0);
        FFT();
        for (int j=L; j<=mid; j++) g[i][j]+=A[R-j-1].X; //!!!!!
    }

    Solve(L,mid);
}

int main()
{
    freopen("E.in","r",stdin);
    freopen("E.out","w",stdout);

    scanf("%d%d%d%d",&n,&m,&t,&cost);
    for (int i=1; i<=m; i++)
    {
        scanf("%d%d%d",&u[i],&v[i],&w[i]);
        sum[i][0]=0.0;
        for (int j=1; j<=t; j++)
        {
            int x;
            scanf("%d",&x);
            p[i][j]=(double)x/100000.0; //!!!!!
            sum[i][j]=sum[i][j-1]+p[i][j];
        }
    }

    Dijkstra();

    for (int i=1; i<=n; i++)
        for (int j=0; j<=t; j++) f[i][j]=6e7;
    for (int j=0; j<=t; j++) f[n][j]=0.0;
    Solve(0,t);
    printf("%.10lf\n",f[1][0]);

    /*for (int i=1; i<=n; i++)
    {
        for (int j=0; j<=t; j++) printf("%.6lf ",f[i][j]);
        printf("\n");
    }
    printf("\n");
    for (int i=1; i<=m; i++)
    {
        for (int j=0; j<=t; j++) printf("%.6lf ",g[i][j]);
        printf("\n");
    }*/

    return 0;
}