cp-library

競技プログラミング用 C++ ライブラリ

View the Project on GitHub qandori/cp-library

:heavy_check_mark: Stern–Brocot Tree

math/stern_brocot_tree.hpp · GitHub で見る

共通ページ

#include "lib/math/stern_brocot_tree.hpp"

using Tree = SternBrocotTree<>;
Tree node(5, 3);
auto path = node.path();
auto ancestor = node.ancestor(2);
auto common = node.lca(Tree(7, 4));
auto [lower, upper] = node.range();
bool ok = node.move_left(2);
auto [low, high] = Tree::approximate(5, 3, 2); // 1/1, 2/1

正の有理数を Stern–Brocot Tree の頂点として扱います。根は 1/1、左の子は親より小さく、右の子は親より大きい値です。上の 5/3 は根から R 1, L 1, R 1 と移動した頂点で、深さ 3、深さ 2 の祖先は 3/2 です。

関数・型 内容・返り値 計算量
SternBrocotTree<Int = std::uint64_t> 頂点。Int は bool を除く 64 bit 以下の標準整数型 —(型・定数)
value_type pair<Int,Int>。分子・分母の順 —(型・定数)
count_type Int に対応する符号なし整数型。深さと移動回数に使う —(型・定数)
distance_type Int が 32 bit 以下なら uint64_t、それ以外なら __uint128_t —(型・定数)
Run / path_type Run{direction,count} / vector<Run>。方向は 'L' または 'R'。Run は等値比較可能 —(型・定数)
SternBrocotTree() 根 1/1 を構築 \(O(1)\)
SternBrocotTree(numerator,denominator) / SternBrocotTree(value) 正の分子・分母から構築。約分されていなくてもよい \(O(\log M)\)
numerator() / denominator() / value() 現在の既約分数の分子・分母・組 \(O(1)\)
lower_bound() / upper_bound() / range() 部分木を表す開区間の下限・上限・両方の組 \(O(1)\)
depth() / is_root() 根からの辺数 / 根であるか \(O(1)\)
operator== / operator!= 同じ頂点であるか / 異なる頂点であるか \(O(1)\)
move_left(k = 1) / move_right(k = 1) 同じ方向へ k 回移動。結果が Int に収まらなければ false、元の頂点を保持 \(O(1)\)
move_up(k = 1) 親へ k 回移動。根を越える場合は false、元の頂点を保持 \(O(\log M)\)
left(k = 1) / right(k = 1) 元の頂点を変更せず、同じ方向へ k 回移動した頂点を optional で返す。表現不能なら nullopt \(O(1)\)
parent() 親を optional で返す。根なら nullopt \(O(1)\)
contains(value) value が自身を含む部分木に属するか。約分されていなくてもよい \(O(1)\)
is_ancestor(other) 自身が other の祖先であるか。同じ頂点でも true \(O(1)\)
path() 根から自身への経路を連長圧縮した path_type で返す \(O(\log M)\)
Tree::encode(numerator,denominator) 根から指定した正の有理数への圧縮経路を返す \(O(\log M)\)
Tree::encode(numerator,denominator,emit) 各区間で emit(direction,count) を呼ぶ。配列を確保しない \(O(\log M)\)
Tree::decode(path) span<const Run> から頂点を optional で返す。方向が不正、または結果が表現不能なら nullopt \(O(s)\)
ancestor(depth) 指定した深さの祖先を optional で返す。自身より深ければ nullopt \(O(\log M)\)
Tree::ancestor(value,depth) 指定した正の有理数の、指定した深さの祖先を optional で返す \(O(\log M)\)
lca(other) / Tree::lca(value1,value2) 最低共通祖先の頂点を返す \(O(\log M)\)
distance(other) 二頂点間の辺数を distance_type で返す \(O(\log M)\)
Tree::binary_search(bound,predicate) 分子・分母の共通上限を指定して有理数二分探索 判定 \(O(\log M)\) 回
Tree::binary_search(max_numerator,max_denominator,predicate) 分子・分母の上限を別々に指定して有理数二分探索 判定 \(O(\log M)\) 回
Tree::approximate(numerator,denominator,bound) 分子・分母がともに bound 以下の分数による最良の下側・上側近似を返す \(O(\log M)\)
Tree::approximate(numerator,denominator,max_numerator,max_denominator) 分子・分母の上限を別々に指定して最良近似を返す \(O(\log M)\)

