From 000d6c64a1f6fb0d749e197f6560db97bd5f6271 Mon Sep 17 00:00:00 2001 From: hitonanode <32937551+hitonanode@users.noreply.github.com> Date: Sun, 24 Aug 2025 23:54:21 +0900 Subject: [PATCH] enhance SMAWK / add monotone minima --- other_algorithms/monotone_minima.hpp | 32 ++++++++ other_algorithms/monotone_minima.md | 33 ++++++++ other_algorithms/smawk.hpp | 78 +++++++++++++++++-- other_algorithms/smawk.md | 13 +++- .../concave_min_plus_convolution.test.cpp | 21 +++++ 5 files changed, 169 insertions(+), 8 deletions(-) create mode 100644 other_algorithms/monotone_minima.hpp create mode 100644 other_algorithms/monotone_minima.md create mode 100644 other_algorithms/test/concave_min_plus_convolution.test.cpp diff --git a/other_algorithms/monotone_minima.hpp b/other_algorithms/monotone_minima.hpp new file mode 100644 index 00000000..70de34e2 --- /dev/null +++ b/other_algorithms/monotone_minima.hpp @@ -0,0 +1,32 @@ +#pragma once +#include +#include +#include + +// finding minimum element for each row of N*M matrix +// Constraints: the solution is monotonically non-decreasing +// Complexity: O(NM logM) +// Reference: +// https://topcoder-g-hatena-ne-jp.jag-icpc.org/spaghetti_source/20120923/1348327542.html +// Verify: https://mofecoder.com/contests/monoxercon_202508/tasks/monoxercon_202508_k +template +std::vector> +MonotoneMinima(int N, int M, const std::function &weight) { + std::vector> minima(N); + + auto rec = [&](auto &&self, int il, int ir, int jl, int jr) -> void { + if (il >= ir or jl >= jr) return; + const int im = (il + ir) / 2; + T w = weight(im, jl); + int j_argmin = jl; + for (int j = jl + 1; j < jr; ++j) { + if (T wt = weight(im, j); wt < w) w = wt, j_argmin = j; + } + minima[im] = {j_argmin, w}; + self(self, il, im, jl, j_argmin + 1); + self(self, im + 1, ir, j_argmin, jr); + }; + rec(rec, 0, N, 0, M); + + return minima; +} diff --git a/other_algorithms/monotone_minima.md b/other_algorithms/monotone_minima.md new file mode 100644 index 00000000..0054db14 --- /dev/null +++ b/other_algorithms/monotone_minima.md @@ -0,0 +1,33 @@ +--- +title: Monotone minima （Monotone な行列の行最小値の効率的な探索） +documentation_of: ./monotone_minima.hpp +--- + +`SMAWK()` 関数は，monotone な $N \times M$ 行列，すなわち各行の最小値の位置が行を下るにつれ右に移動していくような行列について，各行の最小値の位置およびその値を $O((N + M) \log (N + M))$ で取得する． + +[SMAWK](./smawk.md) のライブラリと互換性があり，SMAWK が使用されている箇所は本関数で代替可能（最悪計算量のオーダーはこちらの方に log がつくが，問題によってはこちらの方が実測上速い）． + +## 使用方法 + +例えば辺重みが Monge な $N$ 頂点の DAG で頂点 $0$ から各頂点への最短路重みを求めたいとき， $N$ 行 $N - 1$ 列の行列を $(j, i)$ 成分が「直前に頂点 $i$ を経由し頂点 $j$ に到達する場合の最短路重み」であるようなものとして定めると本関数が適用できる（SMAWK も適用できるが）． + +```cpp +using T = long long; +T inf = 1LL << 60; +int N; +vector dp(N, inf); +dp[0] = 0; + +auto weight = [&](int s, int t) -> T { /* Implement */ }; + +const auto res = MonotoneMinima( + N, N - 1, [&](int j, int i) -> T { return i < j ? dp[i] + weight(i, j) : inf; }); +``` + +## 問題例 + +- [K. Coupon - Monoxer Programming Contest for Engineers](https://mofecoder.com/contests/monoxercon_202508/tasks/monoxercon_202508_k) + +## Links + +- [Totally Monotone Matrix Searching (SMAWK algorithm) - 週刊 spaghetti_source - TopCoder部](https://topcoder-g-hatena-ne-jp.jag-icpc.org/spaghetti_source/20120923/1348327542.html) diff --git a/other_algorithms/smawk.hpp b/other_algorithms/smawk.hpp index 8758f04b..50e275da 100644 --- a/other_algorithms/smawk.hpp +++ b/other_algorithms/smawk.hpp @@ -6,7 +6,7 @@ #include // SMAWK: finding minima of totally monotone function f(i, j) (0 <= i < N, 0 <= j < M) for each i -// Constraints: every submatrix of f(i, j) is monotone. +// Constraints: every submatrix of f(i, j) is monotone (= totally monotone). // Complexity: O(N + M) // Reference: // https://topcoder-g-hatena-ne-jp.jag-icpc.org/spaghetti_source/20120923/1348327542.html @@ -46,8 +46,19 @@ std::vector> SMAWK(int N, int M, const std::function +std::vector> +SMAWKAntiMonotone(int N, int M, const std::function &weight) { + auto minima = SMAWK(N, M, [&](int i, int j) -> T { return weight(i, M - 1 - j); }); + for (auto &p : minima) p.first = M - 1 - p.first; + return minima; +} + // Concave max-plus convolution -// b must be concave +// **b MUST BE CONCAVE** // Complexity: O(n + m) // Verify: https://www.codechef.com/problems/MAXPREFFLIP template @@ -60,10 +71,7 @@ std::vector concave_max_plus_convolution(const std::vector &a, const std:: } return true; }; - - bool a_concave = is_concave(a), b_concave = is_concave(b); - assert(a_concave or b_concave); - if (!b_concave) return concave_max_plus_convolution(b, a); + assert(is_concave(b)); auto select = [&](int i, int j) -> S { int aidx = j, bidx = i - j; @@ -75,3 +83,61 @@ std::vector concave_max_plus_convolution(const std::vector &a, const std:: for (auto x : minima) ret.push_back(-x.second); return ret; } + +// Concave min-plus convolution +// **b MUST BE CONCAVE** +// Complexity: O((n + m)log(n + m)) +template +std::vector concave_min_plus_convolution(const std::vector &a, const std::vector &b) { + const int n = a.size(), m = b.size(); + + auto is_concave = [&](const std::vector &u) -> bool { + for (int i = 1; i + 1 < int(u.size()); ++i) { + if (u[i - 1] + u[i + 1] > u[i] + u[i]) return false; + } + return true; + }; + assert(is_concave(b)); + + std::vector ret(n + m - 1); + std::vector argmin(n + m - 1, -1); + + // mat[i][j] = a[j] + b[i - j] + auto is_valid = [&](int i, int j) { return 0 <= i - j and i - j < m; }; + auto has_valid = [&](int il, int ir, int jl, int jr) { + if (il >= ir or jl >= jr) return false; + return is_valid(il, jl) or is_valid(il, jr - 1) or is_valid(ir - 1, jl) or + is_valid(ir - 1, jr - 1); + }; + + auto rec = [&](auto &&self, int il, int ir, int jl, int jr) -> void { + if (!has_valid(il, ir, jl, jr)) return; + + if (is_valid(il, jr - 1) and is_valid(ir - 1, jl)) { + auto select = [&](int i, int j) -> S { return a[j + jl] + b[(i + il) - (j + jl)]; }; + const auto res = SMAWKAntiMonotone(ir - il, jr - jl, select); + for (int idx = 0; idx < ir - il; ++idx) { + const int i = il + idx; + if (argmin[i] == -1 or res[idx].second < ret[i]) { + ret[i] = res[idx].second; + argmin[i] = res[idx].first + jl; + } + } + } else { + if (const int di = ir - il, dj = jr - jl; di > dj) { + const int im = (il + ir) / 2; + self(self, il, im, jl, jr); + self(self, im, ir, jl, jr); + } else { + const int jm = (jl + jr) / 2; + self(self, il, ir, jl, jm); + self(self, il, ir, jm, jr); + } + } + }; + + rec(rec, 0, n + m - 1, 0, n); + + return ret; + // return argmin; // If you want argmin (0 <= argmin[idx] < len(a)) +} diff --git a/other_algorithms/smawk.md b/other_algorithms/smawk.md index 4a8a3f27..03abb847 100644 --- a/other_algorithms/smawk.md +++ b/other_algorithms/smawk.md @@ -1,9 +1,10 @@ --- -title: Totally Monotone Matrix Searching (SMAWK) +title: Totally Monotone Matrix Searching (SMAWK), concave max-plus / min-plus convolution documentation_of: ./smawk.hpp --- -Totally monotone な $N \times M$ 行列について，各行の最小値の位置を $O(N + M)$ で取得する． +`SMAWK()` 関数は，totally monotone な $N \times M$ 行列について，各行の最小値の位置およびその値を $O(N + M)$ で取得する． +また， totally anti-monotone （不等式の向きが逆，つまり行最小値の位置が行に対して単調非増加などの性質を持つ）な行列に対して同様に最小値を取得する `SMAWKAntiMonotone()` 関数も実装されている． ## 使用方法 @@ -33,10 +34,18 @@ vector a, b; vector c = concave_max_plus_convolution(a, b); ``` +### 応用例：concave min-plus convolution + +上記と同様の状況設定で，逆に min-plus convolution も SMAWK を応用することで $O((n + m) \log (n + m))$ で計算できる． + +このとき，SMAWK を適用したい仮想的な $(n + m - 1) \times n$ 行列は，無効値の位置の都合が悪く totally monotone でも totally anti-monotone でもないため直接扱えない．ここで，有効値が入った平行四辺形の領域をうまく矩形に分割していくことで `SMAWKAntiMonotone()` 関数が適用可能な状況にしている（この分割統治で計算量に log がつく）． + ## 問題例 - [Communication Channel \| CodeChef](https://www.codechef.com/problems/COMMCHA) - [Maximal Prefix After Flip \| CodeChef](https://www.codechef.com/problems/MAXPREFFLIP) +- [Min Plus Convolution (Convex and Arbitrary)](https://judge.yosupo.jp/problem/min_plus_convolution_convex_arbitrary) +- [Min Plus Convolution (Concave and Arbitrary)](https://judge.yosupo.jp/problem/min_plus_convolution_concave_arbitrary) ## Links diff --git a/other_algorithms/test/concave_min_plus_convolution.test.cpp b/other_algorithms/test/concave_min_plus_convolution.test.cpp new file mode 100644 index 00000000..404a33c7 --- /dev/null +++ b/other_algorithms/test/concave_min_plus_convolution.test.cpp @@ -0,0 +1,21 @@ +#define PROBLEM "https://judge.yosupo.jp/problem/min_plus_convolution_concave_arbitrary" + +#include "../smawk.hpp" + +#include +#include +using namespace std; + +int main() { + cin.tie(nullptr), ios::sync_with_stdio(false); + + int N, M; + cin >> N >> M; + vector A(N), B(M); + for (auto &a : A) cin >> a; + for (auto &b : B) cin >> b; + + vector ret = concave_min_plus_convolution(B, A); + + for (int i = 0; i < N + M - 1; ++i) cout << ret[i] << " \n"[i + 1 == N + M - 1]; +}