Problem Link

Statement

有 $n$ 个酿酒点和 $m$ 个存酒点。对于第 $i$ 个酿酒点，生产 $x$ （$x$ 为实数）升酒需要花费 $a_i \times x^2 + b_i \times x$，最多能生产 $c_i$ 升酒。对于第 $i$ 个存酒点，其最多能存储 $d_i$ 升酒。每个酿酒点生产酒后需立即运至该酿酒点可到达的存酒点，每个酿酒点可到达的存酒点已知并给出。求：1. 所有酿酒点最大一共能产生多少酒 2. 在酿酒量达到最大值时花费的最小值为多少，以分数形式给出。

$n, m \le 100,~0 \le a_i,b_i,c_i,d_i \le 3$

Solution

第一问建图后直接最大流即可得出。

酿酒代价为二次函数不易于计算，一个自然的想法是将这个二次函数从定义域上切割为若干段，将每一段内的部分近似为一个一次函数，然后就可以对于每一段建立一条对应的边，跑费用流即可。由于此处二次函数斜率递增，定义域靠前的段在代价上优于靠后的段，因此不需要考虑取这些段的先后顺序而可能带来的不合法取法。

但是此题要求输出分数，也就是不允许误差的存在。我们发现上面的做法分的段数越多，结果误差就越小，当每一段的长度无线接近于零，段数接近正无穷时误差即可消除。此时每一段内的一次函数就成为了该二次函数的导数：$2 a_i x + b_i$。

虽然解决了误差问题，但是段数达到了正无穷，我们无法直接对于每一段建边跑费用流。但是我们发现，若假设存酒点容量无穷，即只考虑产酒点总产量而不考虑存酒点限制的情况下，所有的最优方案一定是这样的：存在一个实数 $l$ 作为阈值，每一个产酒点对应的花费二次函数上导数小于 $l$ 的部分对应的定义域大小即为该产酒点产量。贪心策略下显然成立。

因此我们考虑构建函数 $F(x)$ 表示当阈值为 $x$ 时，不考虑存酒点限制情况下的产酒点总产量。若令 $F_i(x)$ 表示阈值为 $x$ 时第 $i$ 个产酒点的产量，可以发现 $F_i(x) = \frac 1 {2a_i} x - \frac {b_i} {2a_i} (x \in [b_i,~2 a_i c_i + b_i])$。而 $F(x) = \sum_i F_i(x)$，所以 $F(x)$ 应当由 $O(n)$ 个一次函数首尾相连形成。

该阈值与费用流 SPFA 跑出的 $dist[T]$ 功能类似，若要模拟费用流过程，我们可以令阈值从 $0$ 开始不断增长，每次更改阈值后重新通过 $F_i(x)$ 计算每条边的对应流量。费用流过程中还可将存酒点的流量限制加入，从而算出考虑存酒点限制情况下的总产量。

由于不允许存在精度误差因此我们无法做到真正模拟该过程，但是通过这个思路我们发现在阈值不断增加的过程中，存酒点总是先尽可能多地存酒，直至其容量达到上限后对部分酿酒点的产量产生限制。也就是说，我们可以通过 $F(x)$ 结合网络流得到 $G(x)$ 表示考虑存酒点容量限制的情况下，阈值为 $x$ 时的总产量。$G(x)$ 只会在 $F(x)$ 的基础上多 $O(m)$ 个端点，其图像由 $O(n + m)$ 个一次函数首尾相连形成。

但是 $G(x)$ 不像 $F(x)$ 可以直接通过计算得出，$G(x)$ 的每一个点都需要网络流计算得出。不过我们发现 $G(x)$ 的大部分图像呈现上凸壳的形态，我们可以在定义域上分治得到整个图像：令 $solve(l,~r)$ 过程查找横坐标在 $(l,~r)$ 间的所有端点位置，我们在 $l$ 和 $r$ 上分别用网络流算出坐标，并结合残量网络算出斜率，可得两个端点处的一次函数。取这两个一次函数的交点，令其为 $mid$，计算 $mid$ 处的一次函数，若与两端的一次函数相同则说明 $(l,~r)$ 范围内已再无端点，否则向 $(l,~mid)$ 和 $(mid,~r)$ 递归求解。但是注意到 $a_i$ 可能等于零，也就是说一些一次函数可能出现斜率无穷的情况，因此需要在 $x = 0,~1,~2,~3$ 处强制断开，在 $(0,~1),~(1,~2),~(2,~3),~(3,~+\infty)$ 段分开分治求解。

计算出 $G(x)$ 的图像后，通过计算其积分即可得到最终答案。

Code

View on GitHub

/**
 * @file 2138.cpp
 * @author Macesuted (i@macesuted.moe)
 * @date 2022-07-05
 *
 * @copyright Copyright (c) 2022
 *
 */

#include <bits/stdc++.h>
using namespace std;

#ifndef LOCAL
#define endl '\n'
#endif

bool mem1;

#define maxn 105
#define eps 1e-9
#define feps 1e-11

bool con[maxn][maxn];
int n, m, a[maxn], b[maxn], c[maxn], d[maxn];

class Frac {
   private:
    int64_t x, y;

    int64_t gcd(int64_t a, int64_t b) { return b ? gcd(b, a % b) : a; }

