牛ゲー（差分制約系）

graph/difference_constraints.hpp · GitHub で見る

共通ページ

整数変数に対する差分制約を扱います。解の存在判定・復元と、変数間の差や各変数の最大値・最小値を求められます。

#include "lib/graph/difference_constraints.hpp"

DifferenceConstraints<> system(3);
system.add_range(0, 1, 2, 5); // 2 <= x[1] - x[0] <= 5
system.add_lt(1, 2, 4);       // x[2] - x[1] < 4
system.fix(0, 10);            // x[0] = 10

auto result = system.solve();
if (result.feasible()) {
	const auto &x = result.values(); // 全制約を満たす解の一つ
	assert(x[0] == 10);
	assert(result.max_difference(0, 2) == 8);
	assert(result.min_difference(0, 1) == 2);
	assert(result.upper_bound(2) == 18);
	assert(!result.lower_bound(2)); // 下に非有界
	auto upper = result.upper_bounds(0); // x[v] - x[0] の最大値: 0, 5, 8
}

DifferenceConstraints<Cost = std::int64_t>(n = 0) で変数の個数を指定します。変数番号は [0,n)、origin() は定数 0 を表す頂点 n です。二変数の制約・差の問い合わせでは origin() も指定できます。差の向きは全て x[v] - x[u] です。

既定の計算・戻り値は std::int64_t です。128 bit が必要な場合は DifferenceConstraints<__int128_t> を使ってください。Cost は 16 byte 以下の符号付き整数型で、戻り値の value_type は Cost が 8 byte 以下なら std::int64_t、それ以外なら __int128_t です。値と途中計算が型の範囲に収まるようにしてください。扱えない値を検出すると std::overflow_error を投げます。

制約追加	意味
add_le(u,v,c)	x[v] - x[u] <= c
add_ge(u,v,c)	x[v] - x[u] >= c
add_lt(u,v,c)	x[v] - x[u] < c
add_gt(u,v,c)	x[v] - x[u] > c
add_eq(u,v,c)	x[v] - x[u] == c
add_range(u,v,lower,upper)	lower <= x[v] - x[u] <= upper
add_upper_bound(v,c)	x[v] <= c
add_lower_bound(v,c)	x[v] >= c
fix(v,c)	x[v] == c
add_range(v,lower,upper)	lower <= x[v] <= upper

制約追加は ID を std::size_t で返し、等式・区間・固定値は 2 個の ID を std::pair で返します。constraint_count() はこの本数、constraints()[id] は x[to] - x[from] <= weight を表す {from,to,weight} です。reserve(count) にはこの本数の予約量を指定します。

solve() は Result を返します。結果は呼び出した時点の系を保持し、その後の制約追加に影響されません。

結果の操作	内容・返り値
feasible()	全制約を同時に満たせるか
size() / origin()	変数数／定数 0 の補助頂点
values()	実現可能な解を const vector<value_type>& で取得。通常の変数のみ
value(v)	解の 1 変数。value(origin()) は 0
negative_cycle()	矛盾を示す正規化済み制約の閉路。解があれば空
upper_bounds(source)	全変数について x[v] - x[source] の最大値
upper_bounds()	全変数について x[v] の最大値
max_difference(u,v)	x[v] - x[u] の最大値
min_difference(u,v)	x[v] - x[u] の最小値
upper_bound(v) / lower_bound(v)	x[v] の最大値／最小値

解がない場合は解・上下限を取得しないでください。negative_cycle() は元の制約辺を閉路の順に返し、その重みの総和は負になります。

values() は解の一つで、各変数を最大化した解ではありません。最大値・最小値は std::optional<value_type> で返し、std::nullopt はその方向に非有界であることを表します。有限な値は実際に取り得る極値です。upper_bounds は通常の変数に対応する vector<optional<value_type>> を返します。

同じ始点から複数の差を調べる場合は upper_bounds(source) を一度呼び、その結果を参照してください。max_difference などの個別の問い合わせは、呼び出すたびに計算を行います。

牛ゲーでは、制約を x[v] - x[u] <= c の形で表し、変数間の差の上限を最短路として求めます。この形にそろえた全ての制約で c >= 0 なら、Dijkstra 法をそのまま使えます。 add_le(u,v,c) のみで制約を追加する場合は、全ての c が非負であることが条件です。origin() に関する制約も含めます。

この条件を満たす場合、solve() は \(O(N+M)\) 時間で済み、その後の上下限問い合わせは Dijkstra 法と同じ計算量です。負の c を含む一般の場合でも、solve() が成功した後の上下限問い合わせは同じ計算量になります。