分子・分母は正で Int に収まる値とします。構築時には自然に約分され、0 や無限大を頂点として持つことはできません。range() と approximate() の端点では、0 を {0,1}、無限大を {1,0} で表します。range() の区間には自身も含まれ、端点は含まれません。有限の有理数との四則演算には Fraction を使えます。

approximate() は正の有理数 \(x=\displaystyle\frac{\mathrm{numerator}}{\mathrm{denominator}}\) を挟む、最大の \(l\le x\) と最小の \(r\ge x\) を pair<value_type,value_type> で返します。返される分数は既約です。上限はすべて正とし、分子 0 または分母 0 の端点も候補に含みます。入力が上限内で表現できるときは両端ともその値になります。入力は約分されていなくても構いません。例えば approximate(10,6,5,3) は両端とも 5/3 です。

Tree::binary_search(max_numerator,max_denominator,predicate) は、上限内の正の既約分数について、predicate(p,q) が true となる最大の分数と false となる最小の分数を std::pair<value_type,value_type> で返します。Tree::binary_search(bound,predicate) では分子・分母の上限をともに bound とします。上限は正の Int、判定関数の引数は分子・分母の順の Int です。引数の既定値はありません。

判定は有理数の昇順で true から false に一度だけ切り替わる必要があります。すべて true・すべて false でも構いません。該当する正の分数がなければ下端は {0,1}、上端は {1,0} を返します。この二つは境界用の番兵で、判定関数には渡しません。判定関数が呼ばれるのは分子・分母がそれぞれ指定上限以下の正の既約分数だけです。同じ値に対する判定結果は変えないでください。呼び出し順序や回数そのものへの依存は避けます。

approximate() と異なり、返す両端は必ず異なります。境界ちょうどの値をどちらに含めるかは、判定関数の <・<= などで指定します。探索対象は正の有理数で、負の値には対応しません。頂点の状態を変更しない静的関数です。判定関数はコピーせず、例外を投げた場合はそのまま伝播します。

次の例では、浮動小数点を使わずに平方根を挟みます。

#include "lib/math/stern_brocot_tree.hpp"
#include <cassert>
#include <cstdint>

int main() {
	using Tree = SternBrocotTree<>;
	auto [lower, upper] = Tree::binary_search(100,
		[](std::uint64_t p, std::uint64_t q) { return p * p <= 2 * q * q; });
	assert((lower == Tree::value_type{41, 29}));
	assert((upper == Tree::value_type{99, 70}));

	auto [ok, ng] = Tree::binary_search(20, 10,
		[](std::uint64_t p, std::uint64_t q) { return 3 * p <= 5 * q; });
	assert((ok == Tree::value_type{5, 3}));
	assert((ng == Tree::value_type{17, 10}));
}

上限には Int の最大値まで指定できます。判定関数内の演算はオーバーフローしないようにしてください。上の平方根の例は上限 100 に対して安全ですが、一般の 64 bit 上限で p*p や 2*q*q が同じ型に収まるとは限りません。

移動回数・深さは count_type で指定します。0 回の移動は成功し、深さ 0 の祖先は根です。decode() は空の経路を根として扱い、回数 0 や同方向の連続区間も受け付けます。encode()・path() の返り値には回数 0 を含まず、方向は交互に並びます。

分子・分母は元の Int に収まる必要がありますが、中間積には幅の広い整数型を使います。64 bit 型の二頂点間の距離は 64 bit を超える場合があるため、distance_type は 128 bit です。Tree::ancestor() と Tree::lca() は入力の圧縮経路を配列にせず処理し、共通部分だけが必要なときは途中で打ち切ります。

分子・分母の最大値を \(M\)、圧縮経路の区間数を \(s\) として、根の構築・値取得・範囲取得・深さ・比較・包含判定・左右への移動・子・親の取得は \(O(1)\) です。有理数からの構築・encode()・path()・ancestor()・lca()・distance()・approximate()・一括 move_up() は \(O(\log M)\) 時間、decode() は \(O(s)\) 時間です。経路を配列で返す場合の空間は \(O(\log M)\)、それ以外の追加空間と頂点一つの保持は \(O(1)\) です。

有理数二分探索では \(M=\max(\mathrm{max\_numerator},\mathrm{max\_denominator})\) として、判定関数の呼び出し回数は \(O(\log M)\)、判定 1 回の時間を \(T\) とすれば全体で \(O(T\log M)\) 時間、追加空間は \(O(1)\) です。同方向への移動回数を指数探索してから二分探索するため、長い一本道を 1 頂点ずつ進みません。上限 1 の場合は定数時間です。

依存ライブラリ

