[KOI 2025 #2] 序列与查询

· · 题解

子问题 1

当给定一个查询时，可以固定 i 并增加 j，从而枚举所有连续子区间并计算它们的和。输出其中的最大值即可。时间复杂度为 O(QN^2)。

子问题 2

我们定义数组 S_i = \sum_{j=1}^{i-1} A_j。这可以通过递推式 S_{i+1} = S_i + A_i - X 在 O(N) 时间内计算出来。

当给定查询 X 时，问题的答案是 \max_{0 \le j < i \le N} (S_i + iX - S_j - jX) = \max_{1 \le i \le N} (S_i + iX - (\min_{0 \le j \le i-1} S_j + jX))。我们定义数组 D_i = \min_{j<i} (S_j + jX)。这可以通过递推式 D_{i+1} = \min(D_i, S_i + iX) 在 O(N) 时间内计算出来。问题的答案是 \max_{1 \le i \le N} (S_i + iX - D_{i-1})，并且可以在 O(N) 时间内计算。时间复杂度为 O(QN)。

子问题 3

我们定义 M_l 为长度为 l 的子数组中和的最大值。数组 M 可以在 O(N^2) 时间内计算出来。现在，问题的答案是 \max_{l=1}^{N} (M_l + lX)。如果直接计算，每个查询需要 O(N) 的时间，需要进行优化。

当我们在二维平面上绘制 N 个点 (x_l, y_l) = (l, M_l) 时，对于查询 P，我们想计算 x_l \cdot P + y_l 的最大值（输入中用 X 表示，但为了避免混淆，这里用 P 表示）。当我们求出这 N 个点的凸包（convex hull）时，可以发现，取得最大值的点是在凸包上画一条斜率为 -P 的直线时的切点。这样的切点可以通过求出凸包的上半部分（upper convex hull），并对每条直线使用二分搜索在 O(\log N) 时间内找到。时间复杂度为 O(N^2 + Q \log N)。

子问题 6

快速计算数组 M 是很困难的，因此需要其他方法。

让我们考虑一个与子问题 3 的算法等价但效率稍低的方法：当在二维平面上给定 N(N+1)/2 个点 \{(i - j, S_i - S_j)\}_{0 \le j < i \le N} 时，对于查询 P，问题答案就是 x_l \cdot P + y_l 的最大值。在求出所有 N(N+1)/2 个点的上凸包之后，可以使用与子问题 3 相同的二分搜索，在 O(\log N) 时间内计算每个查询。但是，由于点的数量太多，求凸包并不容易。

我们用分治法来解决这个问题。定义函数 f(L, R) 返回点集 \{(i-j, S_i - S_j)\}_{L \le j < i \le R} 的上凸包。定义 M=(L+R)/2，然后调用 f(L, M) 和 f(M+1, R)，就可以得到 i \le M, j \ge L 以及 i > M, j > M 情况下凸包的候选点。因此，我们计算点集 \{(i-j, S_i - S_j)\}_{L \le j \le M < i \le R} 的上凸包，并与递归计算出的凸包候选点合并，然后再次计算这些点的上凸包并返回。

某个点 (i-j, S_i - S_j) 在上凸包上，等价于存在一个斜率 p，使得 (i-j)p + S_i - S_j 的最大值在该点唯一取得。（换句话说，所有其他点的 (i-j)p + S_i - S_j 的值都比该点小）。为此，S_i + ip 在 i \in [M+1, R] 中必须是最大的，而 S_j + jp 在 j \in [L, M] 中必须是最小的。这等价于，我们只需要考虑在 M+1 \le i \le R 中位于上凸包上的点 (i, S_i)，以及在 L \le j \le M 中位于下凸包上的点 (j, S_j)。

当我们求出点集 \{(i, S_i)\}_{i=M+1}^R 的上凸包后，可以使用凸包上的相邻点来计算每个点取得最大值时的 p 的范围。同样地，当我们求出点集 \{(j, S_j)\}_{j=L}^M 的下凸包后，可以计算每个点取得最小值时的 p 的范围。根据上述逻辑，我们只需计算那些 p 的范围重叠的 (i, j) 对即可，这与归并排序中的合并过程类似，可以用非常短的代码来解决。（这种计算被称为 Minkowski Sum。）

分治算法由于在求凸包的过程中需要排序，因此需要 T(N) = 2T(N/2) + O(N \log N) = O(N \log^2 N) 的时间。因此，时间复杂度为 O(N \log^2 N + Q \log N)。

子问题 7

上述分治算法可以在 O(N \log N) 的时间内实现。首先，在求点集 \{(i, S_i)\}_{i=M+1}^R 的上凸包时，可以自然地按 i 递增的顺序处理各点，然后使用 Monotone Chain 算法在 O(N) 时间内计算上凸包。下凸包也同样。此外，没有必要让 f(L, R) 精确地返回上凸包：在合并过程中得到的点可以不按任何顺序收集起来，最后再统一计算一次上凸包即可。最终，在处理查询之前计算上凸包时，由于每个点的 x 坐标都在 [1, N] 范围内，因此可以用将每个 x 坐标对应的 y 坐标最大值记录在数组中的方式来代替排序。时间复杂度为 O((N+Q) \log N)。

Remark 1. 由于输入限制较大，在使用 CCW 计算凸包等时，可能会发生 long long 范围溢出的问题。这个问题有四种不同的解决方法（按便利性递增的顺序列出）。示例代码是用第二种方法实现的。

Remark 2. 关于这个问题，先前已知的解法是 Sakai 在 2024 年的论文 [Sak24]，其预处理时间为 O(N \log^2 N)。本题解的解法在保持相同的空间复杂度和查询时间的同时，将预处理时间改进为 O(N \log N)。此外，与论文中的方法相比，本题解的解法要简单得多。

References

[Sak24] Yoshifumi Sakai. A data structure for the maximum-sum segment problem with offsets. In 35th Annual Symposium on Combinatorial Pattern Matching (CPM 2024), pages 26-1. Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2024.

#include <bits/stdc++.h>
using namespace std;
using lint = long long;
using pi = array<lint, 2>;
#define sz(v) ((int)(v).size())
#define all(v) (v).begin(), (v).end()
#define cr(v, n) (v).clear(), (v).resize(n);

pi operator+(const pi &a, const pi &b) { return pi{a[0] + b[0], a[1] + b[1]}; }
pi operator-(const pi &a, const pi &b) { return pi{a[0] - b[0], a[1] - b[1]}; }
int ccw(pi a, pi b) {
    double apx = 1.0 * a[0] * b[1] - 1.0 * a[1] * b[0];
    if (fabs(apx) > 1e18)
        return apx > 0 ? 1 : -1;
    lint ans = a[0] * b[1] - a[1] * b[0];
    if (ans == 0)
        return 0;
    return ans > 0 ? 1 : -1;
}
int ccw(pi a, pi b, pi c) { return ccw(b - a, c - a); }

struct cht {
    vector<pi> ch;
    bool bad(pi a, pi b, pi c) { return ccw(b - a, b - c) <= 0; }
    void init(vector<pi> &a) {
        for (auto &x : a) {
            if (sz(ch) && ch.back()[0] == x[0])
                continue;
            while (sz(ch) >= 2 && bad(ch[sz(ch) - 2], ch.back(), x))
                ch.pop_back();
            ch.push_back(x);
        }
    }
    lint query(lint x) {
        int s = 0, e = sz(ch) - 1;
        while (s != e) {
            int m = (s + e) / 2;
            if (ch[m][0] * x + ch[m][1] > ch[m + 1][0] * x + ch[m + 1][1])
                e = m;
            else
                s = m + 1;
        }
        return ch[s][0] * x + ch[s][1];
    }
} cht;

vector<lint> a, dap;

void dnc(int l, int r) {
    if (l == r)
        return;
    int m = (l + r) / 2;
    dnc(l, m);
    dnc(m + 1, r);
    vector<pi> L, R;
    for (int j = m + 1; j <= r; j++) {
        pi p{j, a[j]};
        while (sz(R) >= 2 && ccw(R[sz(R) - 2], R.back(), p) >= 0)
            R.pop_back();
        R.push_back(p);
    }
    for (int j = m; j >= l; j--) {
        pi p{-j, -a[j]};
        while (sz(L) >= 2 && ccw(L[sz(L) - 2], L.back(), p) >= 0)
            L.pop_back();
        L.push_back(p);
    }
    pi beg = L[0] + R[0];
    vector<pi> d1, d2;
    for (int k = 0; k + 1 < sz(L); k++) {
        d1.push_back(L[k + 1] - L[k]);
    }
    for (int k = 0; k + 1 < sz(R); k++) {
        d2.push_back(R[k + 1] - R[k]);
    }
    vector<pi> delta(sz(d1) + sz(d2));
    merge(all(d1), all(d2), delta.begin(), [&](const pi &a, const pi &b) { return ccw(a, b) < 0; });
    for (int i = 0; i < sz(delta); i++) {
        dap[beg[0] - 1] = max(dap[beg[0] - 1], beg[1]);
        beg = beg + delta[i];
    }
    dap[beg[0] - 1] = max(dap[beg[0] - 1], beg[1]);
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout.tie(0);
    int n, q;
    cin >> n >> q;
    cr(a, n + 1);
    cr(dap, n);
    for (int i = 0; i < n; i++) {
        lint z;
        cin >> z;
        a[i + 1] = a[i] + z;
    }
    fill(all(dap), -1e18);
    dnc(0, n);
    vector<pi> points;
    for (int i = 0; i < n; i++) {
        points.push_back({i + 1, dap[i]});
    }
    cht.init(points);
    for (int i = 0; i < q; i++) {
        lint x;
        cin >> x;
        cout << cht.query(x) << "\n";
    }
}