   public:
    Frac(int64_t x_ = 0, int64_t y_ = 1) {
        int64_t g = gcd(x_, y_);
        x = x_ / g, y = y_ / g;
    }
    operator double() { return (double)x / y; }
    operator int64_t() { return x / y; }
    Frac operator+(const Frac& o) const { return Frac(x * o.y + y * o.x, y * o.y); }
    Frac operator-(const Frac& o) const { return Frac(x * o.y - y * o.x, y * o.y); }
    Frac operator*(const Frac& o) const { return Frac(x * o.x, y * o.y); }
    Frac operator/(const Frac& o) const { return Frac(x * o.y, y * o.x); }
    bool operator==(const Frac& o) const { return x == o.x && y == o.y; }
    Frac& operator+=(const Frac& o) { return *this = *this + o; }
    Frac& operator-=(const Frac& o) { return *this = *this - o; }
    Frac& operator*=(const Frac& o) { return *this = *this * o; }
    Frac& operator/=(const Frac& o) { return *this = *this / o; }
    void print(void) {
        if (y < 0) x = -x, y = -x;
        return cout << x << '/' << y, void();
    }
};

typedef pair<Frac, Frac> pff;

class Dinic {
   private:
    struct Edge {
        int to, rev;
        double cap;
    };

    vector<vector<Edge>> graph;
    vector<vector<Edge>::iterator> cur;
    vector<int> dist;
    queue<int> que;

    bool bfs(void) {
        for (int i = 1; i <= n + m + 2; i++) dist[i] = -1, cur[i] = graph[i].begin();
        que.push(n + m + 1), dist[n + m + 1] = 0;
        while (!que.empty()) {
            int p = que.front();
            que.pop();
            for (auto i : graph[p])
                if (i.cap > feps && !~dist[i.to]) dist[i.to] = dist[p] + 1, que.push(i.to);
        }
        return ~dist[n + m + 2];
    }
    double dfs(int p, double rest) {
        if (rest < feps || p == n + m + 2) return rest;
        double use = 0, c;
        for (auto i = cur[p]; i != graph[p].end() && rest > feps; i++) {
            cur[p] = i;
            if (dist[p] + 1 != dist[i->to] || i->cap < feps) continue;
            if ((c = dfs(i->to, min(rest, i->cap))) < feps) dist[i->to] = -1;
            i->cap -= c, graph[i->to][i->rev].cap += c, use += c, rest -= c;
        }
        return use;
    }

   public:
    void addEdge(int from, int to, double cap) {
        return graph[from].push_back(Edge{to, (int)graph[to].size(), cap}),
               graph[to].push_back(Edge{from, (int)graph[from].size() - 1, 0});
    }
    pff dinic(double lim) {
        graph.clear(), graph.resize(n + m + 3), cur.resize(n + m + 3), dist.resize(n + m + 3);
        for (int i = 1; i <= n; i++)
            if (a[i] == 0) {
                if (b[i] < lim) addEdge(n + m + 1, i, c[i]);
            } else
                addEdge(n + m + 1, i, min((lim - b[i]) / 2. / a[i], (double)c[i]));
        for (int i = 1; i <= n; i++)
            for (int j = 1; j <= m; j++)
                if (con[i][j]) addEdge(i, n + j, 1e18);
        for (int i = 1; i <= m; i++) addEdge(n + i, n + m + 2, d[i]);
        while (bfs()) dfs(n + m + 1, 1e18);
        Frac K, B;
        for (int i = 1; i <= n; i++)
            if (dist[i] == -1) {
                if (a[i] == 0) {
                    if (b[i] < lim) B += Frac(c[i]);
                } else if (lim > 2 * a[i] * c[i] + b[i])
                    B += Frac(c[i]);
                else if (lim > b[i])
                    K += Frac(1, 2 * a[i]), B -= Frac(b[i], a[i] * 2);
            }
        for (int i = 1; i <= m; i++)
            if (~dist[n + i]) B += Frac(d[i]);
        return {K, B};
    }
} net;

vector<pff> nodes;

void solve(pff l, pff r) {
    if (l == r) return;
    Frac mid = (l.second - r.second) / (r.first - l.first);
    pff fml = net.dinic((double)mid - eps), fmr = net.dinic((double)mid + eps);
    if (l == fml) return nodes.emplace_back(mid, fml.first * mid + fml.second), void();
    return solve(l, fml), solve(fmr, r);
}

void solve(void) {
    cin >> n >> m;
    for (int i = 1; i <= n; i++) cin >> a[i] >> b[i] >> c[i];
    for (int i = 1; i <= m; i++) cin >> d[i];
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= m; j++) cin >> con[i][j];
    cout << int64_t(net.dinic(1e18).second) << endl;
    nodes.emplace_back(Frac(), Frac());
    for (int i = 1; i <= 3; i++) {
        auto l = net.dinic(i - 1 + eps), r = net.dinic(i - eps);
        solve(l, r);
        nodes.emplace_back(Frac(i), Frac(i) * r.first + r.second);
    }
    solve(net.dinic(3 + eps), net.dinic(1e18));
    Frac ans;
    for (int i = 1; i < (int)nodes.size(); i++) {
        pff f1 = net.dinic((double)nodes[i - 1].first + eps), f2 = net.dinic((double)nodes[i].first - eps),
            f3 = net.dinic((double)nodes[i].first + eps);
        ans += Frac(1, 2) * (nodes[i].second - f1.first * nodes[i - 1].first - f1.second) *
                   (nodes[i - 1].first + nodes[i].first) +
               nodes[i].first * (f3.first * nodes[i].first + f3.second - f2.first * nodes[i].first - f2.second);
    }
    ans.print(), cout << endl;
    return;
}

bool mem2;

int main() {
    ios::sync_with_stdio(false), cin.tie(nullptr);
#ifdef LOCAL
    cerr << "Memory Cost: " << abs(&mem1 - &mem2) / 1024. / 1024. << "MB" << endl;
#endif

    int _ = 1;
    while (_--) solve();

#ifdef LOCAL
    cerr << "Time Cost: " << clock() * 1000. / CLOCKS_PER_SEC << "MS" << endl;
#endif
    return 0;
}