Monge グラフの d 辺最短路

graph/monge_d_edge_shortest_path.hpp · GitHub で見る

共通ページ

頂点番号が増える向きに全ての辺を持つ Monge DAG で、始点から終点までちょうど指定本数の辺を使う最短路長を求めます。負の辺コストも扱えます。辺は列挙せず、コストを返す関数で渡します。

#include "lib/graph/monge_d_edge_shortest_path.hpp"
#include <cassert>

auto weight = [](int i, int j) { return 1LL * (j - i) * (j - i); };
auto answer = monge_d_edge_shortest_path(6, 2, weight);
assert(answer && *answer == 13);
assert(!monge_d_edge_shortest_path(6, 0, weight));
assert(monge_d_edge_shortest_path(1, 0, weight) == 0);

monge_d_edge_shortest_path<Distance = std::int64_t>(n, d, weight) は std::optional<Distance> を返します。n は頂点数、d は使用する辺数で、どちらも非負の int です。頂点は [0,n)、始点は 0、終点は n-1 です。経路長だけを返し、経路の復元は行いません。

weight(i,j) は 0 <= i < j < n の全てについて有限な辺コストを返す必要があります。自己辺・逆向きの辺・欠けた辺は扱いません。コールバックは符号付き 64 bit 以下の整数型を返し、同じ引数には同じ値を返してください。

全ての添字 \(0\le a<b<c<e<n\) で、次の Monge 不等式が必要です。条件の検査は行いません。

\[w(a,c)+w(b,e)\le w(a,e)+w(b,c).\]

辺数と空入力の扱いは次のとおりです。

入力	返り値
n == 0	std::nullopt
n == 1 && d == 0	0
n == 1 && d > 0	std::nullopt
n > 1 && (d == 0 || d >= n)	std::nullopt
n > 1 && 1 <= d && d < n	ちょうど d 辺を使う最短路長

上の表の最後の行以外では weight を呼びません。辺コストや複数の経路が同じ値でも使用できます。

Distance は符号付き 128 bit 以下の整数型です。既定の 64 bit に収まらない答えが必要な場合は、monge_d_edge_shortest_path<__int128_t>(n,d,weight) を使えます。内部の距離・ペナルティ・比較は __int128_t で計算し、返り値が Distance の範囲外なら std::overflow_error を送出します。コールバック内の演算はこの拡張の対象外なので、weight 自体がオーバーフローしない型で計算してください。

ちょうど \(k\) 辺の最短路長を \(F(k)\) とすると、Monge 性により \(F\) は離散凸です。各辺に整数ペナルティ \(\lambda\) を加えた最短路を計算し、同じ長さでは辺数が少ない経路を選んで、必要な傾きを二分探索します。最適辺数が飛ぶ場合も、同じ傾き上にある \(F(d)\) を求められます。

探索区間は次の値から自動で計算するので、コストの上限やペナルティの範囲を指定する必要はありません。

\[F(n-2)-F(n-1)\le\lambda\le F(1)-F(2).\]

ここで \(F(1)=w(0,n-1)\)、\(F(2)=\displaystyle\min_{1\le i<n-1}(w(0,i)+w(i,n-1))\) です。\(F(n-1)\) は隣接する全ての辺の和、\(F(n-2)\) はその経路から内部頂点を 1 個除いたときの最小値で、それぞれ線形時間で求められます。ペナルティを固定した最短路計算には Online Monge Minima を使います。

参考: Nyaan’s Library、maspypy library。

weight の評価を \(O(1)\)、全ての辺について \(\lvert w(i,j)\rvert\le W\)、\(W\ge2\) とすると、時間は \(O(n\log n\log W)\)、追加空間は \(O(n)\) です。ペナルティの探索幅を \(R\) とすると、\(R\ge2\) で時間は \(O(n\log n\log R)\)、幅が 0 または 1 の場合は \(O(n\log n)\) です。d == 1 は \(O(1)\) 時間・空間、d == n-1・d == 2・d == n-2 は \(O(n)\) 時間・空間で直接計算します。コールバックを呼ばない入力は \(O(1)\) 時間・空間です。

依存ライブラリ

• Online Monge Minima

コード

#ifndef MONGE_D_EDGE_SHORTEST_PATH_HPP
#define MONGE_D_EDGE_SHORTEST_PATH_HPP

#include "../components/online_monge_minima.hpp"
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <functional>
#include <limits>
#include <optional>
#include <stdexcept>
#include <type_traits>
#include <vector>

template <class Distance = std::int64_t, class Weight>
std::optional<Distance> monge_d_edge_shortest_path(int n, int d, Weight weight) {
	using T = std::remove_cvref_t<std::invoke_result_t<Weight &, int, int>>;
	static_assert(std::numeric_limits<T>::is_integer && std::numeric_limits<T>::is_signed && sizeof(T) <= 8);
	static_assert(std::numeric_limits<Distance>::is_integer && std::numeric_limits<Distance>::is_signed && sizeof(Distance) <= 16);
	using Wide = __int128_t;
	assert(n >= 0 && d >= 0);
	if (!n || d >= n || (!d && n > 1)) return std::nullopt;
	if (n == 1) return Distance(0);
	auto result = [](Wide value) -> std::optional<Distance> {
		if (value < Wide(std::numeric_limits<Distance>::lowest()) || value > Wide(std::numeric_limits<Distance>::max()))
			throw std::overflow_error("Monge shortest path distance");
		return Distance(value);
	};
	Wide single = weight(0, n - 1);
	if (d == 1) return result(single);
	std::vector<Wide> adjacent(n - 1);
	Wide full = 0;
	for (int i = 0; i < n - 1; ++i) full += adjacent[i] = weight(i, i + 1);
	if (d == n - 1) return result(full);
	Wide two = Wide(weight(0, 1)) + Wide(weight(1, n - 1));
	Wide removed = Wide(weight(0, 2)) - adjacent[0] - adjacent[1];
	for (int i = 2; i < n - 1; ++i) {
		two = std::min(two, Wide(weight(0, i)) + Wide(weight(i, n - 1)));
		removed = std::min(removed, Wide(weight(i - 1, i + 1)) - adjacent[i - 1] - adjacent[i]);
	}
	if (d == 2) return result(two);
	if (d == n - 2) return result(full + removed);
	Wide low = removed, high = single - two;
	assert(low <= high);
	std::vector<Wide> distance(n);
	std::vector<int> edges(n);
	auto solve = [&](Wide penalty) {
		online_monge_minima(n, [&](int row, int old_col, int new_col) {
			Wide a = distance[old_col] + Wide(weight(old_col, row));
			Wide b = distance[new_col] + Wide(weight(new_col, row));
			return b < a || (b == a && edges[new_col] < edges[old_col]);
		}, [&](int row, int column) {
			distance[row] = distance[column] + Wide(weight(column, row)) + penalty;
			edges[row] = edges[column] + 1;
		});
	};
	Wide best = single + high;
	while (low < high) {
		Wide mid = low + (high - low) / 2;
		solve(mid);
		if (edges[n - 1] == d) return result(distance[n - 1] - mid * d);
		if (edges[n - 1] < d) high = mid, best = distance[n - 1];
		else low = mid + 1;
	}
	return result(best - high * d);
}

#endif

Back to top page