Merge pull request #348 from hitonanode/monge-shortest-path

Monge shortest path

other_algorithms/monge_shortest_path.hpp

@@ -0,0 +1,84 @@
+#pragma once
+#include <cassert>
+#include <vector>
+
+// Shortest path of Monge-weighted graph
+// Variant of LARSCH Algorithm: https://noshi91.hatenablog.com/entry/2023/02/18/005856
+// Complexity: O(n log n)
+//
+// Given a directed graph with n vertices and weighted edges
+// (w(i, j) = cost_callback(i, j) (i < j)),
+// this class calculates the shortest path from vertex 0 to all other vertices.
+template <class Cost> struct monge_shortest_path {
+    std::vector<Cost> dist; // dist[i] = shortest distance from 0 to i
+    std::vector<int> amin;  // amin[i] = previous vertex of i in the shortest path
+
+    template <class F> void _check(int i, int k, F cost_callback) {
+        if (i <= k) return;
+        if (Cost c = dist[k] + cost_callback(k, i); c < dist[i]) dist[i] = c, amin[i] = k;
+    }
+
+    template <class F> void _rec_solve(int l, int r, F cost_callback) {
+        if (r - l == 1) return;
+
+        const int m = (l + r) / 2;
+        for (int k = amin[l]; k <= amin[r]; ++k) _check(m, k, cost_callback);
+
+        _rec_solve(l, m, cost_callback);
+        for (int k = l + 1; k <= m; ++k) _check(r, k, cost_callback);
+        _rec_solve(m, r, cost_callback);
+    }
+
+    template <class F> Cost solve(int n, F cost_callback) {
+        assert(n > 0);
+        dist.resize(n);
+        amin.assign(n, 0);
+
+        dist[0] = Cost();
+        for (int i = 1; i < n; ++i) dist[i] = cost_callback(0, i);
+
+        _rec_solve(0, n - 1, cost_callback);
+
+        return dist.back();
+    }
+
+    template <class F> int num_edges() const {
+        int ret = 0;
+        for (int c = (int)amin.size() - 1; c >= 0; c = amin[c]) ++ret;
+        return ret;
+    }
+};
+
+// Find shortest path length from 0 to n - 1 with k edges, min_edges <= k <= max_edges
+// https://noshi91.hatenablog.com/entry/2022/01/13/001217
+template <class Cost, class F>
+Cost monge_shortest_path_with_specified_edges(int n, int min_edges, int max_edges,
+                                              Cost max_abs_cost, F cost_callback) {
+
+    assert(1 <= n);
+    assert(0 <= min_edges);
+    assert(min_edges <= max_edges);
+    assert(max_edges <= n - 1);
+
+    monge_shortest_path<Cost> msp;
+
+    auto eval = [&](Cost p) -> Cost {
+        msp.solve(n, [&](int i, int j) { return cost_callback(i, j) - p; });
+        return -msp.dist.back() - p * (p < 0 ? max_edges : min_edges);
+    };
+
+    Cost lo = -max_abs_cost * 3, hi = max_abs_cost * 3;
+
+    while (lo + 1 < hi) {
+        Cost p = (lo + hi) / 2, f = eval(p), df = eval(p + 1) - f;
+        if (df == Cost()) {
+            return -f;
+        } else {
+            (df < Cost() ? lo : hi) = p;
+        }
+    }
+
+    Cost flo = eval(lo), fhi = eval(hi);
+
+    return flo < fhi ? -flo : -fhi;
+}

other_algorithms/monge_shortest_path.md

@@ -0,0 +1,47 @@
+---
+title: Shortest path of DAG with Monge weights
+documentation_of: ./monge_shortest_path.hpp
+---
+
+$n$ 頂点の DAG で辺重みが Monge となるようなものに対して最短路長を高速に求める． [1] で紹介されている簡易版 LARSCH Algorithm が実装されていて，計算量は $O(n \log n)$ ．
+
+また，辺数が `min_edges` 以上 `max_edges` 以下であるようなものに限った最短路長を高速に求めることも可能（計算量にさらに重み二分探索の $\log$ がつく）．
+
+## 使用方法
+
+### 最短路長の計算
+
+```cpp
+auto f = [&](int s, int t) -> Cost {
+    //
+};
+
+monge_shortest_path<Cost> msp;
+Cost ret = msp.solve(n, f);
+```
+
+### 辺の本数の下限・上限を指定した最短路長の計算
+
+```cpp
+auto f = [&](int s, int t) -> Cost {
+    //
+};
+
+int n;  // 頂点数
+int l, r;  // 辺の本数が [l, r] の範囲に収まる最短路を見つけたい
+Cost max_weight;  // f() が返す値の絶対値の上界
+
+Cost ret = monge_shortest_path_with_specified_edges(n, l, r, max_weight, f);
+```
+
+## 問題例
+
+- [No.705 ゴミ拾い Hard - yukicoder](https://yukicoder.me/problems/no/705)
+- [AtCoder Beginner Contest 218 H - Red and Blue Lamps](https://atcoder.jp/contests/abc218/tasks/abc218_h)
+- [東京海上日動プログラミングコンテスト2024（AtCoder Beginner Contest 355） G - Baseball](https://atcoder.jp/contests/abc355/tasks/abc355_g)
+- [東北大学プログラミングコンテスト 2022 K - Lebesgue Integral](https://atcoder.jp/contests/tupc2022/tasks/tupc2022_k)
+
+## Links
+
+- [1] [簡易版 LARSCH Algorithm - noshi91のメモ](https://noshi91.hatenablog.com/entry/2023/02/18/005856)
+- [2] [Aliens DP における二分探索の色々 - noshi91のメモ](https://noshi91.hatenablog.com/entry/2023/11/20/052227#fn-c9578a2a)

other_algorithms/test/concave_max_plus_convolution.test.cpp

@@ -0,0 +1,22 @@
+#define PROBLEM "https://judge.yosupo.jp/problem/min_plus_convolution_convex_arbitrary"
+
+#include "../smawk.hpp"
+
+#include <iostream>
+#include <vector>
+using namespace std;
+
+int main() {
+    cin.tie(nullptr);
+    ios::sync_with_stdio(false);
+
+    int N, M;
+    cin >> N >> M;
+    vector<int> A(N), B(M);
+    for (auto &a : A) cin >> a, a = -a;
+    for (auto &b : B) cin >> b, b = -b;
+
+    auto ret = concave_max_plus_convolution<int, ((1 << 30) - 1) * 2>(B, A);
+
+    for (int i = 0; i < N + M - 1; ++i) cout << -ret[i] << " \n"[i + 1 == N + M - 1];
+}

other_algorithms/test/monge_shortest_path.yuki705.test.cpp

@@ -0,0 +1,28 @@
+#define PROBLEM "https://yukicoder.me/problems/no/705"
+#include "../monge_shortest_path.hpp"
+
+#include <cassert>
+#include <cmath>
+#include <iostream>
+using namespace std;
+
+int main() {
+    cin.tie(nullptr), ios::sync_with_stdio(false);
+
+    int N;
+    cin >> N;
+    vector<int> A(N), X(N), Y(N);
+    for (auto &a : A) cin >> a;
+    for (auto &x : X) cin >> x;
+    for (auto &y : Y) cin >> y;
+
+    auto weight = [&](int j, int i) {
+        assert(j < i);
+        --i;
+        const long long dx = abs(A.at(i) - X.at(j)), dy = Y.at(j);
+        return dx * dx * dx + dy * dy * dy;
+    };
+
+    monge_shortest_path<long long> msp;
+    cout << msp.solve(N + 1, weight) << '\n';
+}