convex/minplus_convolution.hpp

• View this file on GitHub
• Last update: 2023-08-10 03:11:20+09:00
• Include: #include "convex/minplus_convolution.hpp"

Depends on

• convex/monotone_minima.hpp

Required by

• convex/minplus_convolution_of_triples.hpp

Verified with

• test/library_checker/convolution/minplus_convolution_conv_arb.test.cpp
• test/library_checker/convolution/minplus_convolution_conv_conv.test.cpp
• test/mytest/minplus_conv_triple.test.cpp
• test/mytest/minplus_convex.test.cpp

Code

#include "convex/monotone_minima.hpp"

template <typename T>
vc<T> minplus_convolution_convex_convex(vc<T>& A, vc<T>& B) {
  int n = len(A), m = len(B);
  if (n == 0 && m == 0) return {};
  vc<T> C(n + m - 1, infty<T>);
  while (n > 0 && A[n - 1] == infty<T>) --n;
  while (m > 0 && B[m - 1] == infty<T>) --m;
  if (n == 0 && m == 0) return C;
  int a = 0, b = 0;
  while (a < n && A[a] == infty<T>) ++a;
  while (b < m && B[b] == infty<T>) ++b;
  C[a + b] = A[a] + B[b];
  for (int i = a + b + 1; i < n + m - 1; ++i) {
    if (b == m - 1 || (a != n - 1 && A[a + 1] + B[b] < A[a] + B[b + 1])) {
      chmin(C[i], A[++a] + B[b]);
    } else {
      chmin(C[i], A[a] + B[++b]);
    }
  }
  return C;
}

template <typename T>
vc<T> minplus_convolution_arbitrary_convex(vc<T>& A, vc<T>& B) {
  int n = len(A), m = len(B);
  if (n == 0 && m == 0) return {};
  vc<T> C(n + m - 1, infty<T>);
  while (m > 0 && B[m - 1] == infty<T>) --m;
  if (m == 0) return C;
  int b = 0;
  while (b < m && B[b] == infty<T>) ++b;

  auto select = [&](int i, int j, int k) -> bool {
    if (i < k) return false;
    if (i - j >= m - b) return true;
    return A[j] + B[b + i - j] >= A[k] + B[b + i - k];
  };
  vc<int> J = monotone_minima(n + m - b - 1, n, select);
  FOR(i, n + m - b - 1) {
    T x = A[J[i]], y = B[b + i - J[i]];
    if (x < infty<T> && y < infty<T>) C[b + i] = x + y;
  }
  return C;
}

template <typename T, bool convA, bool convB>
vc<T> minplus_convolution(vc<T>& A, vc<T>& B) {
  static_assert(convA || convB);
  if constexpr (convA && convB) return minplus_convolution_convex_convex(A, B);
  if constexpr (convA && !convB)
    return minplus_convolution_arbitrary_convex(B, A);
  if constexpr (convB && !convA)
    return minplus_convolution_arbitrary_convex(A, B);
  return {};
}

#line 1 "convex/monotone_minima.hpp"

// select(i,j,k) : (i,j) -> (i,k) を行うかどうか
template <typename F>
vc<int> monotone_minima(int H, int W, F select) {
  vc<int> min_col(H);
  auto dfs = [&](auto& dfs, int x1, int x2, int y1, int y2) -> void {
    if (x1 == x2) return;
    int x = (x1 + x2) / 2;
    int best_y = y1;
    for (int y = y1 + 1; y < y2; ++y) {
      if (select(x, best_y, y)) best_y = y;
    }
    min_col[x] = best_y;
    dfs(dfs, x1, x, y1, best_y + 1);
    dfs(dfs, x + 1, x2, best_y, y2);
  };
  dfs(dfs, 0, H, 0, W);
  return min_col;
}
#line 2 "convex/minplus_convolution.hpp"

template <typename T>
vc<T> minplus_convolution_convex_convex(vc<T>& A, vc<T>& B) {
  int n = len(A), m = len(B);
  if (n == 0 && m == 0) return {};
  vc<T> C(n + m - 1, infty<T>);
  while (n > 0 && A[n - 1] == infty<T>) --n;
  while (m > 0 && B[m - 1] == infty<T>) --m;
  if (n == 0 && m == 0) return C;
  int a = 0, b = 0;
  while (a < n && A[a] == infty<T>) ++a;
  while (b < m && B[b] == infty<T>) ++b;
  C[a + b] = A[a] + B[b];
  for (int i = a + b + 1; i < n + m - 1; ++i) {
    if (b == m - 1 || (a != n - 1 && A[a + 1] + B[b] < A[a] + B[b + 1])) {
      chmin(C[i], A[++a] + B[b]);
    } else {
      chmin(C[i], A[a] + B[++b]);
    }
  }
  return C;
}

template <typename T>
vc<T> minplus_convolution_arbitrary_convex(vc<T>& A, vc<T>& B) {
  int n = len(A), m = len(B);
  if (n == 0 && m == 0) return {};
  vc<T> C(n + m - 1, infty<T>);
  while (m > 0 && B[m - 1] == infty<T>) --m;
  if (m == 0) return C;
  int b = 0;
  while (b < m && B[b] == infty<T>) ++b;

  auto select = [&](int i, int j, int k) -> bool {
    if (i < k) return false;
    if (i - j >= m - b) return true;
    return A[j] + B[b + i - j] >= A[k] + B[b + i - k];
  };
  vc<int> J = monotone_minima(n + m - b - 1, n, select);
  FOR(i, n + m - b - 1) {
    T x = A[J[i]], y = B[b + i - J[i]];
    if (x < infty<T> && y < infty<T>) C[b + i] = x + y;
  }
  return C;
}

template <typename T, bool convA, bool convB>
vc<T> minplus_convolution(vc<T>& A, vc<T>& B) {
  static_assert(convA || convB);
  if constexpr (convA && convB) return minplus_convolution_convex_convex(A, B);
  if constexpr (convA && !convB)
    return minplus_convolution_arbitrary_convex(B, A);
  if constexpr (convB && !convA)
    return minplus_convolution_arbitrary_convex(A, B);
  return {};
}

Back to top page