検証コード

コード

#ifndef STERN_BROCOT_TREE_HPP
#define STERN_BROCOT_TREE_HPP

#include "../components/stern_brocot_path.hpp"
#include <algorithm>
#include <cassert>
#include <concepts>
#include <cstdint>
#include <limits>
#include <optional>
#include <span>
#include <type_traits>
#include <utility>
#include <vector>

template <std::integral Int = std::uint64_t>
requires (sizeof(Int) <= 8 && !std::is_same_v<Int, bool>)
class SternBrocotTree {
public:
	using value_type = std::pair<Int, Int>;
	using count_type = std::make_unsigned_t<Int>;
	using distance_type = std::conditional_t<(sizeof(Int) <= 4), std::uint64_t, __uint128_t>;
	struct Run {
		char direction;
		count_type count;
		bool operator==(const Run &) const = default;
	};
	using path_type = std::vector<Run>;
private:
	using U = count_type;
	using W = distance_type;
	using Path = stern_brocot_detail::Path<U>;
	Int ln = 0, ld = 1, rn = 1, rd = 0;
	U level = 0;
	void advance(bool right, U k) {
		if (right) { ln += k * U(rn); ld += k * U(rd); }
		else { rn += k * U(ln); rd += k * U(ld); }
		level += k;
	}
	bool move(bool right, U k) {
		W n = W(numerator()) + W(k) * U(right ? rn : ln);
		W d = W(denominator()) + W(k) * U(right ? rd : ld);
		if (n > W(std::numeric_limits<Int>::max()) || d > W(std::numeric_limits<Int>::max())) return false;
		advance(right, k); return true;
	}
public:
	SternBrocotTree() = default;
	SternBrocotTree(Int numerator, Int denominator) {
		assert(numerator > 0 && denominator > 0);
		Path path{U(numerator), U(denominator)};
		while (!path.empty()) { bool right = path.right(); advance(right, path.pop()); }
	}
	explicit SternBrocotTree(value_type value) : SternBrocotTree(value.first, value.second) {}
	Int numerator() const { return ln + rn; }
	Int denominator() const { return ld + rd; }
	value_type value() const { return {numerator(), denominator()}; }
	value_type lower_bound() const { return {ln, ld}; }
	value_type upper_bound() const { return {rn, rd}; }
	std::pair<value_type, value_type> range() const { return {lower_bound(), upper_bound()}; }
	U depth() const { return level; }
	bool is_root() const { return level == 0; }
	bool operator==(const SternBrocotTree &) const = default;
	bool move_left(U k = 1) { return move(false, k); }
	bool move_right(U k = 1) { return move(true, k); }
	bool move_up(U k = 1) {
		if (k > level) return false;
		level -= k;
		while (k) {
			bool right = ln >= rn;
			U step = k == 1 ? 1 : std::min<U>(k, right ? U(ln) / U(rn) : U(rd) / U(ld));
			if (right) { ln -= step * U(rn); ld -= step * U(rd); }
			else { rn -= step * U(ln); rd -= step * U(ld); }
			k -= step;
		}
		return true;
	}
	std::optional<SternBrocotTree> left(U k = 1) const {
		auto result = *this;
		if (!result.move_left(k)) return std::nullopt;
		return result;
	}
	std::optional<SternBrocotTree> right(U k = 1) const {
		auto result = *this;
		if (!result.move_right(k)) return std::nullopt;
		return result;
	}
	std::optional<SternBrocotTree> parent() const {
		auto result = *this;
		if (!result.move_up()) return std::nullopt;
		return result;
	}
	bool contains(value_type x) const {
		assert(x.first > 0 && x.second > 0);
		return W(ln) * U(x.second) < W(ld) * U(x.first) && W(x.first) * U(rd) < W(x.second) * U(rn);
	}
	bool is_ancestor(const SternBrocotTree &other) const { return contains(other.value()); }
	path_type path() const { return encode(numerator(), denominator()); }
	template <class Emit>
	static void encode(Int numerator, Int denominator, Emit &&emit) {
		assert(numerator > 0 && denominator > 0);
		Path path{U(numerator), U(denominator)};
		while (!path.empty()) { char direction = path.right() ? 'R' : 'L'; emit(direction, path.pop()); }
	}
	static path_type encode(Int numerator, Int denominator) {
		path_type result;
		encode(numerator, denominator, [&](char direction, U count) { result.push_back({direction, count}); });
		return result;
	}
	static std::optional<SternBrocotTree> decode(std::span<const Run> path) {
		SternBrocotTree result;
		for (auto [direction, count] : path) {
			if (direction != 'L' && direction != 'R') return std::nullopt;
			if (!result.move(direction == 'R', count)) return std::nullopt;
		}
		return result;
	}
	static std::optional<SternBrocotTree> ancestor(value_type x, U depth) {
		assert(x.first > 0 && x.second > 0);
		SternBrocotTree result;
		Path path{U(x.first), U(x.second)};
		while (depth && !path.empty()) {
			bool right = path.right();
			U step = std::min(depth, path.pop());
			result.advance(right, step); depth -= step;
		}
		if (depth) return std::nullopt;
		return result;
	}
	std::optional<SternBrocotTree> ancestor(U depth) const {
		if (depth > level) return std::nullopt;
		return ancestor(value(), depth);
	}
	static SternBrocotTree lca(value_type x, value_type y) {
		assert(x.first > 0 && x.second > 0 && y.first > 0 && y.second > 0);
		SternBrocotTree result;
		Path a{U(x.first), U(x.second)}, b{U(y.first), U(y.second)};
		while (!a.empty() && !b.empty() && a.right() == b.right()) {
			bool right = a.right();
			U k = a.pop(), l = b.pop();
			result.advance(right, std::min(k, l));
			if (k != l) break;
		}
		return result;
	}
	SternBrocotTree lca(const SternBrocotTree &other) const { return lca(value(), other.value()); }
	W distance(const SternBrocotTree &other) const { return W(level) + other.level - 2 * W(lca(other).level); }
	template <class Predicate>
	static std::pair<value_type, value_type> binary_search(Int max_numerator, Int max_denominator, Predicate &&predicate) {
		assert(max_numerator > 0 && max_denominator > 0);
		W n = U(max_numerator), d = U(max_denominator);
		W a = 0, b = 1, c = 1, e = 0;
		while (a + c <= n && b + e <= d) {
			bool side = bool(predicate(Int(a + c), Int(b + e)));
			W &x = side ? a : c, &y = side ? b : e;
			W u = side ? c : a, v = side ? e : b;
			W limit = u ? (n - x) / u : (d - y) / v;
			if (v) limit = std::min(limit, (d - y) / v);
			auto same = [&](W k) { return bool(predicate(Int(x + k * u), Int(y + k * v))) == side; };
			W low = 1, high = 1;
			while (low < limit) {
				high = low > limit / 2 ? limit : low * 2;
				if (!same(high)) break;
				low = high;
			}
			while (high - low > 1) {
				W mid = low + (high - low) / 2;
				if (same(mid)) low = mid;
				else high = mid;
			}
			x += low * u; y += low * v;
		}
		return {{Int(a), Int(b)}, {Int(c), Int(e)}};
	}
	template <class Predicate>
	static std::pair<value_type, value_type> binary_search(Int bound, Predicate &&predicate) {
		return binary_search(bound, bound, std::forward<Predicate>(predicate));
	}
	static std::pair<value_type, value_type> approximate(Int numerator, Int denominator, Int max_numerator, Int max_denominator) {
		assert(numerator > 0 && denominator > 0 && max_numerator > 0 && max_denominator > 0);
		U a = numerator, b = denominator, ln = 0, ld = 1, rn = 1, rd = 0;
		U n = max_numerator, d = max_denominator;
		while (a != b) {
			bool right = a > b;
			U k = right ? (a - 1) / b : (b - 1) / a;
			U step = k;
			if (right ? W(ln) + W(k) * rn > n || W(ld) + W(k) * rd > d
				: W(rn) + W(k) * ln > n || W(rd) + W(k) * ld > d) {
				step = right ? (n - ln) / rn : (d - rd) / ld;
				if (right && rd) step = std::min<U>(step, (d - ld) / rd);
				if (!right && ln) step = std::min<U>(step, (n - rn) / ln);
			}
			if (right) { ln += step * rn; ld += step * rd; a -= step * b; }
			else { rn += step * ln; rd += step * ld; b -= step * a; }
			if (step < k) break;
		}
		if (a == b && W(ln) + rn <= n && W(ld) + rd <= d) {
			value_type exact{Int(ln + rn), Int(ld + rd)};
			return {exact, exact};
		}
		return {{Int(ln), Int(ld)}, {Int(rn), Int(rd)}};
	}
	static std::pair<value_type, value_type> approximate(Int numerator, Int denominator, Int bound) {
		return approximate(numerator, denominator, bound, bound);
	}
};

#endif
Back to top page