Codeforces553E_CF553EKyoya and Train（概率_CDQ分治_FFT）

Posted 2021-03-04 zyt1253679098

题目

Codeforces 553E

我为什么要写这道题？因为说到 553 ，你有没有想到 ……

翻译

这个 Kyoya Ootori 怎么看都像是日语名字但是我是真查不出来对应的汉字是什么（好像是什么京屋鳳之类的），方便起见直接认为主人公叫张三。

题目名称：张三和火车

描述

张三想坐火车去学校。有 (n) 个火车站和在不同车站间开行的 (m) 条单向火车线路。张三现在在 (1) 号火车站，学校在 (n) 号火车站。他必须买票才能坐火车，并且坐火车也需要花费时间。然而，由于火车存在缺陷，所以到达目的地所需的时间是随机的。如果张三到达学校的时间严格大于 (t) ，他需要支付 (x) 元的罚款。

对每条火车线路会给出票价和花费时间的概率分布。更加形式化地，第 (i) 条线路的票价是 (c_i) ；对于 (1leq kleq t) ，概率分布 (p_{i,k}) 表示这条线路需要花费 (k) 单位时间的概率。张三搭乘火车所需的时间是互相独立的随机变量（此外，如果张三搭乘一条火车线路超过一次，这条线路可能花费不同的时间，这些值也互相独立）。

张三想让去学校的花费的期望最小（包括票价和迟到的罚款）。当然，张三会采取最优方案去学校，并且每当他到达一个火车站，他可以根据目前已经花费的时间重新计算他的最优方案。如果张三采取最优方案，他去学校的花费的期望是多少？

输入

第一行包含四个整数 (n,m,t,x(2leq nleq 50,1leq mleq 100,1leq tleq 20000,0leq xleq 10^6)) 。

接下来 (2m) 行描述 (m) 条火车线路。

第 (2i) 行包含三个整数 (a_i,b_i,c_i) ，表示第 (i) 条单向线路从 (a_i) 到 (b_i) ，票价为 (c_i(1leq a_i,b_ileq n,a_i eq b_i,0leq c_ileq 10^6)) 。

第 ((2i+1)) 行包含 (t) 个整数 (p_{i,1},p_{i,2},dots,p_{i,t}) ，(frac{p_{i,k}}{100000}) 是这条线路花费 (k) 单位时间的概率。 (() 对于 (1leq kleq t,0leq p_{i,k}leq 100000,sum_{k=1}^t p_{i,k}=100000)) 。

保证任意一对车站之间每个方向只有不超过一条线路。

输出

输出一个实数，表示到学校花费的最小期望。和标准答案的绝对或相对误差不超过 (10^6) 即为正确。

分析

其实思路挺显然的。

很明显，如果已经迟到了，因为迟到一分钟和迟到一天没区别，那么就不需要管时间了，直接沿着花钱最少的路线走到学校就可以了。因此只需要考虑还没有迟到时，即时间不超过 (t) 的决策。因为时间 (a) 可以转移到时间 (a+t) ，所以为了方便，实际上要考虑的范围是 ([0,2t]) 。

看到 (ncdot t) 这么小，考虑用 (f_{i,j}) 表示当前在 (i) ，已经过去了 (j) 单位时间。然后就有一个很显然的转移方程：

[f_{u,i}=egin{cases}d_u+x(i>t) min({a|a=c_e+sum_{k=1}^{t}f_{v,i+k}cdot p_{e,k},mathrm{from}_e=u,mathrm{to}_e=v}) mathrm{otherwise}end{cases}]

（其中边 (e) 起点为 (u) ，终点为 (v) ）

虽然图不是 DAG ，但只会从只会从时间晚的状态转移到时间早的状态，所以转移没有环。直接暴力做的复杂度是 (O(mt^2)) 。

仔细观察一下， (sum_{k=1}^tf_{v,i+k}cdot p_{e,k}) 是不是长得一脸多项式卷积的样子（什么？你说 (i+k) 和 (k) 是差为定值不是和为定值？你把其中一个翻转一下不就是和为定值了吗）。但很不幸的一点是，这个式子里用到了 (f) ，也就是卷积的结果。这是一个经典问题，详见 【洛谷4721】【模板】分治FFT（CDQ分治_NTT） 。

