enhance SMAWK / add monotone minima

Author: hitonanode

other_algorithms/monotone_minima.hpp

@@ -0,0 +1,32 @@
+#pragma once
+#include <cassert>
+#include <functional>
+#include <vector>
+
+// 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 <class T>
+std::vector<std::pair<int, T>>
+MonotoneMinima(int N, int M, const std::function<T(int i, int j)> &weight) {
+    std::vector<std::pair<int, T>> 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;
+}

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<T> dp(N, inf);
+dp[0] = 0;
+
+auto weight = [&](int s, int t) -> T { /* Implement */ };
+
+const auto res = MonotoneMinima<T>(
+    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)

other_algorithms/smawk.hpp

@@ -6,7 +6,7 @@
 #include <vector>
 
 // 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<std::pair<int, T>> SMAWK(int N, int M, const std::function<T(int i,
     return minima;
 }
 
+// Find minima of totally ANTI-monotone function f(i, j) (0 <= i < N, 0 <= j < M) for each i
+// Constraints: every submatrix of f(i, j) is anti-monotone.
+// Complexity: O(N + M)
+template <class T>
+std::vector<std::pair<int, T>>
+SMAWKAntiMonotone(int N, int M, const std::function<T(int i, int j)> &weight) {
+    auto minima = SMAWK<T>(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 <class S, S INF>
@@ -60,10 +71,7 @@ std::vector<S> concave_max_plus_convolution(const std::vector<S> &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<S, INF>(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<S> concave_max_plus_convolution(const std::vector<S> &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 <class S>
+std::vector<S> concave_min_plus_convolution(const std::vector<S> &a, const std::vector<S> &b) {
+    const int n = a.size(), m = b.size();
+
+    auto is_concave = [&](const std::vector<S> &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<S> ret(n + m - 1);
+    std::vector<int> 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<S>(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))
+}

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<int> a, b;
 vector<int> c = concave_max_plus_convolution<int, 1 << 30>(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
 

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 <iostream>
+#include <vector>
+using namespace std;
+
+int main() {
+    cin.tie(nullptr), ios::sync_with_stdio(false);
+
+    int N, M;
+    cin >> N >> M;
+    vector<int> A(N), B(M);
+    for (auto &a : A) cin >> a;
+    for (auto &b : B) cin >> b;
+
+    vector<int> ret = concave_min_plus_convolution<int>(B, A);
+
+    for (int i = 0; i < N + M - 1; ++i) cout << ret[i] << " \n"[i + 1 == N + M - 1];
+}