max_plus_convolution.hpp¶
Max-plus convolution of two concave sequences.
Verified by min_plus_convolution_concave_arbitrary, min_plus_convolution_convex_arbitrary, min_plus_convolution_convex_convex.
计算两个凹序列的 max-plus 卷积；用于具有四边形不等式/凹性的 DP 合并。
Implementation¶
#ifndef NOYA_MAXPLUS_CONVOLUTION_HPP
#define NOYA_MAXPLUS_CONVOLUTION_HPP 1
/// @complexity Time: O(n + m) for the supported convex/concave cases.
/// Space: O(n + m) output and monotone-optimum workspace.
#include "noya/smawk.hpp"
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <functional>
#include <limits>
#include <vector>
namespace noya {
/// @brief Max-plus convolution of two concave sequences.
template <class T>
std::vector<T> two_concave_maxplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
if (a.empty())
return b;
if (b.empty())
return a;
const int n = int(a.size());
const int m = int(b.size());
int p = 0, q = 0;
std::vector<T> c(n + m - 1);
c[0] = a[0] + b[0];
for (int i = 1; i < n + m - 1; i++) {
if (p + 1 == n) {
q++;
} else if (q + 1 == m) {
p++;
} else {
if (a[p + 1] - a[p] > b[q + 1] - b[q]) {
p++;
} else {
q++;
}
}
c[i] = a[p] + b[q];
}
return c;
}
/// @brief Min-plus convolution of two convex sequences.
/// The first differences of a convex sequence are nondecreasing. Merging the
/// two difference sequences therefore describes, in order, every step of the
/// lower boundary of their min-plus convolution.
template <class T>
std::vector<T> convex_convex_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
if (a.empty())
return b;
if (b.empty())
return a;
const int n = int(a.size());
const int m = int(b.size());
std::vector<T> negative_a(n);
std::vector<T> negative_b(m);
for (int i = 0; i < n; i++) {
negative_a[i] = -a[i];
}
for (int i = 0; i < m; i++) {
negative_b[i] = -b[i];
}
auto negative_c = two_concave_maxplus_convolution(negative_a, negative_b);
std::vector<T> c(n + m - 1);
for (int i = 0; i < n + m - 1; i++)
c[i] = -negative_c[i];
return c;
}
/// @brief Backward-compatible name for convex-convex min-plus convolution.
template <class T>
std::vector<T> two_concave_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
return convex_convex_minplus_convolution(a, b);
}
/// @brief Max-plus convolution where b is concave.
/// Concavity makes the implicit matrix a[i] + b[row-i] totally monotone, so
/// SMAWK finds all row maxima after only O(n + m) value comparisons.
template <class T>
std::vector<T> concave_maxplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
if (a.empty())
return b;
if (b.empty())
return a;
const int n = int(a.size());
const int m = int(b.size());
const auto get = [&](const int &i, const int &j) -> T {
return a[j] + b[i - j];
};
const auto select = [&](const int &i, const int &j, const int &k) -> bool {
if (i < k)
return false;
if (i - j >= m)
return true;
return get(i, j) <= get(i, k);
};
const auto amax = smawk(n + m - 1, n, select);
std::vector<T> c(n + m - 1);
for (int i = 0; i < n + m - 1; i++)
c[i] = get(i, amax[i]);
return c;
}
/// @brief Min-plus convolution of an arbitrary sequence and a convex sequence.
/// Convexity makes the implicit matrix a[i] + b[row-i] totally monotone.
/// Negating both inputs turns row minima into row maxima, which SMAWK finds in
/// linear time without materializing the matrix.
template <class T>
std::vector<T> arbitrary_convex_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
if (a.empty())
return b;
if (b.empty())
return a;
const int n = int(a.size());
const int m = int(b.size());
std::vector<T> negative_a(n);
std::vector<T> negative_b(m);
for (int i = 0; i < n; i++) {
negative_a[i] = -a[i];
}
for (int i = 0; i < m; i++) {
negative_b[i] = -b[i];
}
auto negative_c = concave_maxplus_convolution(negative_a, negative_b);
std::vector<T> c(n + m - 1);
for (int i = 0; i < n + m - 1; i++)
c[i] = -negative_c[i];
return c;
}
/// @brief Min-plus convolution of a convex sequence and an arbitrary sequence.
template <class T>
std::vector<T> convex_arbitrary_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
return arbitrary_convex_minplus_convolution(b, a);
}
namespace max_plus_convolution_internal {
template <class Value>
std::vector<int> monotone_row_minima(int rows, int columns, Value value) {
std::vector<int> result(rows);
int stride = 1;
while (stride < rows) {
stride <<= 1;
}
for (; stride > 0; stride >>= 1) {
for (int row = stride - 1; row < rows; row += 2 * stride) {
int first = row >= stride ? result[row - stride] : 0;
int last = row + stride < rows ? result[row + stride] : columns - 1;
result[row] = first;
for (int column = first + 1; column <= last; column++) {
if (value(row, column) < value(row, result[row])) {
result[row] = column;
}
}
}
}
return result;
}
template <class T>
std::vector<T>
arbitrary_concave_minplus_convolution(const std::vector<T> &arbitrary,
const std::vector<T> &concave) {
if (arbitrary.empty()) {
return concave;
}
if (concave.empty()) {
return arbitrary;
}
for (int i = 0; i + 2 < int(concave.size()); i++) {
assert(concave[i + 1] - concave[i] >=
concave[i + 2] - concave[i + 1]);
}
const int width = int(arbitrary.size());
const int concave_size = int(concave.size());
const int height = width + concave_size - 1;
std::vector<int> minimum_column(height), maximum_column(height, width - 1);
for (int row = concave_size; row < height; row++) {
minimum_column[row] = row - concave_size + 1;
}
for (int row = 0; row <= height - concave_size; row++) {
maximum_column[row] = row;
}
std::vector<int> minimum_row(width), maximum_row(width);
for (int column = 0; column < width; column++) {
minimum_row[column] = column;
maximum_row[column] = concave_size - 1 + column;
}
std::vector<T> result(height, std::numeric_limits<T>::max());
std::function<void(int, int, int, int)> solve =
[&](int first_row, int last_row, int first_column, int last_column) {
if (maximum_column[first_row] >= last_column &&
first_column >= minimum_column[last_row]) {
auto value = [&](int local_row, int reversed_column) {
int column = last_column - reversed_column;
int row = first_row + local_row;
return arbitrary[column] + concave[row - column];
};
auto minima = monotone_row_minima(last_row - first_row + 1,
last_column - first_column + 1,
value);
for (int row = first_row; row <= last_row; row++) {
result[row] = std::min(result[row],
value(row - first_row,
minima[row - first_row]));
}
return;
}
if (std::int64_t(last_row - first_row) *
(last_column - first_column) <
1024) {
for (int row = first_row; row <= last_row; row++) {
int from = std::max(minimum_column[row], first_column);
int to = std::min(maximum_column[row], last_column);
for (int column = from; column <= to; column++) {
result[row] = std::min(
result[row], arbitrary[column] + concave[row - column]);
}
}
return;
}
if (last_row - first_row > last_column - first_column) {
int middle = (first_row + last_row) / 2;
int new_last = std::min(maximum_column[middle], last_column);
if (first_column <= new_last) {
solve(first_row, middle, first_column, new_last);
}
int new_first = std::max(minimum_column[middle], first_column);
if (new_first <= last_column) {
solve(middle + 1, last_row, new_first, last_column);
}
} else {
int middle = (first_column + last_column) / 2;
int new_last = std::min(maximum_row[middle], last_row);
if (first_row <= new_last) {
solve(first_row, new_last, first_column, middle);
}
int new_first = std::max(minimum_row[middle], first_row);
if (new_first <= last_row) {
solve(new_first, last_row, middle + 1, last_column);
}
}
};
solve(0, height - 1, 0, width - 1);
return result;
}
} // namespace max_plus_convolution_internal
/// @brief Min-plus convolution of a concave sequence and an arbitrary
/// sequence. Valid pairs form a diagonal staircase rather than one rectangular
/// Monge matrix. Recursively splitting that staircase produces fully valid
/// rectangles; after reversing their columns, concavity makes row minima
/// monotone and they are found together. Small boundary rectangles are scanned
/// directly.
template <class T>
std::vector<T> concave_arbitrary_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
return max_plus_convolution_internal::arbitrary_concave_minplus_convolution(
b, a);
}
/// @brief Backward-compatible name for arbitrary-convex min-plus convolution.
template <class T>
std::vector<T> concave_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
return arbitrary_convex_minplus_convolution(a, b);
}
} // namespace noya
#endif // NOYA_MAXPLUS_CONVOLUTION_HPP
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <functional>
#include <limits>
#include <numeric>
#include <vector>
/// @complexity Time: O(n + m) for the supported convex/concave cases.
/// Space: O(n + m) output and monotone-optimum workspace.
/// @complexity Time: O(rows + columns) matrix probes.
/// Space: O(rows + columns).
namespace noya {
/// @brief SMAWK algorithm: compute row minima of a totally monotone matrix.
/// A stack reduction leaves at most one candidate column per row, recursion
/// solves the odd rows, and monotone argmins bound the scan that interpolates
/// each even row. Every row and column is discarded or scanned only O(1) times.
/// @return Vector where ans[i] is the column index of the minimum in row i.
template <class Select>
std::vector<int> smawk(const int row_size, const int col_size,
const Select &select) {
const std::function<std::vector<int>(const std::vector<int> &,
const std::vector<int> &)>
solve = [&](const std::vector<int> &row,
const std::vector<int> &col) -> std::vector<int> {
const int n = int(row.size());
if (n == 0)
return {};
std::vector<int> c2;
for (const int i : col) {
while (!c2.empty() && select(row[c2.size() - 1], c2.back(), i))
c2.pop_back();
if (c2.size() < n)
c2.push_back(i);
}
std::vector<int> r2;
for (int i = 1; i < n; i += 2)
r2.push_back(row[i]);
const std::vector<int> a2 = solve(r2, c2);
std::vector<int> ans(n);
for (int i = 0; i != a2.size(); i += 1)
ans[i * 2 + 1] = a2[i];
int j = 0;
for (int i = 0; i < n; i += 2) {
ans[i] = c2[j];
const int end = i + 1 == n ? c2.back() : ans[i + 1];
while (c2[j] != end) {
j += 1;
if (select(row[i], ans[i], c2[j]))
ans[i] = c2[j];
}
}
return ans;
};
std::vector<int> row(row_size);
std::iota(row.begin(), row.end(), 0);
std::vector<int> col(col_size);
std::iota(col.begin(), col.end(), 0);
return solve(row, col);
}
} // namespace noya
namespace noya {
/// @brief Max-plus convolution of two concave sequences.
template <class T>
std::vector<T> two_concave_maxplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
if (a.empty())
return b;
if (b.empty())
return a;
const int n = int(a.size());
const int m = int(b.size());
int p = 0, q = 0;
std::vector<T> c(n + m - 1);
c[0] = a[0] + b[0];
for (int i = 1; i < n + m - 1; i++) {
if (p + 1 == n) {
q++;
} else if (q + 1 == m) {
p++;
} else {
if (a[p + 1] - a[p] > b[q + 1] - b[q]) {
p++;
} else {
q++;
}
}
c[i] = a[p] + b[q];
}
return c;
}
/// @brief Min-plus convolution of two convex sequences.
/// The first differences of a convex sequence are nondecreasing. Merging the
/// two difference sequences therefore describes, in order, every step of the
/// lower boundary of their min-plus convolution.
template <class T>
std::vector<T> convex_convex_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
if (a.empty())
return b;
if (b.empty())
return a;
const int n = int(a.size());
const int m = int(b.size());
std::vector<T> negative_a(n);
std::vector<T> negative_b(m);
for (int i = 0; i < n; i++) {
negative_a[i] = -a[i];
}
for (int i = 0; i < m; i++) {
negative_b[i] = -b[i];
}
auto negative_c = two_concave_maxplus_convolution(negative_a, negative_b);
std::vector<T> c(n + m - 1);
for (int i = 0; i < n + m - 1; i++)
c[i] = -negative_c[i];
return c;
}
/// @brief Backward-compatible name for convex-convex min-plus convolution.
template <class T>
std::vector<T> two_concave_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
return convex_convex_minplus_convolution(a, b);
}
/// @brief Max-plus convolution where b is concave.
/// Concavity makes the implicit matrix a[i] + b[row-i] totally monotone, so
/// SMAWK finds all row maxima after only O(n + m) value comparisons.
template <class T>
std::vector<T> concave_maxplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
if (a.empty())
return b;
if (b.empty())
return a;
const int n = int(a.size());
const int m = int(b.size());
const auto get = [&](const int &i, const int &j) -> T {
return a[j] + b[i - j];
};
const auto select = [&](const int &i, const int &j, const int &k) -> bool {
if (i < k)
return false;
if (i - j >= m)
return true;
return get(i, j) <= get(i, k);
};
const auto amax = smawk(n + m - 1, n, select);
std::vector<T> c(n + m - 1);
for (int i = 0; i < n + m - 1; i++)
c[i] = get(i, amax[i]);
return c;
}
/// @brief Min-plus convolution of an arbitrary sequence and a convex sequence.
/// Convexity makes the implicit matrix a[i] + b[row-i] totally monotone.
/// Negating both inputs turns row minima into row maxima, which SMAWK finds in
/// linear time without materializing the matrix.
template <class T>
std::vector<T> arbitrary_convex_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
if (a.empty())
return b;
if (b.empty())
return a;
const int n = int(a.size());
const int m = int(b.size());
std::vector<T> negative_a(n);
std::vector<T> negative_b(m);
for (int i = 0; i < n; i++) {
negative_a[i] = -a[i];
}
for (int i = 0; i < m; i++) {
negative_b[i] = -b[i];
}
auto negative_c = concave_maxplus_convolution(negative_a, negative_b);
std::vector<T> c(n + m - 1);
for (int i = 0; i < n + m - 1; i++)
c[i] = -negative_c[i];
return c;
}
/// @brief Min-plus convolution of a convex sequence and an arbitrary sequence.
template <class T>
std::vector<T> convex_arbitrary_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
return arbitrary_convex_minplus_convolution(b, a);
}
namespace max_plus_convolution_internal {
template <class Value>
std::vector<int> monotone_row_minima(int rows, int columns, Value value) {
std::vector<int> result(rows);
int stride = 1;
while (stride < rows) {
stride <<= 1;
}
for (; stride > 0; stride >>= 1) {
for (int row = stride - 1; row < rows; row += 2 * stride) {
int first = row >= stride ? result[row - stride] : 0;
int last = row + stride < rows ? result[row + stride] : columns - 1;
result[row] = first;
for (int column = first + 1; column <= last; column++) {
if (value(row, column) < value(row, result[row])) {
result[row] = column;
}
}
}
}
return result;
}
template <class T>
std::vector<T>
arbitrary_concave_minplus_convolution(const std::vector<T> &arbitrary,
const std::vector<T> &concave) {
if (arbitrary.empty()) {
return concave;
}
if (concave.empty()) {
return arbitrary;
}
for (int i = 0; i + 2 < int(concave.size()); i++) {
assert(concave[i + 1] - concave[i] >=
concave[i + 2] - concave[i + 1]);
}
const int width = int(arbitrary.size());
const int concave_size = int(concave.size());
const int height = width + concave_size - 1;
std::vector<int> minimum_column(height), maximum_column(height, width - 1);
for (int row = concave_size; row < height; row++) {
minimum_column[row] = row - concave_size + 1;
}
for (int row = 0; row <= height - concave_size; row++) {
maximum_column[row] = row;
}
std::vector<int> minimum_row(width), maximum_row(width);
for (int column = 0; column < width; column++) {
minimum_row[column] = column;
maximum_row[column] = concave_size - 1 + column;
}
std::vector<T> result(height, std::numeric_limits<T>::max());
std::function<void(int, int, int, int)> solve =
[&](int first_row, int last_row, int first_column, int last_column) {
if (maximum_column[first_row] >= last_column &&
first_column >= minimum_column[last_row]) {
auto value = [&](int local_row, int reversed_column) {
int column = last_column - reversed_column;
int row = first_row + local_row;
return arbitrary[column] + concave[row - column];
};
auto minima = monotone_row_minima(last_row - first_row + 1,
last_column - first_column + 1,
value);
for (int row = first_row; row <= last_row; row++) {
result[row] = std::min(result[row],
value(row - first_row,
minima[row - first_row]));
}
return;
}
if (std::int64_t(last_row - first_row) *
(last_column - first_column) <
1024) {
for (int row = first_row; row <= last_row; row++) {
int from = std::max(minimum_column[row], first_column);
int to = std::min(maximum_column[row], last_column);
for (int column = from; column <= to; column++) {
result[row] = std::min(
result[row], arbitrary[column] + concave[row - column]);
}
}
return;
}
if (last_row - first_row > last_column - first_column) {
int middle = (first_row + last_row) / 2;
int new_last = std::min(maximum_column[middle], last_column);
if (first_column <= new_last) {
solve(first_row, middle, first_column, new_last);
}
int new_first = std::max(minimum_column[middle], first_column);
if (new_first <= last_column) {
solve(middle + 1, last_row, new_first, last_column);
}
} else {
int middle = (first_column + last_column) / 2;
int new_last = std::min(maximum_row[middle], last_row);
if (first_row <= new_last) {
solve(first_row, new_last, first_column, middle);
}
int new_first = std::max(minimum_row[middle], first_row);
if (new_first <= last_row) {
solve(new_first, last_row, middle + 1, last_column);
}
}
};
solve(0, height - 1, 0, width - 1);
return result;
}
} // namespace max_plus_convolution_internal
/// @brief Min-plus convolution of a concave sequence and an arbitrary
/// sequence. Valid pairs form a diagonal staircase rather than one rectangular
/// Monge matrix. Recursively splitting that staircase produces fully valid
/// rectangles; after reversing their columns, concavity makes row minima
/// monotone and they are found together. Small boundary rectangles are scanned
/// directly.
template <class T>
std::vector<T> concave_arbitrary_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
return max_plus_convolution_internal::arbitrary_concave_minplus_convolution(
b, a);
}
/// @brief Backward-compatible name for arbitrary-convex min-plus convolution.
template <class T>
std::vector<T> concave_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
return arbitrary_convex_minplus_convolution(a, b);
}
} // namespace noya