Implementation of various algorithms
#include "Optimization/monge_enum_d_shortest_paths.hpp"であるような DAG を考える。$f(l, r)$ が Monge であるとする。
このとき、$d = 1, 2, \dots, D$ に対して、頂点 $0$ から各頂点 $v$ への最短パスのうち、パスに含まれる辺の本数がちょうど $d$ 本であるようなものを求める。
次の DP を高速化するアルゴリズムであるとも言える。
dp[j][d] = min_{i < j} (dp[i][d - 1] + f(i, j))
これは通常 $O(N^2 D)$ の計算量を要するが、$f$ が Monge ならば高速に実行できる。
template<class FUNC> vector<vector<pair<VAL, int>>> solve(int N, const FUNC &f, int D, const string &solver = "smawk")
N: 頂点番号が $0, 1, 2, \dots, N$ であるとする (頂点数は $N + 1$)f(l, r): 辺 $(l, r)$ のコストを返す関数solver: “smawk” (SMAWK 法を使用) か “monotone” (Monotone Minima 法を使用) を指定するdp として
dp[v][d].first: 頂点 $0$ から頂点 $v$ への、ちょうど $d$ 本の辺を用いる最短路の重みdp[v][d].second: 頂点 $0$ から頂点 $v$ への、ちょうど $d$ 本の辺を用いる最短路において、$v$ に至る直前の頂点の番号 (存在しない場合は -1)template<class ADD, class DEL, class GETCOST> vector<vector<pair<VAL, int>>> solve
(int N, const ADD &add, const DEL &del, const GETCOST &get, int D, const string &solver = "monotone")
N: 頂点番号が $0, 1, 2, \dots, N$ であるとする (頂点数は $N + 1$)add(i): 辺 $(l, r)$ のコストが、区間 $[l, r)$ についての値に対応するとして、区間に index $i$ の分を追加するときのコストの差分を記述する関数del(i): 辺 $(l, r)$ のコストが、区間 $[l, r)$ についての値に対応するとして、区間から index $i$ の分を削除するときのコストの差分を記述する関数get(): 辺 $(l, r)$ のコストを返す関数solver: “monotone” (Monotone Minima 法を使用) を指定する (“smawk” は不可)dp として
dp[v][d].first: 頂点 $0$ から頂点 $v$ への、ちょうど $d$ 本の辺を用いる最短路の重みdp[v][d].second: 頂点 $0$ から頂点 $v$ への、ちょうど $d$ 本の辺を用いる最短路において、$v$ に至る直前の頂点の番号 (存在しない場合は -1)#pragma once
#include <bits/stdc++.h>
#include <cassert>
using namespace std;
#include "monotone_minima.hpp"
#include "smawk.hpp"
// find the d-edges shortest path for d = 1, 2, ..., D, in O(ND)
// vertex: 0, 1, 2, ..., N
// f(i, j) must be Monge (Then, A(j, i) = dp[i] + f(i, j) is also Monge)
// you can choose SMAWK (random access) or Monotone Minima (slide access)
template<class VAL> struct MongeShortestPathWithDEdges {
VAL INF = numeric_limits<VAL>::max() / 2;
// solver (random access ver)
template<class FUNC> vector<vector<pair<VAL, int>>> solve(int N, const FUNC &f, int D, const string &solver = "smawk") {
assert(D <= N);
vector<vector<pair<VAL, int>>> res(N + 1, vector<pair<VAL, int>>(D + 1, make_pair(INF, -1)));
vector<VAL> dp(N + 1, INF);
res[0][0].first = 0, dp[0] = 0;
auto f2 = [&](int i, int j) -> VAL { return (j < i ? dp[j] + f(j, i) : INF); };
for (int d = 1; d <= D; d++) {
vector<pair<VAL, int>> tmp;
if (solver == "smawk") tmp = SMAWK<VAL>(N + 1, N + 1, f2);
else if (solver == "monotone") tmp = MonotoneMinima<VAL>(N + 1, N + 1, f2);
else {
cout << "no solver" << endl;
assert(true);
}
for (int i = d; i <= N; i++) {
dp[i] = tmp[i].first;
res[i][d] = tmp[i];
}
}
return res;
}
// solver (slide access ver)
template<class ADD, class DEL, class GETCOST> vector<vector<pair<VAL, int>>> solve
(int N, const ADD &add, const DEL &del, const GETCOST &get, int D, const string &solver = "monotone") {
assert(D <= N);
vector<vector<pair<VAL, int>>> res(N + 1, vector<pair<VAL, int>>(D + 1, make_pair(INF, -1)));
vector<VAL> dp(N + 1, INF);
res[0][0].first = 0, dp[0] = 0;
int left = 0, right = 0;
auto f = [&](int i, int j) -> long long {
while (left < i) del(left++);
while (i < left) add(--left);
while (right < j) add(right++);
while (j < right) del(--right);
return get();
};
auto f2 = [&](int i, int j) -> VAL { return (j < i ? dp[j] + f(j, i) : INF); };
for (int d = 1; d <= D; d++) {
vector<pair<VAL, int>> tmp;
if (solver == "monotone") tmp = MonotoneMinima<VAL>(N + 1, N + 1, f2);
else {
cout << "no solver" << endl;
assert(true);
}
for (int i = d; i <= N; i++) {
dp[i] = tmp[i].first;
res[i][d] = tmp[i];
}
}
return res;
}
};#line 2 "Optimization/monge_enum_d_shortest_paths.hpp"
#include <bits/stdc++.h>
#line 4 "Optimization/monge_enum_d_shortest_paths.hpp"
using namespace std;
#line 4 "Optimization/monotone_minima.hpp"
using namespace std;
// find min_j f(i, j) for all i, by Monotone Minima, O(H + W log H)
// f(i, j) must be monotone (argmin is not decreasing)
// slide access OK
template<class VAL, class FUNC> void MonotoneMinimaRec
(int HL, int HR, int WL, int WR, const FUNC &f, vector<pair<VAL, int>> &res) {
if (HR - HL <= 0) return;
int HM = (HL + HR) / 2;
res[HM].second = WL;
for (int i = WL; i < WR; i++) {
VAL val = f(HM, i);
if (res[HM].first > val) res[HM] = make_pair(val, i);
}
MonotoneMinimaRec(HL, HM, WL, res[HM].second + 1, f, res);
MonotoneMinimaRec(HM + 1, HR, res[HM].second, WR, f, res);
}
template<class VAL, class FUNC> vector<pair<VAL, int>> MonotoneMinima(int H, int W, const FUNC &f) {
vector<pair<VAL, int>> res(H, make_pair(numeric_limits<VAL>::max() / 2, -1));
MonotoneMinimaRec(0, H, 0, W, f, res);
return res;
}
#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 8 "Optimization/monge_enum_d_shortest_paths.hpp"
// find the d-edges shortest path for d = 1, 2, ..., D, in O(ND)
// vertex: 0, 1, 2, ..., N
// f(i, j) must be Monge (Then, A(j, i) = dp[i] + f(i, j) is also Monge)
// you can choose SMAWK (random access) or Monotone Minima (slide access)
template<class VAL> struct MongeShortestPathWithDEdges {
VAL INF = numeric_limits<VAL>::max() / 2;
// solver (random access ver)
template<class FUNC> vector<vector<pair<VAL, int>>> solve(int N, const FUNC &f, int D, const string &solver = "smawk") {
assert(D <= N);
vector<vector<pair<VAL, int>>> res(N + 1, vector<pair<VAL, int>>(D + 1, make_pair(INF, -1)));
vector<VAL> dp(N + 1, INF);
res[0][0].first = 0, dp[0] = 0;
auto f2 = [&](int i, int j) -> VAL { return (j < i ? dp[j] + f(j, i) : INF); };
for (int d = 1; d <= D; d++) {
vector<pair<VAL, int>> tmp;
if (solver == "smawk") tmp = SMAWK<VAL>(N + 1, N + 1, f2);
else if (solver == "monotone") tmp = MonotoneMinima<VAL>(N + 1, N + 1, f2);
else {
cout << "no solver" << endl;
assert(true);
}
for (int i = d; i <= N; i++) {
dp[i] = tmp[i].first;
res[i][d] = tmp[i];
}
}
return res;
}
// solver (slide access ver)
template<class ADD, class DEL, class GETCOST> vector<vector<pair<VAL, int>>> solve
(int N, const ADD &add, const DEL &del, const GETCOST &get, int D, const string &solver = "monotone") {
assert(D <= N);
vector<vector<pair<VAL, int>>> res(N + 1, vector<pair<VAL, int>>(D + 1, make_pair(INF, -1)));
vector<VAL> dp(N + 1, INF);
res[0][0].first = 0, dp[0] = 0;
int left = 0, right = 0;
auto f = [&](int i, int j) -> long long {
while (left < i) del(left++);
while (i < left) add(--left);
while (right < j) add(right++);
while (j < right) del(--right);
return get();
};
auto f2 = [&](int i, int j) -> VAL { return (j < i ? dp[j] + f(j, i) : INF); };
for (int d = 1; d <= D; d++) {
vector<pair<VAL, int>> tmp;
if (solver == "monotone") tmp = MonotoneMinima<VAL>(N + 1, N + 1, f2);
else {
cout << "no solver" << endl;
assert(true);
}
for (int i = d; i <= N; i++) {
dp[i] = tmp[i].first;
res[i][d] = tmp[i];
}
}
return res;
}
};