feature/smark min plus convolution

docs/min-plus-convolution-concave-arbitary.md

@@ -0,0 +1,41 @@
+---
+title: Min Plus Convolution (Concave and Arbitary)
+documentation_of: //dp/min-plus-convolution-concave-arbitary.hpp
+---
+
+凹数列と任意の数列の min-plus 畳み込みを、分割統治と SMAWK により準線形時間で計算する。有効な添字の領域を長方形に分割し、各長方形が表す全単調行列の行最小値を SMAWK で求める。
+
+# min_plus_convolution_concave_arbitary
+
+```cpp
+template <typename T>
+vector<T> min_plus_convolution_concave_arbitary(const vector<T>& a,
+                                                const vector<T>& b)
+```
+
+$a$ を凹数列、$b$ を任意の数列として、各 $k$ に対する $\min_{i+j=k}(a_i+b_j)$ を並べた配列を返す。どちらかが空なら空配列を返す。
+
+## テンプレート引数
+
+- `T`: 加算と `<` による比較が可能で、`std::numeric_limits<T>::max()` が利用可能な要素型
+
+## 引数
+
+- `a`: 隣接差分が広義単調減少する凹数列
+- `b`: 任意の数列
+
+## 戻り値
+
+両方の入力が空でない場合、長さ $\lvert a\rvert + \lvert b\rvert - 1$ の min-plus 畳み込みを返す。
+
+## 前提条件
+
+- `a` は凹数列である
+- 配列長と要素の加算結果は、それぞれ `int` と `T` で表現できる
+
+## 計算量
+
+$N = \lvert a\rvert$, $M = \lvert b\rvert$ とする。
+
+- 時間: $O((N + M)\log(N + M))$
+- 空間: $O(N + M)$

docs/min-plus-convolution-convex-arbitary.md

@@ -0,0 +1,41 @@
+---
+title: Min Plus Convolution (Convex and Arbitary)
+documentation_of: //dp/min-plus-convolution-convex-arbitary.hpp
+---
+
+凸数列と任意の数列の min-plus 畳み込みを SMAWK により線形時間で計算する。
+
+# min_plus_convolution_convex_arbitary
+
+```cpp
+template <typename T>
+vector<T> min_plus_convolution_convex_arbitary(const vector<T>& a,
+                                               const vector<T>& b)
+```
+
+$a$ を凸数列、$b$ を任意の数列として、各 $k$ に対する $\min_{i+j=k}(a_i+b_j)$ を並べた配列を返す。どちらかが空なら空配列を返す。
+
+## テンプレート引数
+
+- `T`: 加算と `<` による比較が可能な要素型
+
+## 引数
+
+- `a`: 隣接差分が広義単調増加する凸数列
+- `b`: 任意の数列
+
+## 戻り値
+
+両方の入力が空でない場合、長さ $\lvert a\rvert + \lvert b\rvert - 1$ の min-plus 畳み込みを返す。
+
+## 前提条件
+
+- `a` は凸数列である
+- 配列長と要素の加算結果は、それぞれ `int` と `T` で表現できる
+
+## 計算量
+
+$N = \lvert a\rvert$, $M = \lvert b\rvert$ とする。
+
+- 時間: $O(N + M)$
+- 空間: $O(N + M)$

docs/monotone-minima.md

@@ -3,18 +3,68 @@ title: Monotone Minima
 documentation_of: //dp/monotone-minima.hpp
 ---
 
-$2$ 変数関数 $f(i, j) (0 \leq i \lt H, 0 \leq j \lt W)$ が Monotone であるとは、すべての $k$ に対して $\mathrm{argmin} f(k, *) \leq \mathrm{argmin} f(k + 1, *)$ を満たすことをいう。つまり各行の最小値をとる位置が右下に単調に下がっていることを意味する。
+$H \times W$ 行列の各行について、最適な列を分割統治で求める。最適列の位置が行番号に対して広義単調増加する行列に利用できる。
 
 Monge $\Rightarrow$ Totally Monotone(TM) $\Rightarrow$ Monotone なので、Monotone は弱い条件である。
 
 # monotone_minima
 
 ```cpp
-vector<pair<int, T> > monotone_minima(int H, int W, const function<T(int, int)>& f, const Compare& comp = Compare())
+template <typename F>
+vector<int> monotone_minima(int H, int W, F comp)
 ```
 