変数と origin() を合わせた頂点数を \(N\)、constraint_count() を \(M\) とします。64 bit の計算量では \(w=64\) とします。

操作・条件	64 bit	128 bit
制約追加	償却 \(O(1)\)	償却 \(O(1)\)
solve()：一般の場合	最悪 \(O(N+NM)\)	最悪 \(O(N+NM+M\log M)\)
solve()：上記の全ての c が非負の場合	\(O(N+M)\)	\(O(N+M)\)
上下限・差の最大値・最小値の取得（1 回）	\(O(N+Mw)\)	\(O(N+M\log M)\)

solve() と上下限の取得は \(O(N+M)\) 追加空間を使います。同じ変数同士の差の取得と、解・負閉路の参照は \(O(1)\) 時間です。

依存ライブラリ

• Radix Heap
• Dijkstra の内部キュー
• Dijkstra Topological
• Dijkstra の一意経路の探索
• 有向グラフの CSR
• Shortest Path SCC
• Bellman-Ford
• Dijkstra

検証コード

• test/difference_constraints.test.cpp

コード

#ifndef DIFFERENCE_CONSTRAINTS_HPP
#define DIFFERENCE_CONSTRAINTS_HPP

#include "bellman_ford.hpp"
#include "dijkstra.hpp"
#include <optional>

template<class Cost = std::int64_t>
class DifferenceConstraints {
	static_assert(std::numeric_limits<Cost>::is_integer && std::numeric_limits<Cost>::is_signed && sizeof(Cost) <= 16);

  public:
	using value_type = std::conditional_t<(sizeof(Cost) <= 8), std::int64_t, __int128_t>;
	using Constraint = DirectedCSR<value_type>::Edge;
	class Result {
		friend class DifferenceConstraints;
		int n_;
		DirectedCSR<value_type> graph_;
		std::vector<value_type> values_;
		std::vector<Constraint> negative_cycle_;
		explicit Result(int n) : n_(n) {}
		void check([[maybe_unused]] int v) const { assert(0 <= v && v <= n_); }
		value_type potential(int v) const { return v == n_ ? 0 : values_[v]; }
		template<class Function>
		auto with_distances(int source, Function &&function) const {
			assert(feasible());
			check(source);
			if (!n_) return function(std::vector<value_type>{0});
			auto distance = dijkstra(graph_, source, std::numeric_limits<value_type>::max());
			if (std::find(distance.begin(), distance.end(), value_type(-1)) != distance.end())
				for (int v = 0; v <= n_; ++v) if (distance[v] != -1)
					for (auto [to, weight] : graph_[v])
						if (distance[to] == -1) throw std::overflow_error("DifferenceConstraints distance");
			return function(distance);
		}

	  public:
		Result(const Result &) = default;
		Result &operator=(const Result &) = default;
		Result(Result &&other) noexcept : n_(std::exchange(other.n_, 0)), graph_(std::move(other.graph_)),
			values_(std::move(other.values_)), negative_cycle_(std::move(other.negative_cycle_)) {}
		Result &operator=(Result &&other) noexcept {
			if (this != &other) {
				n_ = std::exchange(other.n_, 0);
				graph_ = std::move(other.graph_);
				values_ = std::move(other.values_);
				negative_cycle_ = std::move(other.negative_cycle_);
			}
			return *this;
		}
		int size() const { return n_; }
		int origin() const { return n_; }
		bool feasible() const { return negative_cycle_.empty(); }
		const std::vector<value_type> &values() const { assert(feasible()); return values_; }
		value_type value(int v) const { assert(feasible()); check(v); return potential(v); }
		const std::vector<Constraint> &negative_cycle() const { return negative_cycle_; }
		std::vector<std::optional<value_type>> upper_bounds(int source) const {
			return with_distances(source, [&](const auto &reduced) {
				std::vector<std::optional<value_type>> result(n_);
				for (int v = 0; v < n_; ++v)
					if (reduced[v] != -1) result[v] = add_sub(reduced[v], potential(v), potential(source));
				return result;
			});
		}
		std::vector<std::optional<value_type>> upper_bounds() const { return upper_bounds(origin()); }
		std::optional<value_type> max_difference(int u, int v) const {
			assert(feasible());
			check(u), check(v);
			if (u == v) return 0;
			return with_distances(u, [&](const auto &reduced) -> std::optional<value_type> {
				if (reduced[v] == -1) return std::nullopt;
				return add_sub(reduced[v], potential(v), potential(u));
			});
		}
		std::optional<value_type> min_difference(int u, int v) const {
			auto result = max_difference(v, u);
			if (result) *result = checked_sub(0, *result);
			return result;
		}
		std::optional<value_type> upper_bound(int v) const { return max_difference(origin(), v); }
		std::optional<value_type> lower_bound(int v) const { return min_difference(origin(), v); }
	};

