diff --git a/docs/min-plus-convolution-concave-arbitary.md b/docs/min-plus-convolution-concave-arbitary.md new file mode 100644 index 00000000..837f7876 --- /dev/null +++ b/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 +vector min_plus_convolution_concave_arbitary(const vector& a, + const vector& b) +``` + +$a$ を凹数列、$b$ を任意の数列として、各 $k$ に対する $\min_{i+j=k}(a_i+b_j)$ を並べた配列を返す。どちらかが空なら空配列を返す。 + +## テンプレート引数 + +- `T`: 加算と `<` による比較が可能で、`std::numeric_limits::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)$ diff --git a/docs/min-plus-convolution-convex-arbitary.md b/docs/min-plus-convolution-convex-arbitary.md new file mode 100644 index 00000000..8dd092db --- /dev/null +++ b/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 +vector min_plus_convolution_convex_arbitary(const vector& a, + const vector& 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)$ diff --git a/docs/monotone-minima.md b/docs/monotone-minima.md index 9b73386e..79a762ee 100644 --- a/docs/monotone-minima.md +++ b/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 > monotone_minima(int H, int W, const function& f, const Compare& comp = Compare()) +template +vector 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 +vector 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)$ diff --git a/docs/smawk.md b/docs/smawk.md new file mode 100644 index 00000000..9fde6b5d --- /dev/null +++ b/docs/smawk.md @@ -0,0 +1,40 @@ +--- +title: SMAWK +documentation_of: //dp/smawk.hpp +--- + +全単調行列の各行について、最適な列を線形時間で求める。行列の要素そのものを保持せず、2 列の優劣を判定する関数だけを受け取る。 + +# smawk + +```cpp +template +vector 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 diff --git a/dp/divide-and-conquer-optimization.hpp b/dp/divide-and-conquer-optimization.hpp index b83f6d9d..a3ab6825 100644 --- a/dp/divide-and-conquer-optimization.hpp +++ b/dp/divide-and-conquer-optimization.hpp @@ -16,8 +16,10 @@ std::vector > 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; } diff --git a/dp/min-plus-convolution-concave-arbitary.hpp b/dp/min-plus-convolution-concave-arbitary.hpp new file mode 100644 index 00000000..70665f08 --- /dev/null +++ b/dp/min-plus-convolution-concave-arbitary.hpp @@ -0,0 +1,72 @@ +#pragma once + +#include +#include +#include + +#include "smawk.hpp" + +template +std::vector min_plus_convolution_concave_arbitary(const std::vector& a, + const std::vector& b) { + if (a.empty() || b.empty()) return {}; + int N = static_cast(a.size()); + int M = static_cast(b.size()); + int H = N + M - 1; + + std::vector 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 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 result(H, std::numeric_limits::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; +} diff --git a/dp/min-plus-convolution-convex-arbitary.hpp b/dp/min-plus-convolution-convex-arbitary.hpp new file mode 100644 index 00000000..03dc70f8 --- /dev/null +++ b/dp/min-plus-convolution-convex-arbitary.hpp @@ -0,0 +1,24 @@ +#pragma once + +#include + +#include "smawk.hpp" + +template +std::vector min_plus_convolution_convex_arbitary(const std::vector& a, + const std::vector& b) { + if (a.empty() || b.empty()) return {}; + int H = static_cast(a.size()); + int W = static_cast(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 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; +} diff --git a/dp/monotone-minima.hpp b/dp/monotone-minima.hpp index 71a73ba5..fbf4ecc8 100644 --- a/dp/monotone-minima.hpp +++ b/dp/monotone-minima.hpp @@ -1,31 +1,31 @@ #pragma once -#include -#include #include -template > -std::vector > monotone_minima( - int H, int W, const std::function& f, - const Compare& comp = Compare()) { - std::vector > dp(H); - std::function dfs = [&](int top, int bottom, - int left, int right) { +template +std::vector monotone_minima_select(int H, int W, Select select) { + std::vector 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 +std::vector 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; + }); } diff --git a/dp/online-offline-dp.hpp b/dp/online-offline-dp.hpp index d2cdf6d4..dda17577 100644 --- a/dp/online-offline-dp.hpp +++ b/dp/online-offline-dp.hpp @@ -19,11 +19,15 @@ std::vector online_offline_dp(int W, const std::function& 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; } } }; diff --git a/dp/smawk.hpp b/dp/smawk.hpp new file mode 100644 index 00000000..54cf3049 --- /dev/null +++ b/dp/smawk.hpp @@ -0,0 +1,58 @@ +#pragma once + +#include +#include +#include + +template +std::vector smawk(int H, int W, F comp) { + std::vector ret(H, -1); + if (H == 0 || W == 0) return ret; + + auto dfs = [&](auto&& self, const std::vector& rows, + const std::vector& cols) -> void { + if (rows.empty()) return; + std::vector reduced; + reduced.reserve(std::min(rows.size(), cols.size())); + for (int c : cols) { + while (!reduced.empty()) { + int r = rows[reduced.size() - 1]; + int old_c = reduced.back(); + if (comp(r, old_c, c)) { + reduced.pop_back(); + } else { + break; + } + } + if (reduced.size() < rows.size()) reduced.emplace_back(c); + } + + std::vector odd_rows; + odd_rows.reserve(rows.size() / 2); + for (int i = 1; i < static_cast(rows.size()); i += 2) { + odd_rows.emplace_back(rows[i]); + } + self(self, odd_rows, reduced); + + int left = 0; + for (int i = 0; i < static_cast(rows.size()); i += 2) { + int right = static_cast(reduced.size()) - 1; + if (i + 1 < static_cast(rows.size())) { + right = left; + while (reduced[right] != ret[rows[i + 1]]) ++right; + } + int best = left; + for (int p = left + 1; p <= right; ++p) { + if (comp(rows[i], reduced[best], reduced[p])) best = p; + } + ret[rows[i]] = reduced[best]; + left = right; + } + }; + + std::vector rows(H), cols(W); + std::iota(rows.begin(), rows.end(), 0); + std::iota(cols.begin(), cols.end(), 0); + dfs(dfs, rows, cols); + return ret; +} diff --git a/test/unittest/monotone-minima.test.cpp b/test/unittest/monotone-minima.test.cpp new file mode 100644 index 00000000..f136a405 --- /dev/null +++ b/test/unittest/monotone-minima.test.cpp @@ -0,0 +1,61 @@ +// competitive-verifier: STANDALONE + +#include "../../dp/monotone-minima.hpp" + +#include +#include +#include + +int main() { + const std::vector> matrix = { + {0, 1, 4, 9, 16}, + {4, 1, 0, 1, 4}, + {16, 9, 4, 1, 0}, + }; + auto argmin = + monotone_minima(3, 5, [&](int row, int old_column, int new_column) { + return matrix[row][new_column] < matrix[row][old_column]; + }); + assert((argmin == std::vector{0, 2, 4})); + auto selected = + monotone_minima_select(3, 5, [&](int row, int left, int right) { + int best = left; + for (int column = left + 1; column < right; ++column) { + if (matrix[row][column] < matrix[row][best]) best = column; + } + return best; + }); + assert(selected == argmin); + + auto empty_rows = monotone_minima(0, 5, [](int, int, int) { return false; }); + assert(empty_rows.empty()); + auto empty_columns = + monotone_minima(3, 0, [](int, int, int) { return false; }); + assert((empty_columns == std::vector{-1, -1, -1})); + int select_calls = 0; + auto empty_select = monotone_minima_select( + 3, 0, [&](int, int, int) { return ++select_calls; }); + assert((empty_select == std::vector{-1, -1, -1})); + assert(select_calls == 0); + + std::mt19937 random(123456789); + for (int height = 1; height <= 30; ++height) { + for (int width = 1; width <= 30; ++width) { + std::vector center(height); + for (int row = 1; row < height; ++row) { + center[row] = center[row - 1] + random() % 3; + } + for (int& column : center) column %= width; + for (int row = 1; row < height; ++row) { + if (center[row] < center[row - 1]) center[row] = center[row - 1]; + } + auto result = monotone_minima( + height, width, [&](int row, int old_column, int new_column) { + int old_distance = old_column - center[row]; + int new_distance = new_column - center[row]; + return new_distance * new_distance < old_distance * old_distance; + }); + assert(result == center); + } + } +} diff --git a/test/verify/yosupo-min-plus-convolution-concave-arbitrary.test.cpp b/test/verify/yosupo-min-plus-convolution-concave-arbitrary.test.cpp new file mode 100644 index 00000000..752031fc --- /dev/null +++ b/test/verify/yosupo-min-plus-convolution-concave-arbitrary.test.cpp @@ -0,0 +1,22 @@ +// clang-format off +// competitive-verifier: PROBLEM https://judge.yosupo.jp/problem/min_plus_convolution_concave_arbitrary +// clang-format on + +#include +#include + +#include "../../dp/min-plus-convolution-concave-arbitary.hpp" + +int main() { + int N, M; + std::cin >> N >> M; + std::vector A(N), B(M); + for (int& a : A) std::cin >> a; + for (int& b : B) std::cin >> b; + auto C = min_plus_convolution_concave_arbitary(A, B); + for (int i = 0; i < static_cast(C.size()); ++i) { + if (i) std::cout << ' '; + std::cout << C[i]; + } + std::cout << '\n'; +} diff --git a/test/verify/yosupo-min-plus-convolution-convex-arbitrary.test.cpp b/test/verify/yosupo-min-plus-convolution-convex-arbitrary.test.cpp new file mode 100644 index 00000000..a834656c --- /dev/null +++ b/test/verify/yosupo-min-plus-convolution-convex-arbitrary.test.cpp @@ -0,0 +1,22 @@ +// clang-format off +// competitive-verifier: PROBLEM https://judge.yosupo.jp/problem/min_plus_convolution_convex_arbitrary +// clang-format on + +#include +#include + +#include "../../dp/min-plus-convolution-convex-arbitary.hpp" + +int main() { + int N, M; + std::cin >> N >> M; + std::vector A(N), B(M); + for (long long& a : A) std::cin >> a; + for (long long& b : B) std::cin >> b; + auto C = min_plus_convolution_convex_arbitary(A, B); + for (int i = 0; i < static_cast(C.size()); ++i) { + if (i) std::cout << ' '; + std::cout << C[i]; + } + std::cout << '\n'; +}