具体到这道题来说。因为 (f) 不是和式，而是对一堆和式取 min ，所以要设 (g_{e,i}) 表示 (e) 这条边对应的 (sum_{k=1}^tf_{v,i+k}cdot p_{e,k}) 。

对于每个分治的区间，先向下分治 ((mid,r]) ，再枚举每一条边，用 FFT 计算 ((mid,r]) 里的 (f) 对 ([l,mid]) 里的 (g) 的贡献，最后向下分治 ([l,mid]) 。因此，当分治到 ([l,l]) 时，对 (l) 有贡献的区间（即 ((l,+infin)) ）已经全部算完了，所以直接用 (g_{e,l}) 更新 (f_{u,l}) 。

用 FFT 计算 ((mid,r]) 对 ([l,mid]) 的贡献时，因为先分治了右区间，所以已经知道了 ((mid,r]) 的 (f) 。把整个区间翻转，这样就变成了传统的已知左边计算右边的模型，同时也把上面式子里的 (i+k) 变成了 (i-k) 。

我感觉我说了一团浆糊。其实思路是比较好理解的，上面两段没说清楚的话直接看代码吧。

分治的每一层的多项式长度之和是 (O(T)) 的，因此时间复杂度是 (O(mTlog^2T)) 。

代码

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cctype>
#include <vector>
#include <functional>
#include <queue>
#include <cmath>
using namespace std;

