algorithm
Implementation of various algorithms
tests/yukicoder/yukicoder_0705_monge_shortest_path.test.cpp
- View this file on GitHub
- Last update: 2026-11-12 11:56:29+09:00
- Problem: https://yukicoder.me/problems/no/705
Depends on
-
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)
-
Monge 単一始点最短路問題 (by LARSCH 法, in O(N))
(Optimization/larsch.hpp)
-
Monge 単一始点最短路アルゴリズム群のラッパー
(Optimization/monge_shortest_path.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
#define PROBLEM "https://yukicoder.me/problems/no/705"
#include "monge_shortest_path.hpp"
// yukicoder No.705 ゴミ拾い Hard
/*
A, X, Y: N 個
これらを区間に分割していく
dp[j] = min_{0 ≦ i < j}(dp[i] + |A[j-1] - X[i]|^3 + |-Y[i]|^3)
i -> j のコスト：|A[j-1] - X[i]|^3 + |-Y[i]|^3 ...... 差の凸関数 (Monge) + 縞々 (Monge) -> Monge
スタート: 0, ゴール: N
*/
int main() {
int N;
cin >> N;
vector<long long> A(N), X(N), Y(N);
for (int i = 0; i < N; i++) cin >> A[i];
for (int i = 0; i < N; i++) cin >> X[i];
for (int i = 0; i < N; i++) cin >> Y[i];
auto func = [&](int i, int j) -> long long {
long long dx = abs(A[j-1] - X[i]), dy = abs(Y[i]);
return dx * dx * dx + dy * dy * dy;
};
MongeShortestPath<long long> opt;
auto res = opt.solve(N, func, "simple_larsch");
cout << res[N].first << endl;
}#line 1 "tests/yukicoder/yukicoder_0705_monge_shortest_path.test.cpp"
#define PROBLEM "https://yukicoder.me/problems/no/705"
#line 2 "Optimization/monge_shortest_path.hpp"
#include <bits/stdc++.h>
#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 "tests/yukicoder/yukicoder_0705_monge_shortest_path.test.cpp"
// yukicoder No.705 ゴミ拾い Hard
/*
A, X, Y: N 個
これらを区間に分割していく
dp[j] = min_{0 ≦ i < j}(dp[i] + |A[j-1] - X[i]|^3 + |-Y[i]|^3)
i -> j のコスト：|A[j-1] - X[i]|^3 + |-Y[i]|^3 ...... 差の凸関数 (Monge) + 縞々 (Monge) -> Monge
スタート: 0, ゴール: N
*/
int main() {
int N;
cin >> N;
vector<long long> A(N), X(N), Y(N);
for (int i = 0; i < N; i++) cin >> A[i];
for (int i = 0; i < N; i++) cin >> X[i];
for (int i = 0; i < N; i++) cin >> Y[i];
auto func = [&](int i, int j) -> long long {
long long dx = abs(A[j-1] - X[i]), dy = abs(Y[i]);
return dx * dx * dx + dy * dy * dy;
};
MongeShortestPath<long long> opt;
auto res = opt.solve(N, func, "simple_larsch");
cout << res[N].first << endl;
}