Implementation of various algorithms
#include "Optimization/smawk.hpp"$H \times W$ 行列 $A$ が Totally Monotone であるとき、 計算量 $O(H + W)$ で、$A$ の各行に対する最小値を求める。
template<class VAL, class FUNC> vector<pair<VAL, int>> SMAWK(int H, int W, const FUNC &f)
res として
res[l].first: $l$ 行目の最小値res[l].second: $l$ 行目の値が最小となる列の番号#pragma once
#include <bits/stdc++.h>
#include <cassert>
using namespace std;
// find min_j f(i, j) for all i, by SMAWK, O(H + W)
// f(i, j) must be totally monotone
// slide access NG
template<class VAL, class FUNC> void SMAWKRec
(const vector<int> &X, const vector<int> &Y, const FUNC &f, vector<pair<VAL, int>> &res) {
if (X.empty()) return;
// Reduce Step
vector<int> X2, Y2;
for (auto y : Y) {
while (!Y2.empty()) {
int py = Y2.back(), x = X[(int)Y2.size() - 1];
if (f(x, y) >= f(x, py)) break;
Y2.pop_back();
}
if (Y2.size() < X.size()) Y2.emplace_back(y);
}
// Recurse Step
for (int i = 1; i < (int)X.size(); i += 2) X2.emplace_back(X[i]);
SMAWKRec(X2, Y2, f, res);
// Interpolate Step
int p = 0;
for (int i = 0; i < (int)X.size(); i += 2) {
int lim = (i + 1 < (int)X.size() ? res[X[i + 1]].second : Y.back()), best = Y[p];
while (Y[p] < lim) {
p++;
if (f(X[i], Y[p]) < f(X[i], best)) best = Y[p];
}
res[X[i]] = {f(X[i], best), best};
}
}
template<class VAL, class FUNC> vector<pair<VAL, int>> SMAWK(int H, int W, const FUNC &f) {
if (H == 0) return {};
assert(W > 0);
vector<pair<VAL, int>> res(H, make_pair(numeric_limits<VAL>::max() / 2, -1));
vector<int> X(H), Y(W);
for (int i = 0; i < H; i++) X[i] = i;
for (int j = 0; j < W; j++) Y[j] = j;
SMAWKRec(X, Y, f, res);
return res;
}#line 2 "Optimization/smawk.hpp"
#include <bits/stdc++.h>
#line 4 "Optimization/smawk.hpp"
using namespace std;
// find min_j f(i, j) for all i, by SMAWK, O(H + W)
// f(i, j) must be totally monotone
// slide access NG
template<class VAL, class FUNC> void SMAWKRec
(const vector<int> &X, const vector<int> &Y, const FUNC &f, vector<pair<VAL, int>> &res) {
if (X.empty()) return;
// Reduce Step
vector<int> X2, Y2;
for (auto y : Y) {
while (!Y2.empty()) {
int py = Y2.back(), x = X[(int)Y2.size() - 1];
if (f(x, y) >= f(x, py)) break;
Y2.pop_back();
}
if (Y2.size() < X.size()) Y2.emplace_back(y);
}
// Recurse Step
for (int i = 1; i < (int)X.size(); i += 2) X2.emplace_back(X[i]);
SMAWKRec(X2, Y2, f, res);
// Interpolate Step
int p = 0;
for (int i = 0; i < (int)X.size(); i += 2) {
int lim = (i + 1 < (int)X.size() ? res[X[i + 1]].second : Y.back()), best = Y[p];
while (Y[p] < lim) {
p++;
if (f(X[i], Y[p]) < f(X[i], best)) best = Y[p];
}
res[X[i]] = {f(X[i], best), best};
}
}
template<class VAL, class FUNC> vector<pair<VAL, int>> SMAWK(int H, int W, const FUNC &f) {
if (H == 0) return {};
assert(W > 0);
vector<pair<VAL, int>> res(H, make_pair(numeric_limits<VAL>::max() / 2, -1));
vector<int> X(H), Y(W);
for (int i = 0; i < H; i++) X[i] = i;
for (int j = 0; j < W; j++) Y[j] = j;
SMAWKRec(X, Y, f, res);
return res;
}