namespace zyt
{
    template<typename T>
    inline bool read(T &x)
    {
        char c;
        bool f = false;
        x = 0;
        do
            c = getchar();
        while (c != EOF && c != '-' && !isdigit(c));
        if (c == EOF)
            return false;
        if (c == '-')
            f = true, c = getchar();
        do
            x = x * 10 + c - '0', c = getchar();
        while (isdigit(c));
        if (f)
            x = -x;
        return true;
    }
    inline bool read(double &x)
    {
        return ~scanf("%lf", &x);
    }
    template<typename T>
    inline void write(T x)
    {
        char buf[20];
        char *pos = buf;
        if (x < 0)
            putchar('-'), x = -x;
        do
            *pos++ = x % 10 + '0';
        while (x /= 10);
        while (pos > buf)
            putchar(*--pos);
    }
    inline void write(const double x, const int fixed = 9)
    {
        printf("%.*f", fixed, x);
    }
    typedef pair<int, int> pii;
    const int N = 60, M = 110, T = 2e4 + 10, B = 20, INF = 0x3f3f3f3f;
    int n, m, t, x, dis[N], head[N], ecnt;
    struct edge
    {
        int from, to, w, next;
    }e[M];
    double p[M][T], f[N][T << 1], g[M][T << 1];
    void add(const int a, const int b, const int c)
    {
        e[ecnt] = (edge){a, b, c, head[a]}, head[a] = ecnt++;
    }
    void Dijkstra()
    {
        static bool vis[N];
        priority_queue<pii, vector<pii>, greater<pii> > q;
        memset(dis, INF, sizeof(int[n + 1]));
        memset(vis, 0, sizeof(bool[n + 1]));
        dis[n] = 0;
        q.push(pii(0, n));
        while (!q.empty())
        {
            int u = q.top().second;
            q.pop();
            if (vis[u])
                continue;
            vis[u] = true;
            for (int i = head[u]; ~i; i = e[i].next)
            {
                int v = e[i].to;
                if (dis[v] > dis[u] + e[i].w)
                    dis[v] = dis[u] + e[i].w, q.push(pii(dis[v], v));
            }
        }
    }
    namespace Polynomial
    {
        const int LEN = T * 8;
        const double PI = acos(-1);
        struct cpx
        {
            double x, y;
            cpx(const double _x = 0, const double _y = 0)
                : x(_x), y(_y) {}
            cpx operator + (const cpx &b) const
            {
                return cpx(x + b.x, y + b.y);
            }
            cpx operator - (const cpx &b) const
            {
                return cpx(x - b.x, y - b.y);
            }
            cpx operator * (const cpx &b) const
            {
                return cpx(x * b.x - y * b.y, x * b.y + y * b.x);
            }
        };
        cpx omega[B][LEN], winv[B][LEN];
        int rev[B][LEN];
        bool flag[B];
        void init(const int n, const int lg2)
        {
            if (flag[lg2])
                return;
            flag[lg2] = true;
            cpx w = cpx(cos(2 * PI / n), sin(2 * PI / n)), wi = cpx(cos(2 * PI / n), -sin(2 * PI / n));
            omega[lg2][0] = winv[lg2][0] = cpx(1, 0);
            for (int i = 1; i < n; i++)
            {
                omega[lg2][i] = omega[lg2][i - 1] * w;
                winv[lg2][i] = winv[lg2][i - 1] * wi;
            }
            for (int i = 0; i < n; i++)
                rev[lg2][i] = ((rev[lg2][i >> 1] >> 1) | ((i & 1) << (lg2 - 1)));
        }
        void fft(cpx *const a, const cpx *const w, const int *const rev, const int n)
        {
            for (int i = 0; i < n; i++)
                if (i < rev[i])
                    swap(a[i], a[rev[i]]);
            for (int l = 1; l < n; l <<= 1)
                for (int i = 0; i < n; i += (l << 1))
                    for (int k = 0; k < l; k++)
                    {
                        cpx x = a[i + k], y = a[i + l + k] * w[n / (l << 1) * k];
                        a[i + k] = x + y;
                        a[i + l + k] = x - y;
                    }
        }
        void mul(const double *const a, const double *const b, double *const c, const int n)
        {
            static cpx x[LEN], y[LEN];
            int m = 1, lg2 = 0;
            while (m < (n + n - 1))
                m <<= 1, ++lg2;
            init(m, lg2);
            for (int i = 0; i < n; i++)
                x[i] = a[i], y[i] = b[i];
            memset(x + n, 0, sizeof(cpx[m - n]));
            memset(y + n, 0, sizeof(cpx[m - n]));
            fft(x, omega[lg2], rev[lg2], m), fft(y, omega[lg2], rev[lg2], m);
            for (int i = 0; i < m; i++)
                x[i] = x[i] * y[i];
            fft(x, winv[lg2], rev[lg2], m);
            for (int i = 0; i < n; i++)
                c[i] = x[i].x / m;
        }
    }
    void solve(const int l, const int r)
    {
        if (l > t)
            return;
        if (l == r)
        {
            for (int i = 0; i < ecnt; i++)
                f[e[i].to][l] = min(f[e[i].to][l], g[i][l] + e[i].w);
            return;
        }
        int mid = (l + r) >> 1, len = r - l + 1;
        solve(mid + 1, r);
        for (int i = 0; i < ecnt; i++)
        {
            static double a[T << 1], b[T << 1];
            a[0] = 0;
            for (int j = 0; j < len; j++)
            {
                a[j] = (l + len - j - 1 > mid ? f[e[i].from][l + len - j - 1] : 0);
                b[j] = (j <= t ? p[i][j] : 0);
            }
            Polynomial::mul(a, b, a, len);
            for (int j = 0; j <= mid - l; j++)
                g[i][j + l] += a[len - j - 1];
        }
        solve(l, mid);
    }
    int work()
    {
        read(n), read(m), read(t), read(x);
        memset(head, -1, sizeof(int[n + 1]));
        for (int i = 0; i < m; i++)
        {
            int a, b, c;
            read(a), read(b), read(c);
            add(b, a, c);
            for (int j = 1; j <= t; j++)
                read(p[i][j]), p[i][j] /= 1e5;
        }
        Dijkstra();
        for (int i = 1; i <= n; i++)
        {
            for (int j = 0; j <= t; j++)
                f[i][j] = (i == n ? 0 : INFINITY);
            for (int j = t + 1; j <= (t << 1); j++)
                f[i][j] = dis[i] + x;
        }
        solve(0, t << 1);
        write(f[1][0]);
        return 0;
    }
}
int main()
{
    return zyt::work();
}