SMAWK Algorithm
(algorithm/smawk.cpp)

• View this file on GitHub
• Last update: 2020-03-03 23:12:43+09:00

入力

$N \times M$ 行列 $A$

• $A$ は totally monotone

出力

各 $i$ について $\mathop {\rm arg \, max} \limits _ j A _ {i , j}$

計算量

• 時間計算量 $\Theta ( N + M )$
• 空間計算量 $\Theta ( N + M )$

説明

$N \times M$ 行列 $A$ について、いずれかの行で最大値を取り得る高々 $N$ 個の列を $\Theta ( N + M )$ で取り出すアルゴリズムが存在する。(reduce)
$1$ つ飛ばしに $N / 2$ 個程度の行を取り出し、reduce することで $N , M$ 共に半分になった問題を解く。
monotone 性から、飛ばした行の $\mathop {\rm arg \, max}$ を復元する。(interpolate)

文献

• Totally Monotone Matrix Searching (SMAWK algorithm)
  週刊 spaghetti_source
  日本語

• The SMAWK Algorithm (pdf)
  具体例が図示されている

• Geometric applications of a matrix-searching algorithm
  SMAWK Algorithm が提案された論文

語

monotone

$A$ が monotone であるとは、任意の $i \lt j$ に対して

$\mathop {\rm arg \, max} \limits _ k A _ {i , k} \le \mathop {\rm arg \, max} \limits _ k A _ {j , k}$

が成立することを指す。
ただし適切な tie-break により $\mathop {\rm arg \, max}$ は $1$ つに定まるとする。
$\mathop {\rm arg \, min}$ による定義も存在するので注意。

totally monotone

$A$ が totally monotone であるとは、その任意の部分行列が monotone であることを指す。
同値な定義に「$A$ の任意の $2 \times 2$ 部分行列が monotone」がある。

その他

• 実装例の引数 $\rm select$ の定義
  ${\rm select} ( i , j , k ) \mathrel := A _ {i , j}$ と $A _ {i , k} \ ( j \lt k )$ を比較して、前者を選ぶ場合 false、後者を選ぶ場合 true

Depends on

• other/int_alias.cpp

Required by

• Axiotis-Tzamos Knapsack (algorithm/axiotis_tzamos_knapsack.cpp)
• Concave Max Plus Convolution (algorithm/concave_max_plus_convolution.cpp)

Verified with

• test/axiotis_tzamos_knapsack.test.cpp

Code

#include "other/int_alias.cpp"

#include <functional>
#include <numeric>
#include <vector>

template <class Select>
std::vector<usize> smawk(const usize row_size, const usize col_size,
                         const Select &select) {
  using vec_zu = std::vector<usize>;

  const std::function<vec_zu(const vec_zu &, const vec_zu &)> solve =
      [&](const vec_zu &row, const vec_zu &col) -> vec_zu {
    const usize n = row.size();
    if (n == 0)
      return {};
    vec_zu c2;
    for (const usize i : col) {
      while (!c2.empty() && select(row[c2.size() - 1], c2.back(), i))
        c2.pop_back();
      if (c2.size() < n)
        c2.push_back(i);
    }
    vec_zu r2;
    for (usize i = 1; i < n; i += 2)
      r2.push_back(row[i]);
    const vec_zu a2 = solve(r2, c2);
    vec_zu ans(n);
    for (usize i = 0; i != a2.size(); i += 1)
      ans[i * 2 + 1] = a2[i];
    usize j = 0;
    for (usize i = 0; i < n; i += 2) {
      ans[i] = c2[j];
      const usize 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;
  };
  vec_zu row(row_size);
  std::iota(row.begin(), row.end(), 0);
  vec_zu col(col_size);
  std::iota(col.begin(), col.end(), 0);
  return solve(row, col);
}

/**
 * @brief SMAWK Algorithm
 * @docs docs/smawk.md
 */

#line 2 "other/int_alias.cpp"

#include <cstddef>
#include <cstdint>

using i32 = std::int32_t;
using i64 = std::int64_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;
using isize = std::ptrdiff_t;
using usize = std::size_t;
#line 2 "algorithm/smawk.cpp"

#include <functional>
#include <numeric>
#include <vector>

template <class Select>
std::vector<usize> smawk(const usize row_size, const usize col_size,
                         const Select &select) {
  using vec_zu = std::vector<usize>;

  const std::function<vec_zu(const vec_zu &, const vec_zu &)> solve =
      [&](const vec_zu &row, const vec_zu &col) -> vec_zu {
    const usize n = row.size();
    if (n == 0)
      return {};
    vec_zu c2;
    for (const usize i : col) {
      while (!c2.empty() && select(row[c2.size() - 1], c2.back(), i))
        c2.pop_back();
      if (c2.size() < n)
        c2.push_back(i);
    }
    vec_zu r2;
    for (usize i = 1; i < n; i += 2)
      r2.push_back(row[i]);
    const vec_zu a2 = solve(r2, c2);
    vec_zu ans(n);
    for (usize i = 0; i != a2.size(); i += 1)
      ans[i * 2 + 1] = a2[i];
    usize j = 0;
    for (usize i = 0; i < n; i += 2) {
      ans[i] = c2[j];
      const usize 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;
  };
  vec_zu row(row_size);
  std::iota(row.begin(), row.end(), 0);
  vec_zu col(col_size);
  std::iota(col.begin(), col.end(), 0);
  return solve(row, col);
}

/**
 * @brief SMAWK Algorithm
 * @docs docs/smawk.md
 */

Back to top page