  private:
	int n_;
	std::vector<Constraint> constraints_;
	void check([[maybe_unused]] int v) const { assert(0 <= v && v <= n_); }
	std::size_t add(int u, int v, value_type bound) {
		check(u), check(v);
		constraints_.push_back({u, v, bound});
		return constraints_.size() - 1;
	}
	static value_type checked_add(value_type a, value_type b) {
		value_type result;
		if (__builtin_add_overflow(a, b, &result)) throw std::overflow_error("DifferenceConstraints addition");
		return result;
	}
	static value_type checked_sub(value_type a, value_type b) {
		value_type result;
		if (__builtin_sub_overflow(a, b, &result)) throw std::overflow_error("DifferenceConstraints subtraction");
		return result;
	}
	static value_type add_sub(value_type a, value_type b, value_type c) {
		if ((b < 0) == (c < 0)) return checked_add(a, checked_sub(b, c));
		return checked_sub(checked_add(a, b), c);
	}

  public:
	explicit DifferenceConstraints(int n = 0) : n_(n) { assert(0 <= n && n <= std::numeric_limits<int>::max() - 2); }
	int size() const { return n_; }
	int origin() const { return n_; }
	std::size_t constraint_count() const { return constraints_.size(); }
	const std::vector<Constraint> &constraints() const { return constraints_; }
	void reserve(std::size_t count) { constraints_.reserve(count); }
	std::size_t add_le(int u, int v, Cost bound) { return add(u, v, value_type(bound)); }
	std::size_t add_ge(int u, int v, Cost bound) { return add(v, u, checked_sub(0, value_type(bound))); }
	std::size_t add_lt(int u, int v, Cost bound) { return add(u, v, checked_sub(value_type(bound), 1)); }
	std::size_t add_gt(int u, int v, Cost bound) { return add(v, u, checked_sub(-1, value_type(bound))); }
	std::pair<std::size_t, std::size_t> add_eq(int u, int v, Cost value) {
		value_type reverse = checked_sub(0, value_type(value));
		return {add(u, v, value), add(v, u, reverse)};
	}
	std::pair<std::size_t, std::size_t> add_range(int u, int v, Cost lower, Cost upper) {
		value_type reverse = checked_sub(0, value_type(lower));
		return {add(v, u, reverse), add(u, v, upper)};
	}
	std::size_t add_upper_bound(int v, Cost bound) { return add_le(origin(), v, bound); }
	std::size_t add_lower_bound(int v, Cost bound) { return add_ge(origin(), v, bound); }
	std::pair<std::size_t, std::size_t> fix(int v, Cost value) { return add_eq(origin(), v, value); }
	std::pair<std::size_t, std::size_t> add_range(int v, Cost lower, Cost upper) { return add_range(origin(), v, lower, upper); }
	Result solve() const {
		Result result(n_);
		std::vector<Constraint> edges = constraints_;
		edges.reserve(constraints_.size() + std::size_t(n_) + 1);
		const int source = n_ + 1;
		for (int v = 0; v <= n_; ++v) edges.push_back({source, v, 0});
		{
			DirectedCSR<value_type> graph(n_ + 2, edges);
			auto found = bellman_ford<value_type>(graph, source);
			if (found.has_negative_cycle()) {
				for (std::size_t i = 0; i < found.negative_cycle_edges.size(); ++i) {
					int from = found.negative_cycle[i];
					auto [to, weight] = graph[from][found.negative_cycle_edges[i]];
					assert(from <= n_ && to <= n_);
					result.negative_cycle_.push_back({from, to, weight});
				}
				return result;
			}
			result.values_.resize(n_);
			for (int v = 0; v < n_; ++v) result.values_[v] = checked_sub(found.distance[v], found.distance[n_]);
		}
		edges.resize(constraints_.size());
		for (auto &edge : edges) {
			edge.weight = add_sub(edge.weight, result.potential(edge.from), result.potential(edge.to));
			assert(edge.weight >= 0);
			if (edge.weight == std::numeric_limits<value_type>::max()) throw std::overflow_error("DifferenceConstraints reduced weight");
		}
		result.graph_ = DirectedCSR<value_type>(n_ + 1, edges);
		return result;
	}
};

#endif

Back to top page