Implementation of various algorithms
#include "Optimization/dd_smawk.hpp"SMAWK 法に対して、オンライン・オフライン変換を施して得られる分割統治 SMAWK 法。単一始点最短路問題を解く。
DP が
dp[j] = min_{i < j} (dp[i] + func(i, j))
の形で表されるとする。ここで、
g(j, i) = dp[i] + func(i, j)
が Totally Monotone であるならば、上記の DP が分割統治 SMAWK 法によって、$O(N \log N)$ の計算量で DP テーブルを求めることができる。
func(i, j) が Totally Monotone であっても、g(j, i) = dp[i] + func(i, j) が Totally Monotone であるとは限らないfunc(i, j) が Monge であるとき、g(j, i) = dp[i] + func(i, j) も Monge である上記のことから、多くの場合、コストが Monge 性を満たすときに適用できるアルゴリズムであるとされる。
template<class FUNC> vector<pair<VAL, int>> solve(int N, const FUNC &f)
構造体 DDSMAWK<VAL> を宣言して、関数 solve(N, f) を呼び出す。
N: 頂点番号が $0, 1, 2, \dots, N$ であるとする (頂点数は $N + 1$)f(l, r): 辺 $(l, r)$ のコストを返す関数res として
res[v].first: 頂点 $0$ から頂点 $v$ への最短路の重みres[v].second: 頂点 $0$ から頂点 $v$ への最短路において、$v$ に至る直前の頂点の番号 (存在しない場合は -1)dp[v]: 頂点 $0$ から頂点 $v$ への最短路の重みcnt[v]: 頂点 $0$ から頂点 $v$ への最短路に含まれる辺の本数prev[v]: 頂点 $0$ から頂点 $v$ への最短路において、$v$ に至る直前の頂点の番号 (存在しない場合は -1)vector<int> reconstruct()
solve(N, f) 実行後に呼び出すことで、頂点 $0$ から頂点 $N$ への最短路が復元される#pragma once
#include <bits/stdc++.h>
#include <cassert>
using namespace std;
#include "smawk.hpp"
// find shortest path on DAG by D&D SMAWK, O(N log N)
// vertex: 0, 1, 2, ..., N
// A(j, i) = dp[i] + f(i, j) must be totally monotone
// Attention: Even if f is totally monotone, A is not necessarily totally monotone
// slide access NG
template<class VAL> struct DDSMAWK {
VAL INF = numeric_limits<VAL>::max() / 2;
int CNT_INF = numeric_limits<int>::max() / 2;
// results
vector<VAL> dp;
vector<int> cnt, prev;
// solver
template<class FUNC> vector<pair<VAL, int>> solve(int N, const FUNC &f) {
dp.assign(N + 1, INF);
cnt.assign(N + 1, CNT_INF);
prev.assign(N + 1, 0);
vector<pair<VAL, int>> res(N + 1, make_pair(numeric_limits<VAL>::max() / 2, -1));
vector<pair<VAL, int>> tmp(N + 1);
dp[0] = 0, cnt[0] = 0, res[0].first = VAL(0);
auto f2 = [&](int i, int j) -> VAL { return res[j].first + f(j, i); };
auto rec = [&](auto &&rec, int left, int right) -> void {
if (right - left <= 1) return;
int mid = (left + right) / 2;
vector<int> X(right - mid), Y(mid - left);
for (int i = mid; i < right; i++) X[i - mid] = i;
for (int j = left; j < mid; j++) Y[j - left] = j;
rec(rec, left, mid);
SMAWKRec(X, Y, f2, tmp);
for (auto x : X) if (tmp[x].first < res[x].first) res[x] = tmp[x];
rec(rec, mid, right);
};
rec(rec, 0, N + 1);
for (int r = 1; r <= N; r++) {
dp[r] = res[r].first;
prev[r] = res[r].second;
if (prev[r] != -1) cnt[r] = cnt[prev[r]] + 1;
}
return res;
}
// reconstruct the shortest path
vector<int> reconstruct() {
int N = (int)dp.size() - 1;
vector<int> path;
for (int v = N; v > 0; v = prev[v]) path.emplace_back(v);
path.emplace_back(0);
reverse(path.begin(), path.end());
return path;
}
};#line 2 "Optimization/dd_smawk.hpp"
#include <bits/stdc++.h>
#line 4 "Optimization/dd_smawk.hpp"
using namespace std;
#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;
}
#line 7 "Optimization/dd_smawk.hpp"
// find shortest path on DAG by D&D SMAWK, O(N log N)
// vertex: 0, 1, 2, ..., N
// A(j, i) = dp[i] + f(i, j) must be totally monotone
// Attention: Even if f is totally monotone, A is not necessarily totally monotone
// slide access NG
template<class VAL> struct DDSMAWK {
VAL INF = numeric_limits<VAL>::max() / 2;
int CNT_INF = numeric_limits<int>::max() / 2;
// results
vector<VAL> dp;
vector<int> cnt, prev;
// solver
template<class FUNC> vector<pair<VAL, int>> solve(int N, const FUNC &f) {
dp.assign(N + 1, INF);
cnt.assign(N + 1, CNT_INF);
prev.assign(N + 1, 0);
vector<pair<VAL, int>> res(N + 1, make_pair(numeric_limits<VAL>::max() / 2, -1));
vector<pair<VAL, int>> tmp(N + 1);
dp[0] = 0, cnt[0] = 0, res[0].first = VAL(0);
auto f2 = [&](int i, int j) -> VAL { return res[j].first + f(j, i); };
auto rec = [&](auto &&rec, int left, int right) -> void {
if (right - left <= 1) return;
int mid = (left + right) / 2;
vector<int> X(right - mid), Y(mid - left);
for (int i = mid; i < right; i++) X[i - mid] = i;
for (int j = left; j < mid; j++) Y[j - left] = j;
rec(rec, left, mid);
SMAWKRec(X, Y, f2, tmp);
for (auto x : X) if (tmp[x].first < res[x].first) res[x] = tmp[x];
rec(rec, mid, right);
};
rec(rec, 0, N + 1);
for (int r = 1; r <= N; r++) {
dp[r] = res[r].first;
prev[r] = res[r].second;
if (prev[r] != -1) cnt[r] = cnt[prev[r]] + 1;
}
return res;
}
// reconstruct the shortest path
vector<int> reconstruct() {
int N = (int)dp.size() - 1;
vector<int> path;
for (int v = N; v > 0; v = prev[v]) path.emplace_back(v);
path.emplace_back(0);
reverse(path.begin(), path.end());
return path;
}
};