凸・凹 min-plus / max-plus 畳み込み

convolution/min_plus_convolution.hpp · GitHub で見る

共通ページ

#include "lib/convolution/min_plus_convolution.hpp"

std::vector<long long> a{3, 1, 0, 3}, b{5, 1, 3, 3, 2};
auto c = min_plus_convolution_convex(a, b); // {8,4,2,1,3,3,2,5}

min-plus は \(\displaystyle \min_{i+j=k}(a_i+b_j)\)、max-plus は \(\displaystyle \max_{i+j=k}(a_i+b_j)\) を計算します。

関数	\(a\) の条件	時間計算量
min_plus_convolution_convex(a,b)	隣接差が単調非減少	\(O(L)\)
max_plus_convolution_concave(a,b)	隣接差が単調非増加	\(O(L)\)
min_plus_convolution_concave(a,b)	隣接差が単調非増加	\(O(L \log L)\)
max_plus_convolution_convex(a,b)	隣接差が単調非減少	\(O(L \log L)\)

凸か凹かの条件は第 1 引数 a に必要です。例えば上の a の隣接差は {-2,-1,3} で、順に小さくならないので凸です。条件を満たす列を第 1 引数に渡します。両方が非空なら返り値は長さ a.size()+b.size()-1 の配列、どちらかが空なら空配列です。入力は変更しません。

\(b\) は任意です。要素と和は有限で、係数型に収まるものとします。\(L = \lvert a\rvert + \lvert b\rvert\) として、追加空間はいずれも \(O(L)\) です。

依存ライブラリ

• 凹列と任意列の畳み込み（補助部品）
• min-plus / max-plus 畳み込みの整数用共通実装
• SMAWK

検証コード

• test/min_plus_concave.test.cpp
• test/min_plus_convex.test.cpp
• test/min_plus_convex_convex.test.cpp

コード

#ifndef MIN_PLUS_CONVOLUTION_HPP
#define MIN_PLUS_CONVOLUTION_HPP

#include "../components/smawk.hpp"
#include "../components/min_plus_convolution_integral.hpp"
#include "../components/min_plus_convolution_envelope.hpp"
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <functional>
#include <limits>
#include <numeric>
#include <type_traits>
#include <utility>
#include <vector>

namespace min_plus_convolution_detail {

template <class T, class Compare>
std::vector<T> multiply(const std::vector<T> &a, const std::vector<T> &b, Compare compare) {
	if (a.empty() || b.empty()) return {};
	assert(a.size() <= std::size_t(std::numeric_limits<int>::max()) && b.size() <= std::size_t(std::numeric_limits<int>::max()) - a.size());
	if constexpr (min_plus_convolution_integral_detail::supported<T>
		&& (std::is_same_v<Compare, std::less<T>> || std::is_same_v<Compare, std::greater<T>>))
		if (auto result = min_plus_convolution_integral_detail::multiply<T, std::is_same_v<Compare, std::less<T>>>(a, b)) return std::move(*result);
	int n = a.size(), m = b.size(), length = n + m - 1;
	std::vector<int> columns(m);
	std::iota(columns.begin(), columns.end(), 0);
	auto better = [&](int row, int x, int y) {
		if (row < y) return false;
		if (row - x >= n) return true;
		return compare(a[row - y] + b[y], a[row - x] + b[x]);
	};
	auto best = smawk(length, columns, better);
	std::vector<T> result;
	result.reserve(length);
	for (int i = 0; i < length; ++i) result.push_back(a[i - best[i]] + b[best[i]]);
	return result;
}
template <class T, class Compare>
std::vector<T> opposite(const std::vector<T> &a, const std::vector<T> &b, Compare compare) {
	if (a.empty() || b.empty()) return {};
	assert(a.size() <= std::size_t(std::numeric_limits<int>::max()) && b.size() <= std::size_t(std::numeric_limits<int>::max()) - a.size());
	if constexpr (min_plus_convolution_integral_detail::supported<T>
		&& (std::is_same_v<Compare, std::less<T>> || std::is_same_v<Compare, std::greater<T>>))
		if (auto result = min_plus_convolution_envelope_detail::multiply<T, std::is_same_v<Compare, std::less<T>>>(a, b)) return std::move(*result);
	int n = a.size(), m = b.size(), length = n + m - 1;
	std::vector<T> result(length, a[0] + b[0]);
	std::vector<std::uint8_t> done(length);
	auto solve = [&](auto &&self, int l, int r, int parent_l, int parent_r) -> void {
		int low = std::max(0, r - n), high = std::min(m, l + 1);
		std::vector<int> columns;
		for (int i = high - 1; i >= low;) {
			if (parent_l <= i && i < parent_r) i = parent_l - 1;
			else columns.push_back(i--);
		}
		if (!columns.empty()) {
			auto best = smawk(r - l, columns, [&](int row, int x, int y) {
				return compare(a[l + row - y] + b[y], a[l + row - x] + b[x]);
			});
			for (int i = l; i < r; ++i) {
				T value = a[i - best[i - l]] + b[best[i - l]];
				if (!done[i] || compare(value, result[i])) result[i] = value;
				done[i] = 1;
			}
		}
		if (r - l > 1) {
			int mid = l + (r - l) / 2;
			self(self, l, mid, low, high);
			self(self, mid, r, low, high);
		}
	};
	solve(solve, 0, length, 0, 0);
	return result;
}

}

template <class T>
std::vector<T> min_plus_convolution_convex(const std::vector<T> &a, const std::vector<T> &b) {
	return min_plus_convolution_detail::multiply(a, b, std::less<T>{});
}
template <class T>
std::vector<T> max_plus_convolution_concave(const std::vector<T> &a, const std::vector<T> &b) {
	return min_plus_convolution_detail::multiply(a, b, std::greater<T>{});
}

template <class T>
std::vector<T> min_plus_convolution_concave(const std::vector<T> &a, const std::vector<T> &b) {
	return min_plus_convolution_detail::opposite(a, b, std::less<T>{});
}
template <class T>
std::vector<T> max_plus_convolution_convex(const std::vector<T> &a, const std::vector<T> &b) {
	return min_plus_convolution_detail::opposite(a, b, std::greater<T>{});
}

#endif

Back to top page