-各行について、最小値をとる位置と最小値をペアで返す。`f` は $2$ 変数関数、`comp` は比較関数。
+各行の最適な列番号を返す。`comp(i, j, k)` は、行 `i` において列 `k` が列 `j` より真に良いとき `true` を返すものとする。このインターフェースは `smawk` と共通である。同値な候補では左側の列を選ぶ。
+
+## 引数
+
+- `H`: 行数
+- `W`: 列数
+- `comp`: 2 列の優劣を判定する関数
+
+## 戻り値
+
+長さ $H$ の配列を返し、その第 $i$ 要素は行 $i$ の最適な列番号である。$W = 0$ の場合は、すべての要素が $-1$ となる。
+
+## 制約
+
+- $0 \leq H$
+- $0 \leq W$
+- 各行の最適列が広義単調増加する
+
+## 計算量
+
+- 時間: $O(W \log H + H)$ 回の `comp` 呼び出し
+- 空間: $O(H)$
+
+# monotone_minima_select
+
+```cpp
+template <typename Select>
+vector<int> monotone_minima_select(int H, int W, Select select)
+```
+
+各行の最適な列番号を返す。`select(i, l, r)` は、行 $i$ の半開区間 $[l, r)$ に含まれる最適な列番号を返すものとする。各行について候補列が一つの連続区間として渡されるため、列を進めながら評価値を更新できる場合に利用できる。
+
+## 引数
+
+- `H`: 行数
+- `W`: 列数
+- `select`: 指定された行と列区間から最適列を求める関数
+
+## 戻り値
+
+長さ $H$ の配列を返し、その第 $i$ 要素は行 $i$ の最適な列番号である。$W = 0$ の場合は、`select` を呼ばず、すべての要素が $-1$ の配列を返す。
+
+## 制約
+
+- $0 \leq H$
+- $0 \leq W$
+- 各行の最適列が広義単調増加する
+- `select(i, l, r)` は $l \leq j < r$ を満たす最適列 $j$ を返す
 
 ## 計算量
 
-- $O(N \log N)$
+- `select` を $H$ 回呼び出す
+- 渡される区間長の総和は $O(W \log H + H)$
+- 空間: $O(H)$

docs/smawk.md

@@ -0,0 +1,40 @@
+---
+title: SMAWK
+documentation_of: //dp/smawk.hpp
+---
+
+全単調行列の各行について、最適な列を線形時間で求める。行列の要素そのものを保持せず、2 列の優劣を判定する関数だけを受け取る。
+
+# smawk
+
+```cpp
+template <typename F>
+vector<int> smawk(int H, int W, F comp)
+```
+
+各行の最適な列番号を返す。`comp(i, j, k)` は、行 `i` において列 `k` が列 `j` より真に良いとき `true` を返すものとする。同値な候補では左側の列を選ぶ。
+
+## 引数
+
+- `H`: 行数
+- `W`: 列数
+- `comp`: 2 列の優劣を判定する関数
+
+## 戻り値
+
+長さ $H$ の配列を返し、その第 $i$ 要素は行 $i$ の最適な列番号である。$W = 0$ の場合は、すべての要素が $-1$ となる。
+
+## 制約
+
+- $0 \leq H$
+- $0 \leq W$
+- 行列が `comp` の定める順序について全単調である
+
+## 計算量
+
+- 時間: $O(H + W)$ 回の `comp` 呼び出し
+- 空間: $O(H + W)$
+
+# 参考文献
+
+- Aggarwal, Klawe, Moran, Shor, Wilber, Geometric Applications of a Matrix-Searching Algorithm

dp/divide-and-conquer-optimization.hpp

@@ -16,8 +16,10 @@ std::vector<std::vector<T> > divide_and_conquer_optimization(
       if (x >= y) return INF;
       return dp[i - 1][x] + f(x, y);
     };
-    auto ret = monotone_minima(W + 1, W + 1, get_cost, comp);
-    for (int j = 0; j <= W; j++) dp[i][j] = ret[j].second;
+    auto ret = monotone_minima(W + 1, W + 1, [&](int j, int old_k, int new_k) {
+      return comp(get_cost(j, new_k), get_cost(j, old_k));
+    });
+    for (int j = 0; j <= W; j++) dp[i][j] = get_cost(j, ret[j]);
   }
   return dp;
 }

dp/min-plus-convolution-concave-arbitary.hpp

@@ -0,0 +1,72 @@
+#pragma once
+
+#include <algorithm>
+#include <limits>
+#include <vector>
+
+#include "smawk.hpp"
+
+template <typename T>
+std::vector<T> min_plus_convolution_concave_arbitary(const std::vector<T>& a,
+                                                     const std::vector<T>& b) {
+  if (a.empty() || b.empty()) return {};
+  int N = static_cast<int>(a.size());
+  int M = static_cast<int>(b.size());
+  int H = N + M - 1;
+
+  std::vector<int> column_min(H, 0), column_max(H, M - 1);
+  for (int row = N; row < H; ++row) column_min[row] = row - N + 1;
+  for (int row = 0; row <= H - N; ++row) column_max[row] = row;
+
+  std::vector<int> row_min(M), row_max(M);
+  for (int column = 0; column < M; ++column) {
+    row_min[column] = column;
+    row_max[column] = N - 1 + column;
+  }
+
+  std::vector<T> result(H, std::numeric_limits<T>::max());
+  auto divide = [&](auto&& self, int row_left, int row_right, int column_left,
+                    int column_right) -> void {
+    if (column_max[row_left] >= column_right &&
+        column_left >= column_min[row_right]) {
+      auto value = [&](int row, int column) {
+        int j = column_right - column;
+        return b[j] + a[row_left + row - j];
+      };
+      auto argmin =
+          smawk(row_right - row_left + 1, column_right - column_left + 1,
+                [&](int row, int old_column, int new_column) {
+                  return value(row, new_column) < value(row, old_column);
+                });
+      for (int row = row_left; row <= row_right; ++row) {
+        result[row] = std::min(result[row],
+                               value(row - row_left, argmin[row - row_left]));
+      }
+      return;
+    }
+
+    if (row_right - row_left > column_right - column_left) {
+      int row_middle = (row_left + row_right) / 2;
+      int next_column_right = std::min(column_max[row_middle], column_right);
+      if (column_left <= next_column_right) {
+        self(self, row_left, row_middle, column_left, next_column_right);
+      }
+      int next_column_left = std::max(column_min[row_middle], column_left);
+      if (next_column_left <= column_right) {
+        self(self, row_middle + 1, row_right, next_column_left, column_right);
+      }
+    } else {
+      int column_middle = (column_left + column_right) / 2;
+      int next_row_right = std::min(row_max[column_middle], row_right);
+      if (row_left <= next_row_right) {
+        self(self, row_left, next_row_right, column_left, column_middle);
+      }
+      int next_row_left = std::max(row_min[column_middle], row_left);
+      if (next_row_left <= row_right) {
+        self(self, next_row_left, row_right, column_middle + 1, column_right);
+      }
+    }
+  };
+  divide(divide, 0, H - 1, 0, M - 1);
+  return result;
+}

