Min Plus Convolution (Concave and Arbitary) (dp/min-plus-convolution-concave-arbitary.hpp)

• View this file on GitHub
• View document part on GitHub
• Last update: 2026-08-23 09:03:28+00:00
• Include: #include "dp/min-plus-convolution-concave-arbitary.hpp"

凹数列と任意の数列の min-plus 畳み込みを、分割統治と SMAWK により準線形時間で計算する。有効な添字の領域を長方形に分割し、各長方形が表す全単調行列の行最小値を SMAWK で求める。

min_plus_convolution_concave_arbitary

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)$

Depends on

• SMAWK (dp/smawk.hpp)

Verified with

• test/verify/yosupo-min-plus-convolution-concave-arbitrary.test.cpp

Code

#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;
}

#line 2 "dp/min-plus-convolution-concave-arbitary.hpp"

#include <algorithm>
#include <limits>
#include <vector>

#line 2 "dp/smawk.hpp"

#line 4 "dp/smawk.hpp"
#include <numeric>
#line 6 "dp/smawk.hpp"

template <typename F>
std::vector<int> smawk(int H, int W, F comp) {
  std::vector<int> ret(H, -1);
  if (H == 0 || W == 0) return ret;

  auto dfs = [&](auto&& self, const std::vector<int>& rows,
                 const std::vector<int>& cols) -> void {
    if (rows.empty()) return;
    std::vector<int> 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<int> odd_rows;
    odd_rows.reserve(rows.size() / 2);
    for (int i = 1; i < static_cast<int>(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<int>(rows.size()); i += 2) {
      int right = static_cast<int>(reduced.size()) - 1;
      if (i + 1 < static_cast<int>(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<int> rows(H), cols(W);
  std::iota(rows.begin(), rows.end(), 0);
  std::iota(cols.begin(), cols.end(), 0);
  dfs(dfs, rows, cols);
  return ret;
}
#line 8 "dp/min-plus-convolution-concave-arbitary.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;
}

Back to top page