algorithm
Implementation of various algorithms
Monge 全部乗せ
(Optimization/monge_all.hpp)
- View this file on GitHub
- Last update: 2026-11-12 11:56:29+09:00
- Include:
#include "Optimization/monge_all.hpp"
概要
Monge に関するアルゴリズムをまとめたもの。
参考文献
- tatyam: Monge の手引き書
Depends on
-
Monge 単一始点 d-辺最短路問題 (by Alien DP)
(Optimization/alien_dp.hpp)
-
Alien Trick --- min f(x) s.t. g(x) = K (by ラグランジュ緩和)
(Optimization/alien_trick.hpp)
-
anti-Monge 単一始点最短路問題 (by D&D SMAWK 法, in O(N log N))
(Optimization/anti_monge_dd_smawk.hpp)
-
Monotone 単一始点最短路問題 (D&D Monotone Minima, in O(N (log N)^2))
(Optimization/dd_monotone_minima.hpp)
-
TM 単一始点最短路問題 (D&D SMAWK 法, in O(N log N))
(Optimization/dd_smawk.hpp)
-
General Li Chao Tree (一般関数対応)
(Optimization/general_li_chao_tree.hpp)
-
Monge 性 & 単調性を満たす区間 DP の Knuth-Yao Speedup (in O(N^2))
(Optimization/knuth_yao.hpp)
-
Monge 単一始点最短路問題 (by LARSCH 法, in O(N))
(Optimization/larsch.hpp)
-
Monge 性・anti-Monge 性の判定 (in O(HW))
(Optimization/monge_checker.hpp)
-
Monge 単一始点 d-辺最短路の d = 1, 2, ..., D における列挙
(Optimization/monge_enum_d_shortest_paths.hpp)
-
Monge 行列の min-plus 合成 (in O(N^2))
(Optimization/monge_min_plus_composition.hpp)
-
Monge 単一始点最短路アルゴリズム群のラッパー
(Optimization/monge_shortest_path.hpp)
-
Monge 単一始点最短路問題 (by Convex Hull Trick, in O(N log N))
(Optimization/monge_shortest_path_by_cht.hpp)
-
Monotone 行最小値問題 (Monotone Minima, in O(H + W log H))
(Optimization/monotone_minima.hpp)
-
Monge 単一始点最短路問題 (by noshi's 簡易 LARSCH 法, in O(N log N))
(Optimization/simple_larsch.hpp)
-
TM 行最小値問題 (SMAWK 法, in O(H + W))
(Optimization/smawk.hpp)
Code
#pragma once
#include <bits/stdc++.h>
#include <cassert>
using namespace std;
#include "monge_checker.hpp"
#include "monge_shortest_path.hpp"
#include "monge_shortest_path_by_cht.hpp"
#include "monge_enum_d_shortest_paths.hpp"
#include "general_li_chao_tree.hpp"
#include "alien_dp.hpp"
#include "anti_monge_dd_smawk.hpp"
#include "knuth_yao.hpp"
#include "monge_min_plus_composition.hpp"#line 2 "Optimization/monge_all.hpp"
#include <bits/stdc++.h>
#line 4 "Optimization/monge_all.hpp"
using namespace std;
#line 4 "Optimization/monge_checker.hpp"
using namespace std;
// check whether N x M matrix f is Monge or not, in O(NM)
template<class FUNC> bool is_monge
(int N, int M, const FUNC &f, bool upper_triangular = false, bool output_detail = true) {
assert(N > 0 && M > 0);
for (int i = 0; i + 1 < N; i++) {
for (int j = (upper_triangular ? i+2 : 0); j + 1 < M; j++) {
auto f00 = f(i, j), f01 = f(i, j + 1), f10 = f(i + 1, j), f11 = f(i + 1, j + 1);
if (f00 + f11 > f10 + f01) {
if (output_detail) {
cout << "Not Monge! " << endl;
cout << "f(" << i << ", " << j << ") = " << f(i, j) << ", ";
cout << "f(" << i << ", " << j + 1 << ") = " << f(i, j + 1) << endl;
cout << "f(" << i + 1 << ", " << j << ") = " << f(i + 1, j) << ", ";
cout << "f(" << i + 1 << ", " << j + 1 << ") = " << f(i + 1, j + 1) << endl;
}
return false;
}
}
}
return true;
}
// check whether N x M matrix f is anti-Monge or not, in O(NM)
template<class FUNC> bool is_anti_monge
(int N, int M, const FUNC &f, bool upper_triangular = false, bool output_detail = true) {
assert(N > 0 && M > 0);
for (int i = 0; i + 1 < N; i++) {
for (int j = (upper_triangular ? i+2 : 0); j + 1 < M; j++) {
auto f00 = f(i, j), f01 = f(i, j + 1), f10 = f(i + 1, j), f11 = f(i + 1, j + 1);
if (f00 + f11 < f10 + f01) {
if (output_detail) {
cout << "Not anti-Monge! " << endl;
cout << "f(" << i << ", " << j << ") = " << f(i, j) << ", ";
cout << "f(" << i << ", " << j + 1 << ") = " << f(i, j + 1) << endl;
cout << "f(" << i + 1 << ", " << j << ") = " << f(i + 1, j) << ", ";
cout << "f(" << i + 1 << ", " << j + 1 << ") = " << f(i + 1, j + 1) << endl;
}
return false;
}
}
}
return true;
}
#line 4 "Optimization/monge_shortest_path.hpp"
using namespace std;
#line 4 "Optimization/simple_larsch.hpp"
using namespace std;
// noshi's simplified LARSCH
// find shortest path from vertex 0 on DAG in O(N log N)
// vertex: 0, 1, 2, ..., N
// f(i, j) must be Monge (Then, A(j, i) = dp[i] + f(i, j) is also Monge)
// slide access OK
template<class VAL> struct SimpleLARSCH {
VAL INF = numeric_limits<VAL>::max() / 2;
int CNT_INF = numeric_limits<int>::max() / 2;
// results
vector<VAL> dp;
vector<int> cnt, prev;
// solver (random access ver)
template<class FUNC> vector<pair<VAL, int>> solve(int N, const FUNC &f, bool minimize_cnt = true) {
dp.assign(N + 1, INF);
cnt.assign(N + 1, CNT_INF);
prev.assign(N + 1, 0);
dp[0] = 0, cnt[0] = 0;
auto relax = [&](int l, int r) -> void {
VAL val = dp[l] + f(l, r);
int cn = cnt[l] + 1;
if (dp[r] > val || (dp[r] == val && (minimize_cnt ? cn < cnt[r] : cn > cnt[r]))) {
dp[r] = val;
cnt[r] = cn;
prev[r] = l;
}
};
auto rec = [&](auto &&rec, int l, int r) -> void {
if (r - l <= 1) return;
int m = (l + r) / 2;
for (int k = prev[l]; k <= prev[r]; k++) relax(k, m);
rec(rec, l, m);
for (int k = l + 1; k <= m; k++) relax(k, r);
rec(rec, m, r);
};
if (N > 0) relax(0, N), rec(rec, 0, N);
vector<pair<VAL, int>> res(N + 1, make_pair(numeric_limits<VAL>::max() / 2, -1));
res[0].first = VAL(0);
for (int i = 1; i <= N; i++) res[i] = {dp[i], prev[i]};
return res;
}
// solver (slide access ver)
template<class STATE, class ADD, class DEL, class GETCOST> vector<pair<VAL, int>> solve
(int N, const STATE &ini, const ADD &add, const DEL &del, const GETCOST &get, bool minimize_cnt = true) {
dp.assign(N + 1, INF);
cnt.assign(N + 1, CNT_INF);
prev.assign(N + 1, 0);
dp[0] = 0, cnt[0] = 0;
STATE win[2] = { ini, ini };
int cl[2] = {0, 0}, cr[2] = {0, 0};
auto move_cursor = [&](int c, int l, int r) {
while (cr[c] < r) add(win[c], cr[c]++);
while (cl[c] > l) add(win[c], --cl[c]);
while (cr[c] > r) del(win[c], --cr[c]);
while (cl[c] < l) del(win[c], cl[c]++);
};
auto relax = [&](int c, int l, int r) {
move_cursor(c, l, r);
VAL val = dp[l] + get(win[c]);
int cn = cnt[l] + 1;
if (dp[r] > val || (dp[r] == val && (minimize_cnt ? cn < cnt[r] : cn > cnt[r]))) {
dp[r] = val;
cnt[r] = cn;
prev[r] = l;
}
};
auto rec = [&](auto &&rec, int l, int r) -> void {
if (r - l <= 1) return;
int m = (l + r) / 2;
for (int k = prev[l]; k <= prev[r]; k++) relax(0, k, m);
rec(rec, l, m);
for (int k = l + 1; k <= m; k++) relax(1, k, r);
rec(rec, m, r);
};
if (N > 0) relax(0, 0, N), rec(rec, 0, N);
vector<pair<VAL, int>> res(N + 1, make_pair(numeric_limits<VAL>::max() / 2, -1));
res[0].first = VAL(0);
for (int i = 1; i <= N; i++) res[i] = {dp[i], prev[i]};
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 4 "Optimization/larsch.hpp"
using namespace std;
// LARSCH, in O(N)
// vertex: 0, 1, 2, ..., N
// f(i, j) must be Monge (Then, A(j, i) = dp[i] + f(i, j) is also Monge)
template<class VAL, class FUNC> struct SubLARSCH {
struct RowReduce;
struct ColReduce;
struct ColMap {
const ColMap* parent = nullptr;
const vector<int> *v = nullptr;
constexpr int mapping(int i) const {
int x = (v ? (*v)[i] : i);
return (parent ? parent->mapping(x) : x);
}
};
struct Eval {
const FUNC *f = nullptr;
long long a = 1, b = 0; // row = a * i + b
const ColMap *cm = nullptr;
constexpr VAL operator () (int i, int j) const {
int i2 = int(a * i + b), j2 = (cm ? cm -> mapping(j) : j);
return (*f)(i2, j2);
}
};
struct RowReduce {
int N;
Eval e;
int cur_row = 0, state = 0;
unique_ptr<ColReduce> rec;
constexpr RowReduce(int N_, const Eval &e_) : N(N_), e(e_) {
int M = N / 2;
if (M > 0) {
Eval e2 = e;
e2.a = e.a * 2, e2.b = e.a + e.b;
rec = make_unique<ColReduce>(M, e2);
}
}
constexpr void reset() {
cur_row = 0, state = 0;
if (rec) rec->reset();
}
int get_argmin() {
int i = cur_row++;
if (!(i & 1)) {
int prev = state, next = (i + 1 == N ? N - 1 : rec->get_argmin());
state = next;
int res = prev;
for (int j = prev + 1; j <= next; j++) {
if (e(i, res) > e(i, j)) res = j;
}
return res;
} else {
return (e(i, state) <= e(i, i) ? state : i);
}
}
};
struct ColReduce {
int N;
Eval e;
int cur_row = 0;
vector<int> cols;
ColMap cm;
RowReduce rec;
constexpr ColReduce(int N_, const Eval &e_)
: N(N_), e(e_), cols(), cm{e.cm, &cols}, rec(N_, Eval{e.f, e.a, e.b, &cm}) {
cols.reserve(N);
}
constexpr void reset() {
cur_row = 0;
cols.clear();
rec.reset();
}
constexpr void push_col(int i, int j) {
while (!cols.empty()) {
int siz = (int)cols.size();
if (siz == i) break;
int last = cols.back();
if (e(siz - 1, last) > e(siz - 1, j)) cols.pop_back();
else break;
}
if ((int)cols.size() != N) cols.emplace_back(j);
}
constexpr int get_argmin() {
int i = cur_row++;
if (i == 0) {
cols.clear();
cols.emplace_back(0);
} else {
push_col(i, i * 2 - 1);
push_col(i, i * 2);
}
return cols[rec.get_argmin()];
}
};
FUNC f;
ColMap root_cm;
Eval root_eval;
unique_ptr<RowReduce> base;
explicit SubLARSCH(int N_, const FUNC &f_)
: f(std::move(f_)), root_cm{nullptr, nullptr}, root_eval{&f, 1, 0, &root_cm} {
base = make_unique<RowReduce>(N_, root_eval);
}
constexpr void reset() {
base->reset();
}
constexpr int get_argmin() {
return base->get_argmin();
}
};
template<class VAL> struct LARSCH {
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));
dp[0] = 0, cnt[0] = 0, res[0].first = 0;
auto f2 = [&](int i, int j) -> VAL {
i++;
if (i <= j) return INF;
else return dp[j] + f(j, i);
};
SubLARSCH<VAL, decltype(f2)> lar(N, f2);
for (int r = 1; r <= N; r++) {
int l = lar.get_argmin();
dp[r] = dp[l] + f(l, r);
prev[r] = l;
cnt[r] = cnt[l] + 1;
}
for (int r = 1; r <= N; r++) res[r] = {dp[r], prev[r]};
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 4 "Optimization/dd_monotone_minima.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 7 "Optimization/dd_monotone_minima.hpp"
// find shortest path on DAG by D&D Monotone Minima, O(N (log N)^2)
// vertex: 0, 1, 2, ..., N
// A(j, i) = dp[i] + f(i, j) must be monotone (argmin is not decreasing)
// Attention: Even if f is monotone, A is not necessarily monotone
// slide access OK
template<class VAL> struct DDMonotoneMinima {
VAL INF = numeric_limits<VAL>::max() / 2;
int CNT_INF = numeric_limits<int>::max() / 2;
// results
vector<VAL> dp;
vector<int> cnt, prev;
// solver (random access ver)
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);
dp[0] = 0, cnt[0] = 0;
vector<pair<VAL, int>> res(N + 1, make_pair(INF, -1));
vector<pair<VAL, int>> tmp(N + 1, make_pair(INF, -1));
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;
rec(rec, left, mid);
for (int i = mid; i < right; i++) tmp[i] = {INF, -1};
MonotoneMinimaRec(mid, right, left, mid, f2, tmp);
for (int i = mid; i < right; i++) if (tmp[i].first < res[i].first) res[i] = tmp[i];
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;
}
// solver (slide access ver)
template<class ADD, class DEL, class GETCOST> vector<pair<VAL, int>> solve
(int N, const ADD &add, const DEL &del, const GETCOST &get) {
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); };
return solve(N, f2);
}
// 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 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;
}
};
#line 10 "Optimization/monge_shortest_path.hpp"
// Monge Shortest Path Wrapper
// vertex: 0, 1, 2, ..., N
// f(i, j) must be Monge (Then, A(j, i) = dp[i] + f(i, j) is also Monge)
template<class VAL> struct MongeShortestPath {
VAL INF = numeric_limits<VAL>::max() / 2;
int CNT_INF = numeric_limits<int>::max() / 2;
// results
vector<VAL> dp;
vector<int> cnt, prev;
// solver (random access ver)
template<class FUNC> vector<pair<VAL, int>> solve(int N, const FUNC &f
, const string solver = "simple_larsch") {
vector<pair<VAL, int>> res(N + 1, make_pair(numeric_limits<VAL>::max() / 2, -1));
if (solver == "simple_larsch") {
SimpleLARSCH<VAL> opt;
res = opt.solve(N, f);
dp = opt.dp, cnt = opt.cnt, prev = opt.prev;
} else if (solver == "larsch") {
LARSCH<VAL> opt;
res = opt.solve(N, f);
dp = opt.dp, cnt = opt.cnt, prev = opt.prev;
} else if (solver == "monotone") {
DDMonotoneMinima<VAL> opt;
res = opt.solve(N, f);
dp = opt.dp, cnt = opt.cnt, prev = opt.prev;
} else if (solver == "smawk") {
DDSMAWK<VAL> opt;
res = opt.solve(N, f);
dp = opt.dp, cnt = opt.cnt, prev = opt.prev;
} else {
cout << "no solver" << endl;
assert(true);
}
return res;
}
// solver (slide access ver)
template<class STATE, class ADD, class DEL, class GETCOST> vector<pair<VAL, int>> solve
(int N, const STATE &ini, const ADD &add, const DEL &del, const GETCOST &get, const string solver = "simple_larsch") {
vector<pair<VAL, int>> res;
if (solver == "simple_larsch") {
SimpleLARSCH<VAL> opt;
res = opt.solve(N, ini, add, del, get);
dp = opt.dp, cnt = opt.cnt, prev = opt.prev;
} else if (solver == "monotone") {
DDMonotoneMinima<VAL> opt;
STATE st = ini;
int left = 0, right = 0;
auto f = [&](int i, int j) -> long long {
while (left < i) del(st, left++);
while (i < left) add(st, --left);
while (right < j) add(st, right++);
while (j < right) del(st, --right);
return get(st);
};
res = opt.solve(N, f);
dp = opt.dp, cnt = opt.cnt, prev = opt.prev;
}else {
cout << "no solver" << endl;
assert(true);
}
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 4 "Optimization/monge_shortest_path_by_cht.hpp"
using namespace std;
#line 4 "Optimization/general_li_chao_tree.hpp"
using namespace std;
// General Li-Chao Tree (Dynamic ver.) for Pairwise Single Crossing Functions
template<class VAL, class FUNC, class X = long long> struct DynamicGeneralLiChaoTree {
struct Node {
int id = -1; // function id
Node *left = nullptr, *right = nullptr;
};
// inner values using chunk-based memory pool
vector<unique_ptr<Node[]>> blocks;
vector<FUNC> funcs;
Node *root = nullptr;
X XL, XR;
VAL INF = numeric_limits<VAL>::max() / 2;
// constructors
explicit DynamicGeneralLiChaoTree(X xl, X xr, VAL inf = numeric_limits<VAL>::max() / 2)
: XL(xl), XR(xr), INF(inf) {
assert(xl < xr);
}
// basic functions
void reserve_functions(size_t n) {
funcs.reserve(n);
}
VAL eval(int id, X x) const {
return funcs[id](x);
}
bool better(int fid, int gid, X x) const {
VAL a = eval(fid, x), b = eval(gid, x);
return (a < b) || (a == b && fid < gid);
}
void clear() {
funcs.clear();
blocks.clear();
used = BLOCK;
root = nullptr;
}
// insert(f), in O(log X)
int insert(FUNC f) {
int id = (int)funcs.size();
funcs.push_back(std::move(f));
insert_at(root, XL, XR, id);
return id;
}
// insert_segment(l, r, f), in O(log^2 X)
int insert_segment(FUNC f, X l, X r) {
assert(XL <= l && l <= r && r <= XR);
int id = (int)funcs.size();
funcs.push_back(std::move(f));
if (l != r) insert_segment_at(root, XL, XR, l, r, id);
return id;
}
// query(x), in O(log X)
pair<VAL, int> query(X x) const {
assert(XL <= x && x < XR);
VAL best = INF;
int best_id = -1;
Node *p = root;
X left = XL, right = XR;
while (p) {
if (p->id != -1) {
int fid = p->id;
VAL y = eval(fid, x);
if (best_id == -1 || y < best || (y == best && fid < best_id)) {
best = y, best_id = fid;
}
}
if (right - left == 1) break;
X mid = (left + right) / 2;
if (x < mid) p = p->left, right = mid;
else p = p->right, left = mid;
}
return {best, best_id};
}
private:
// chunk-based node allocation
static constexpr int BLOCK = 4096;
int used = BLOCK;
Node *new_node() {
if (used == BLOCK) {
blocks.emplace_back(make_unique<Node[]>(BLOCK));
used = 0;
}
return &blocks.back()[used++];
}
// inner insert
void insert_at(Node *&p, X left, X right, int fid) {
Node **cur = &p;
while (true) {
if (!*cur) *cur = new_node();
Node *v = *cur;
int gid = v->id;
if (gid == -1) {
v->id = fid;
return;
}
bool left_better = better(fid, gid, left), right_better = better(fid, gid, right - 1);
if (!left_better && !right_better) return;
if (left_better && right_better) {
v->id = fid;
return;
}
X mid = (left + right) / 2;
bool mid_better = better(fid, gid, mid);
if (mid_better) v->id = fid, fid = gid;
if (left_better != mid_better) cur = &v->left, right = mid;
else cur = &v->right, left = mid;
}
}
// inner insert_segment
void insert_segment_at(Node *&p, X left, X right, X l, X r, int fid) {
if (r <= left || right <= l) return;
if (l <= left && right <= r) {
insert_at(p, left, right, fid);
return;
}
if (!p) p = new_node();
X mid = (left + right) / 2;
if (l < mid) insert_segment_at(p->left, left, mid, l, r, fid);
if (mid < r) insert_segment_at(p->right, mid, right, l, r, fid);
}
};
// General Li-Chao Tree (Coordinate Compression ver.) for Pairwise Single Crossing Functions
template<class VAL, class FUNC, class X = long long> struct GeneralLiChaoTree {
// inner values
int N, XL, SIZ;
bool compressed = true;
VAL INF = numeric_limits<VAL>::max() / 2;
vector<X> xs;
vector<int> ids;
vector<FUNC> funcs;
// constructors
explicit GeneralLiChaoTree(vector<X> xs_, VAL inf = numeric_limits<VAL>::max() / 2, bool sorted = true)
: N((int)xs_.size()), XL(0), SIZ(1), compressed(true), INF(inf), xs(std::move(xs_)) {
if (!sorted) {
sort(xs.begin(), xs.end());
xs.erase(unique(xs.begin(), xs.end()), xs.end());
}
while (SIZ < N) SIZ <<= 1;
ids.assign(SIZ << 1, -1);
}
explicit GeneralLiChaoTree(int xl, int xr, VAL inf = numeric_limits<VAL>::max() / 2)
: N(xr - xl), XL(xl), SIZ(1), compressed(false), INF(inf) {
assert(N >= 0);
while (SIZ < N) SIZ <<= 1;
ids.assign(SIZ << 1, -1);
}
// basic functions
void reserve_functions(size_t n) {
funcs.reserve(n);
}
VAL eval(int id, X x) const {
return funcs[id](x);
}
bool better(int fid, int gid, X x) const {
VAL a = eval(fid, x), b = eval(gid, x);
return (a < b) || (a == b && fid < gid);
}
void clear() {
funcs.clear();
fill(ids.begin(), ids.end(), -1);
}
// insert(f), in O(log N)
int insert(FUNC f) {
int id = (int)funcs.size();
funcs.push_back(std::move(f));
if (N) insert_at(1, 0, SIZ, id);
return id;
}
// insert_segment(f, l, r), in O(log^2 N)
int insert_segment(FUNC f, X l, X r) {
assert(l <= r);
int id = (int)funcs.size();
funcs.push_back(std::move(f));
if (N == 0 || l == r) return id;
int L = lower_bound_index(l);
int R = lower_bound_index(r);
if (L == R) return id;
int a = L + SIZ, b = R + SIZ;
while (a < b) {
if (a & 1) {
int v = a++;
int depth = 31 - __builtin_clz((unsigned)v);
int length = SIZ >> depth;
int left = (v - (1 << depth)) * length;
insert_at(v, left, left + length, id);
}
if (b & 1) {
int v = --b;
int depth = 31 - __builtin_clz((unsigned)v);
int length = SIZ >> depth;
int left = (v - (1 << depth)) * length;
insert_at(v, left, left + length, id);
}
a >>= 1, b >>= 1;
}
return id;
}
// query(x), x must be in xs, in O(log N)
pair<VAL, int> query(X x) const {
int pos = 0;
if (compressed) {
auto it = lower_bound(xs.begin(), xs.end(), x);
assert(it != xs.end() && *it == x);
pos = int(it - xs.begin());
} else {
assert(X(XL) <= x && x < X(XL + N));
pos = int(x - XL);
}
VAL best = INF;
int best_id = -1;
int v = pos + SIZ;
while (v) {
int fid = ids[v];
if (fid != -1) {
VAL y = eval(fid, x);
if (best_id == -1 || y < best || (y == best && fid < best_id)) {
best = y, best_id = fid;
}
}
v >>= 1;
}
return {best, best_id};
}
private:
// coordinate compression
int lower_bound_index(X x) const {
if (compressed) {
return int(std::lower_bound(xs.begin(), xs.end(), x) - xs.begin());
} else {
if (x <= X(XL)) return 0;
else if (x >= X(XL) + X(N)) return N;
else return int(x - X(XL));
}
}
// inner insert
void insert_at(int v, int left, int right, int fid) {
auto coord = [&](int i) -> X {
i = min(i, N - 1);
return (compressed ? xs[i] : X(i) + X(XL));
};
while (true) {
int gid = ids[v];
if (gid == -1) {
ids[v] = fid;
return;
}
bool left_better = better(fid, gid, coord(left));
bool right_better = better(fid, gid, coord(right - 1));
if (!left_better && !right_better) return;
if (left_better && right_better) {
ids[v] = fid;
return;
}
X mid = (left + right) / 2;
bool mid_better = better(fid, gid, coord(mid));
if (mid_better) ids[v] = fid, fid = gid;
if (left_better != mid_better) {
v <<= 1;
right = mid;
} else {
v = v << 1 | 1;
left = mid;
}
}
}
};
#line 7 "Optimization/monge_shortest_path_by_cht.hpp"
// find shortest path on DAG with Pairwise Single Crossing cost, by Li Chao Tree, O(N log N)
// vertex: 0, 1, 2, ..., N
// A(j, i) = dp[i] + f(i, j) must be Pairwise Single Crossing (must be in all 0 <= i, j <= N, including i >= j)
// Attention: Even if f is Pairwise Single Crossing, A is not necessarily Pairwise Single Crossing
// edges: (i, j), j >= start[i] (i + 1 <= start[i] <= N + 1), start is non-decreasing
template<class VAL, class FUNC> vector<pair<VAL, int>> MongeShortestPathByCHT
(int N, const FUNC &f, const vector<int> &start, bool use_segment_insertion = true) {
assert((int)start.size() == N);
for (int i = 0; i < N; i++) {
assert(i + 1 <= start[i] && start[i] <= N + 1);
if (i > 0) assert(start[i - 1] <= start[i]);
}
VAL INF = numeric_limits<VAL>::max() / 2;
vector<pair<VAL, int>> dp(N + 1, make_pair(numeric_limits<VAL>::max() / 2, -1));
auto make_func = [&](int i) {
return [&, i](long long j) -> VAL { return dp[i].first + f(i, j); };
};
GeneralLiChaoTree<VAL, decltype(make_func(0))> cht(1, N + 1, INF);
dp[0].first = 0;
int i = 0;
for (int j = 1; j <= N; j++) {
while (i < N && j >= start[i]) {
if (use_segment_insertion) {
// 挿入する関数群が j >= start[i] の領域でのみ PSC であるとき、use_segment_insertion = true にする
cht.insert_segment(make_func(i), start[i], N + 1);
} else {
// 挿入する関数群が全ての i, j に対して PSC であるとき、use_segment_insertion = false にする
cht.insert(make_func(i));
}
i++;
}
dp[j] = cht.query(j);
};
return dp;
}
template<class VAL, class FUNC> vector<pair<VAL, int>> MongeShortestPathByCHT
(int N, const FUNC &f, bool use_segment_insertion = true) {
vector<int> start(N);
for (int i = 0; i < N; i++) start[i] = i + 1;
return MongeShortestPathByCHT<VAL>(N, f, start, use_segment_insertion);
}
#line 4 "Optimization/monge_enum_d_shortest_paths.hpp"
using namespace std;
#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;
}
};
#line 4 "Optimization/alien_dp.hpp"
using namespace std;
#line 4 "Optimization/alien_trick.hpp"
using namespace std;
// Alien Trick (by ラグランジュ緩和)
// lag(λ): x✳︎ = argmin_{x} {f(x) + λ g(x)} としたときのペア値 {f(x✳︎), g(x✳︎)} を返す関数 (最大化の方では argmax)
// 返り値: {最適値, そのときの λ}
template<class VAL> struct AlienTrick {
// solver
template<class FUNC> pair<VAL, VAL> solve(const FUNC &lag, VAL K) {
// zero check
auto [f0, g0] = lag(VAL(0));
if (g0 == K) return {f0, VAL(0)};
// find an initial value of binary search (by doubling)
VAL lo, hi, flo, fhi, glo, ghi;
if (g0 > K) {
lo = 0, flo = f0, glo = g0;
tie(hi, fhi, ghi) = find_upper_lambda(lag, K);
} else {
hi = 0, fhi = f0, ghi = g0;
tie(lo, flo, glo) = find_lower_lambda(lag, K);
}
if (glo == K) return {flo, lo};
if (ghi == K) return {fhi, hi};
assert(glo > K && ghi < K);
// binary search
// invariant: lo < hi, g(lo) > K, g(hi) < K
while (hi - lo > VAL(1)) {
VAL mid = (lo + hi) / 2;
auto [fmid, gmid] = lag(mid);
if (gmid == K) return {fmid, mid};
if (gmid > K) lo = mid, flo = fmid, glo = gmid;
else hi = mid, fhi = fmid, ghi = gmid;
}
// When the minimizer skips g(x) = K due to ties, recover the optimum from the Lagrangian dual value
VAL reslo = flo + lo * (glo - K), reshi = fhi + hi * (ghi - K);
if (reslo >= reshi) return {reslo, lo};
else return {reshi, hi};
}
template<class FUNC> pair<VAL, VAL> minimize(const FUNC &lag, VAL K) {
return solve(lag, K);
}
template<class FUNC> pair<VAL, VAL> maximize(const FUNC &lag, VAL K) {
auto reverse_lag = [&](VAL lambda) -> pair<VAL, VAL> {
auto [f, g] = lag(-lambda);
return {-f, g};
};
auto [f, g] = solve(reverse_lag, K);
return {-f, -g};
}
// in case: L <= g(x) <= R
template<class FUNC> pair<VAL, VAL> solve(const FUNC &lag, VAL L, VAL R) {
assert(L <= R);
auto [f0, g0] = lag(VAL(0));
if (L <= g0 && g0 <= R) return {f0, VAL(0)};
if (g0 > R) return solve(lag, R);
else return solve(lag, L);
}
template<class FUNC> pair<VAL, VAL> minimize(const FUNC &lag, VAL L, VAL R) {
return solve(lag, L, R);
}
template<class FUNC> pair<VAL, VAL> maximize(const FUNC &lag, VAL L, VAL R) {
auto reverse_lag = [&](VAL lambda) -> pair<VAL, VAL> {
auto [f, g] = lag(-lambda);
return {-f, g};
};
auto [f, g] = solve(reverse_lag, L, R);
return {-f, -g};
}
// find lo such that g(lo) >= R
template<class FUNC> tuple<VAL, VAL, VAL> find_lower_lambda(const FUNC &lag, VAL K) {
auto [f0, g0] = lag(VAL(0));
VAL lo = 0, hi = 0, flo = f0, glo = g0;
if (glo >= K) return tie(lo, flo, glo);
lo = VAL(-1);
tie(flo, glo) = lag(lo);
while (glo < K) {
assert(lo < hi);
VAL dif = hi - lo;
hi = lo, lo -= VAL(2) * dif;
tie(flo, glo) = lag(lo);
}
return tie(lo, flo, glo);
}
// find hi such that g(hi) <= L
template<class FUNC> tuple<VAL, VAL, VAL> find_upper_lambda(const FUNC &lag, VAL K) {
auto [f0, g0] = lag(VAL(0));
VAL lo = 0, hi = 0, fhi = f0, ghi = g0;
if (ghi <= K) return tie(hi, fhi, ghi);
hi = VAL(1);
tie(fhi, ghi) = lag(hi);
while (ghi > K) {
assert(lo < hi);
VAL dif = hi - lo;
lo = hi, hi += VAL(2) * dif;
tie(fhi, ghi) = lag(hi);
}
return tie(hi, fhi, ghi);
}
};
#line 8 "Optimization/alien_dp.hpp"
// Alien DP
// find the d-edges shortest path from vertex 0 on DAG by Lagrange relaxation
// vertex: 0, 1, 2, ..., N, f(i, j) must be Monge
template<class VAL> struct AlienDP {
// inner values
MongeShortestPath<VAL> msp;
AlienTrick<VAL> at;
// solver (random access ver.)
template<class FUNC> pair<VAL, VAL> solve(int N, const FUNC &f, int D, const string &solver = "simple_larsch") {
auto lag = [&](VAL lambda) -> pair<VAL, VAL> {
auto fg = [&](int i, int j) -> VAL { return f(i, j) + lambda; };
vector<pair<VAL, int>> obj = msp.solve(N, fg, solver);
auto cnt = msp.cnt[N];
auto res = obj.back().first - lambda * cnt;
return {res, cnt};
};
return at.solve(lag, D);
}
// solver (slide access ver.)
template<class STATE, class ADD, class DEL, class GETCOST> pair<VAL, int> solve
(int N, const STATE &ini, const ADD &add, const DEL &del, const GETCOST &get, int D, const string &solver = "simple_larsch") {
auto lag = [&](VAL lambda) -> pair<VAL, VAL> {
auto get_lag = [&](STATE &st) -> VAL { return get(st) + lambda; };
vector<pair<VAL, int>> obj = msp.solve(N, ini, add, del, get_lag, solver);
auto cnt = msp.cnt[N];
auto res = obj.back().first - lambda * cnt;
return {res, cnt};
};
return at.solve(lag, D);
}
// debugger
template<class FUNC> void debug(int N, const FUNC &f, VAL start = 0, VAL goal = 10) {
auto lag = [&](VAL lambda) -> void {
auto fg = [&](int i, int j) -> VAL { return f(i, j) + lambda; };
const auto &dp = msp.solve(N, fg);
auto cnt = msp.cnt[N];
auto res = dp.back().first - lambda * cnt;
cout << lambda << ": " << res << "(" << cnt << "); (";
for (int i = 0; i < (int)dp.size(); i++) {
if (i) cout << " ";
cout << "<" << dp[i].first << ", " << dp[i].second << ">";
}
cout << ")" << endl;
};
for (VAL lambda = start; lambda <= goal; lambda++) lag(lambda);
}
};
#line 4 "Optimization/anti_monge_dd_smawk.hpp"
using namespace std;
#line 7 "Optimization/anti_monge_dd_smawk.hpp"
// find shortest path on DAG with anti-Monge cost, by D&D SMAWK, O(NK log N)
// vertex: 0, 1, 2, ..., N
// f(i, j) must be anti-Monge in j >= i + K
// Then, A(j, i) = dp[i] + f(i, j) is also anti-Monge
// slide access NG
template<class VAL, class FUNC> vector<pair<VAL, int>> antiMongeDDSMAWK(int N, const FUNC &f, int K = 1) {
assert(K >= 1);
const VAL INF = numeric_limits<VAL>::max() / 4;
vector<pair<VAL, int>> res(N + 1, make_pair(INF, -1));
vector<pair<VAL, int>> tmp(N + 1, {INF, -1});
res[0].first = VAL(0);
auto relax = [&](int i, int j) -> void {
if (i != 0 && res[i].second == -1) return;
VAL val = res[i].first + f(i, j);
if (val < res[j].first) res[j] = {val, i};
};
auto rec = [&](auto &&rec, int left, int right) -> void {
if (right - left <= 1) return;
int mid = (left + right) / 2, border = min(right, mid + K - 1);
rec(rec, left, mid);
for (int j = mid; j < border; j++) {
for (int i = left; i < mid && i <= j - K; i++) relax(i, j);
}
if (border < right) {
vector<int> X, Y;
for (int j = border; j < right; j++) X.push_back(j);
for (int i = mid - 1; i >= left; i--) {
if (i == 0 || res[i].second != -1) Y.push_back(i);
}
if (!Y.empty()) {
vector<int> cols(Y.size());
iota(cols.begin(), cols.end(), 0);
auto f2 = [&](int j, int k) -> VAL { return res[Y[k]].first + f(Y[k], j); };
SMAWKRec(X, cols, f2, tmp);
for (auto x : X) {
if (tmp[x].first < res[x].first) res[x] = {tmp[x].first, Y[tmp[x].second]};
}
}
}
rec(rec, mid, right);
};
rec(rec, 0, N + 1);
return res;
}
#line 4 "Optimization/knuth_yao.hpp"
using namespace std;
// f: Monge, interval monotone ((N+1) x (N+1))
// dp[i][j] = min_{i <= k < j} (dp[i][k] + dp[k][j]) + f(i, j), in O(N^2)
template<class VAL, class FUNC> VAL KnuthYao(int N, const FUNC &f) {
VAL INF = numeric_limits<VAL>::max() / 2;
vector<vector<VAL>> dp(N+1, vector<VAL>(N+1, INF));
vector<vector<int>> K(N+1, vector<int>(N+1, -1));
for (int i = 0; i <= N; i++) {
dp[i][i] = 0, K[i][i] = i;
if (i < N) dp[i][i+1] = f(i, i+1), K[i][i+1] = i;
}
for (int between = 2; between <= N; between++) {
for (int i = 0; i + between <= N; i++) {
int j = i + between;
for (int k = K[i][j-1]; k <= K[i+1][j]; k++) {
VAL tmp = dp[i][k] + dp[k][j] + f(i, j);
if (dp[i][j] > tmp) {
dp[i][j] = tmp;
K[i][j] = k;
}
}
}
}
return dp[0][N];
}
#line 4 "Optimization/monge_min_plus_composition.hpp"
using namespace std;
// C[i][j] = min_{i <= k < j} (A[i][k] + B[k][j])
// A, B: 上三角 (N+1) x (N+1) Monge
template<class VAL, class FUNCA, class FUNCB> vector<vector<VAL>> min_plus_composition
(int N, const FUNCA &A, const FUNCB &B) {
VAL INF = numeric_limits<VAL>::max() / 2;
vector<vector<VAL>> res(N+1, vector<VAL>(N+1, INF));
vector<int> argmin(N+1);
for (int i = 0; i <= N; i++) res[i][i] = A(i, i) + B(i, i), argmin[i] = i;
for (int s = 1; s <= N; s++) {
vector<int> argmin2(N + 1 - s);
for (int i = 0; i <= N - s; i++) {
int j = i + s, p = argmin[i], q = argmin[i + 1];
for (int k = p; k <= q; k++) {
if (res[i][j] > A(i, k) + B(k, j)) {
res[i][j] = A(i, k) + B(k, j);
argmin2[i] = k;
}
}
}
swap(argmin, argmin2);
}
return res;
}
#line 15 "Optimization/monge_all.hpp"