dp/min-plus-convolution-convex-arbitary.hpp

@@ -0,0 +1,24 @@
+#pragma once
+
+#include <vector>
+
+#include "smawk.hpp"
+
+template <typename T>
+std::vector<T> min_plus_convolution_convex_arbitary(const std::vector<T>& a,
+                                                    const std::vector<T>& b) {
+  if (a.empty() || b.empty()) return {};
+  int H = static_cast<int>(a.size());
+  int W = static_cast<int>(b.size());
+  const auto c = smawk(H + W - 1, W, [&](int i, int j, int k) {
+    if (i < k) return false;
+    if (i - j >= H) return true;
+    return b[k] + a[i - k] < b[j] + a[i - j];
+  });
+  std::vector<T> ret;
+  ret.reserve(H + W - 1);
+  for (int i = 0; i < H + W - 1; ++i) {
+    ret.emplace_back(b[c[i]] + a[i - c[i]]);
+  }
+  return ret;
+}

dp/monotone-minima.hpp

@@ -1,31 +1,31 @@
 #pragma once
 
-#include <functional>
-#include <utility>
 #include <vector>
 
-template <typename T, typename Compare = std::less<T> >
-std::vector<std::pair<int, T> > monotone_minima(
-    int H, int W, const std::function<T(int, int)>& f,
-    const Compare& comp = Compare()) {
-  std::vector<std::pair<int, T> > dp(H);
-  std::function<void(int, int, int, int)> dfs = [&](int top, int bottom,
-                                                    int left, int right) {
+template <typename Select>
+std::vector<int> monotone_minima_select(int H, int W, Select select) {
+  std::vector<int> ret(H, -1);
+  if (H == 0 || W == 0) return ret;
+  auto dfs = [&](auto&& self, int top, int bottom, int left,
+                 int right) -> void {
     if (top > bottom) return;
     int line = (top + bottom) / 2;
-    T ma;
-    int mi = -1;
-    for (int i = left; i <= right; i++) {
-      T cst = f(line, i);
-      if (mi == -1 || comp(cst, ma)) {
-        ma = cst;
-        mi = i;
-      }
-    }
-    dp[line] = std::make_pair(mi, ma);
-    dfs(top, line - 1, left, mi);
-    dfs(line + 1, bottom, mi, right);
+    int best = select(line, left, right + 1);
+    ret[line] = best;
+    self(self, top, line - 1, left, best);
+    self(self, line + 1, bottom, best, right);
   };
-  dfs(0, H - 1, 0, W - 1);
-  return dp;
+  dfs(dfs, 0, H - 1, 0, W - 1);
+  return ret;
+}
+
+template <typename F>
+std::vector<int> monotone_minima(int H, int W, F comp) {
+  return monotone_minima_select(H, W, [&](int row, int left, int right) {
+    int best = left;
+    for (int column = left + 1; column < right; ++column) {
+      if (comp(row, best, column)) best = column;
+    }
+    return best;
+  });
 }

dp/online-offline-dp.hpp

@@ -19,11 +19,15 @@ std::vector<T> online_offline_dp(int W, const std::function<T(int, int)>& f,
       [&](int l, int m,
           int r) {  // dp[l, m) -> dp[m, r)
         x_base = l, y_base = m;
-        auto ret = monotone_minima(r - m, m - l, get_cost, comp);
+        auto ret =
+            monotone_minima(r - m, m - l, [&](int i, int old_j, int new_j) {
+              return comp(get_cost(i, new_j), get_cost(i, old_j));
+            });
         for (int i = 0; i < ret.size(); i++) {
-          if (!isset[m + i] || comp(ret[i].second, dp[m + i])) {
+          T cost = get_cost(i, ret[i]);
+          if (!isset[m + i] || comp(cost, dp[m + i])) {
             isset[m + i] = true;
-            dp[m + i] = ret[i].second;
+            dp[m + i] = cost;
           }
         }
       };