max_plus_convolution.hpp¶
#include "noya/max_plus_convolution.hpp"
#ifndef NOYA_MAXPLUS_CONVOLUTION_HPP
#define NOYA_MAXPLUS_CONVOLUTION_HPP 1
#include "noya/smawk.hpp"
#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 concave sequences.
template <class T>
std::vector<T> two_concave_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
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 Max-plus convolution where b is concave.
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 where b is concave.
template <class T>
std::vector<T> concave_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
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;
}
} // namespace noya
#endif // NOYA_MAXPLUS_CONVOLUTION_HPP
#include <functional>
#include <numeric>
#include <vector>
// https://noshi91.github.io/Library/algorithm/smawk.cpp.html
namespace noya {
/// @brief SMAWK algorithm: compute row minima of a totally monotone matrix.
/// @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 concave sequences.
template <class T>
std::vector<T> two_concave_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
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 Max-plus convolution where b is concave.
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 where b is concave.
template <class T>
std::vector<T> concave_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
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;
}
} // namespace noya