Stern-Brocot Tree
(math/stern_brocot_tree.hpp)

View this file on GitHub
Last update: 2026-10-05 22:23:07+09:00
Include: #include "math/stern_brocot_tree.hpp"

Overview

The Stern-Brocot tree contains every positive reduced rational exactly once. It is useful when rational numbers must be handled by their position in the tree instead of by floating point values: comparing continued-fraction-like paths, finding the lowest common ancestor of two fractions, walking to a parent/child, or representing a whole interval of fractions by one node.

The root is 1/1. Each node represents the mediant of two boundary fractions L and R; the left child replaces the right boundary by the current fraction, and the right child replaces the left boundary by the current fraction.

This library stores a path from the root in run-length form, such as L^3 R^2, rather than one direction per edge. Encoding a fraction is essentially the Euclidean algorithm, so very deep paths such as 1/10^18 are represented by one run instead of 10^18 - 1 individual moves.

All fraction-returning operations use Rational<T>.

Model

The root path is empty and decodes to 1/1.

For a path, stern_brocot_bounds(path) returns the open interval containing exactly the fractions in the subtree rooted at that path. Boundaries are stored as (numerator, denominator) pairs because the right boundary of the whole tree is positive infinity, represented by 1/0.

For example, the path to 5/3 is R L R:

1/1 --R--> 2/1 --L--> 3/2 --R--> 5/3

Path Types

SternBrocotDirection has values Left and Right.

SternBrocotRun stores a direction and positive repetition count.

SternBrocotPath exposes:

Method	Description
`bool empty() const`	Whether the path denotes the root.
`uint64_t depth() const`	Expanded edge count.
`void push(direction, count)`	Appends moves and merges equal adjacent directions.
`bool move_up(count)`	Removes edges; returns false if the ancestor is above the root.
`SternBrocotPath ancestor(count) const`	Returns an ancestor path; requires sufficient depth.

push accepts count = 0 and does nothing. A valid stored run has positive count; decoding or bound construction asserts this condition.

Functions

All functions are in namespace m1une::math. Template parameter T must be a signed integer type supported by Rational<T>, such as int or long long.

Function signature	Description	Complexity
`template<class T>` `SternBrocotPath stern_brocot_path(T p, T q)`	Encodes the positive reduced fraction `p/q`.	`O(log max(p,q))` runs
`template<class T = long long>` `Rational<T> stern_brocot_decode(const SternBrocotPath& path)`	Decodes a path.	`O(number of runs)`
`template<class T>` `uint64_t stern_brocot_depth(T p, T q)`	Returns the expanded tree depth.	`O(log max(p,q))`
`SternBrocotPath stern_brocot_lca_path(const SternBrocotPath& a, const SternBrocotPath& b)`	Returns the compressed common path prefix.	`O(number of runs)`
`template<class T>` `Rational<T> stern_brocot_lca(T a, T b, T c, T d)`	Returns the LCA of `a/b` and `c/d`.	`O(log max values)`
`template<class T>` `Rational<T> stern_brocot_ancestor(T p, T q, uint64_t up)`	Returns the ancestor `up` edges above `p/q`.	`O(number of runs)`
`template<class T>` `Rational<T> stern_brocot_parent(T p, T q)`	Returns the parent.	`O(log max(p,q))`
`template<class T>` `Rational<T> stern_brocot_move(T p, T q, SternBrocotDirection direction, uint64_t count = 1)`	Descends repeatedly.	`O(log max(p,q))`
`template<class T = long long>` `SternBrocotBounds<T> stern_brocot_bounds(const SternBrocotPath& path)`	Returns the open interval boundaries of the path subtree.	`O(number of runs)`

Here O(log max(p,q)) means the number of Euclidean divisions, not the expanded tree depth. For example, stern_brocot_path(1, 1000000000000000000) creates one run.

SternBrocotBounds<T> stores:

Member	Type	Meaning
`left`	`std::pair<T, T>`	Left boundary fraction `(numerator, denominator)`.
`right`	`std::pair<T, T>`	Right boundary fraction `(numerator, denominator)`.

Contracts

Inputs to stern_brocot_path and functions that encode a fraction must satisfy:

p > 0
q > 0
gcd(p, q) = 1

These conditions are asserted in debug builds.

Additional preconditions:

Decoded numerators, denominators, and interval boundaries must fit in T.
stern_brocot_ancestor(p, q, up) requires up <= stern_brocot_depth(p, q).
stern_brocot_parent(p, q) requires p/q != 1/1.
A manually created SternBrocotPath must contain only positive-count runs.

These are CP-style direct-return functions. Invalid calls are programmer errors and are checked with assertions where practical.

Example

#include "math/stern_brocot_tree.hpp"
#include <iostream>

int main() {
    auto path = m1une::math::stern_brocot_path(5LL, 3LL);
    std::cout << path.depth() << '\n';  // 3

    auto same = m1une::math::stern_brocot_decode<long long>(path);
    std::cout << same << '\n';  // 5/3

    auto parent = m1une::math::stern_brocot_parent(5LL, 3LL);
    std::cout << parent << '\n';  // 3/2

    auto lca = m1une::math::stern_brocot_lca(5LL, 3LL, 7LL, 4LL);
    std::cout << lca << '\n';  // 5/3

    auto bounds = m1une::math::stern_brocot_bounds<long long>(path);
    std::cout << bounds.left.first << '/' << bounds.left.second << '\n';    // 3/2
    std::cout << bounds.right.first << '/' << bounds.right.second << '\n';  // 2/1
}

Notes

The order of paths is not lexicographic order of fractions. Use this library when tree structure matters; use Rational<T> comparisons when numeric order matters.
The subtree interval is open: the boundary fractions themselves are not inside the subtree.
Run-length paths make ancestor and LCA operations fast even for fractions with enormous expanded depth.

Depends on

Rational Number (math/rational.hpp)

Required by

Math All (math/all.hpp)

Verified with

Code

#ifndef M1UNE_MATH_STERN_BROCOT_TREE_HPP
#define M1UNE_MATH_STERN_BROCOT_TREE_HPP 1

#include <cassert>
#include <concepts>
#include <cstdint>
#include <limits>
#include <utility>
#include <vector>

#include "rational.hpp"

namespace m1une {
namespace math {

enum class SternBrocotDirection {
    Left,
    Right,
};

struct SternBrocotRun {
    SternBrocotDirection direction;
    uint64_t count;

    friend bool operator==(const SternBrocotRun&, const SternBrocotRun&) = default;
};

struct SternBrocotPath {
    std::vector<SternBrocotRun> runs;

    bool empty() const {
        return runs.empty();
    }

    uint64_t depth() const {
        uint64_t result = 0;
        for (const SternBrocotRun& run : runs) {
            assert(run.count <= std::numeric_limits<uint64_t>::max() - result);
            result += run.count;
        }
        return result;
    }

    void push(SternBrocotDirection direction, uint64_t count = 1) {
        if (count == 0) return;
        if (!runs.empty() && runs.back().direction == direction) {
            assert(count <= std::numeric_limits<uint64_t>::max() - runs.back().count);
            runs.back().count += count;
        } else {
            runs.push_back(SternBrocotRun{direction, count});
        }
    }

    bool move_up(uint64_t count = 1) {
        if (depth() < count) return false;
        while (count > 0) {
            SternBrocotRun& run = runs.back();
            uint64_t removed = run.count < count ? run.count : count;
            run.count -= removed;
            count -= removed;
            if (run.count == 0) runs.pop_back();
        }
        return true;
    }

    SternBrocotPath ancestor(uint64_t count) const {
        SternBrocotPath result = *this;
        [[maybe_unused]] bool valid = result.move_up(count);
        assert(valid);
        return result;
    }

    friend bool operator==(const SternBrocotPath&, const SternBrocotPath&) = default;
};

template <std::signed_integral T = long long>
struct SternBrocotBounds {
    std::pair<T, T> left;
    std::pair<T, T> right;
};

template <std::signed_integral T>
SternBrocotPath stern_brocot_path(T numerator, T denominator) {
    assert(0 < numerator);
    assert(0 < denominator);
    SternBrocotPath result;
    while (numerator != denominator) {
        if (numerator < denominator) {
            T count = (denominator - 1) / numerator;
            result.push(SternBrocotDirection::Left, uint64_t(count));
            denominator -= count * numerator;
        } else {
            T count = (numerator - 1) / denominator;
            result.push(SternBrocotDirection::Right, uint64_t(count));
            numerator -= count * denominator;
        }
    }
    assert(numerator == 1);
    return result;
}

template <std::signed_integral T = long long>
Rational<T> stern_brocot_decode(const SternBrocotPath& path) {
    using Wide = __int128_t;
    Wide left_numerator = 0;
    Wide left_denominator = 1;
    Wide right_numerator = 1;
    Wide right_denominator = 0;
    Wide numerator = 1;
    Wide denominator = 1;
    [[maybe_unused]] const Wide maximum = std::numeric_limits<T>::max();
    for (const SternBrocotRun& run : path.runs) {
        assert(run.count > 0);
        assert(Wide(run.count) <= maximum);
        Wide count = run.count;
        if (run.direction == SternBrocotDirection::Left) {
            right_numerator = numerator + (count - 1) * left_numerator;
            right_denominator = denominator + (count - 1) * left_denominator;
            numerator += count * left_numerator;
            denominator += count * left_denominator;
        } else {
            left_numerator = numerator + (count - 1) * right_numerator;
            left_denominator = denominator + (count - 1) * right_denominator;
            numerator += count * right_numerator;
            denominator += count * right_denominator;
        }
        assert(numerator <= maximum);
        assert(denominator <= maximum);
        assert(left_numerator <= maximum);
        assert(left_denominator <= maximum);
        assert(right_numerator <= maximum);
        assert(right_denominator <= maximum);
    }
    return Rational<T>(T(numerator), T(denominator));
}

template <std::signed_integral T>
uint64_t stern_brocot_depth(T numerator, T denominator) {
    return stern_brocot_path(numerator, denominator).depth();
}

inline SternBrocotPath stern_brocot_lca_path(
    const SternBrocotPath& first,
    const SternBrocotPath& second
) {
    SternBrocotPath result;
    int limit = int(first.runs.size() < second.runs.size() ? first.runs.size() : second.runs.size());
    for (int i = 0; i < limit; i++) {
        if (first.runs[i].direction != second.runs[i].direction) break;
        uint64_t common =
            first.runs[i].count < second.runs[i].count
                ? first.runs[i].count
                : second.runs[i].count;
        result.push(first.runs[i].direction, common);
        if (first.runs[i].count != second.runs[i].count) break;
    }
    return result;
}

template <std::signed_integral T>
Rational<T> stern_brocot_lca(
    T first_numerator,
    T first_denominator,
    T second_numerator,
    T second_denominator
) {
    SternBrocotPath path = stern_brocot_lca_path(
        stern_brocot_path(first_numerator, first_denominator),
        stern_brocot_path(second_numerator, second_denominator)
    );
    return stern_brocot_decode<T>(path);
}

template <std::signed_integral T>
Rational<T> stern_brocot_ancestor(T numerator, T denominator, uint64_t up) {
    SternBrocotPath path = stern_brocot_path(numerator, denominator);
    [[maybe_unused]] bool valid = path.move_up(up);
    assert(valid);
    return stern_brocot_decode<T>(path);
}

template <std::signed_integral T>
Rational<T> stern_brocot_parent(T numerator, T denominator) {
    return stern_brocot_ancestor(numerator, denominator, 1);
}

template <std::signed_integral T>
Rational<T> stern_brocot_move(
    T numerator,
    T denominator,
    SternBrocotDirection direction,
    uint64_t count = 1
) {
    SternBrocotPath path = stern_brocot_path(numerator, denominator);
    path.push(direction, count);
    return stern_brocot_decode<T>(path);
}

template <std::signed_integral T = long long>
SternBrocotBounds<T> stern_brocot_bounds(const SternBrocotPath& path) {
    using Wide = __int128_t;
    Wide left_numerator = 0;
    Wide left_denominator = 1;
    Wide right_numerator = 1;
    Wide right_denominator = 0;
    Wide numerator = 1;
    Wide denominator = 1;
    [[maybe_unused]] const Wide maximum = std::numeric_limits<T>::max();

    for (const SternBrocotRun& run : path.runs) {
        assert(run.count > 0);
        assert(Wide(run.count) <= maximum);
        Wide count = run.count;
        if (run.direction == SternBrocotDirection::Left) {
            right_numerator = numerator + (count - 1) * left_numerator;
            right_denominator = denominator + (count - 1) * left_denominator;
            numerator += count * left_numerator;
            denominator += count * left_denominator;
        } else {
            left_numerator = numerator + (count - 1) * right_numerator;
            left_denominator = denominator + (count - 1) * right_denominator;
            numerator += count * right_numerator;
            denominator += count * right_denominator;
        }
        assert(numerator <= maximum);
        assert(denominator <= maximum);
        assert(left_numerator <= maximum);
        assert(left_denominator <= maximum);
        assert(right_numerator <= maximum);
        assert(right_denominator <= maximum);
    }
    SternBrocotBounds<T> result;
    result.left = {T(left_numerator), T(left_denominator)};
    result.right = {T(right_numerator), T(right_denominator)};
    return result;
}

}  // namespace math
}  // namespace m1une

#endif  // M1UNE_MATH_STERN_BROCOT_TREE_HPP

#line 1 "math/stern_brocot_tree.hpp"



#include <cassert>
#include <concepts>
#include <cstdint>
#include <limits>
#include <utility>
#include <vector>

#line 1 "math/rational.hpp"



#include <algorithm>
#line 6 "math/rational.hpp"
#include <cmath>
#include <compare>
#line 9 "math/rational.hpp"
#include <iostream>
#line 11 "math/rational.hpp"
#include <sstream>
#include <string>
#include <type_traits>
#line 15 "math/rational.hpp"

namespace m1une {
namespace math {

namespace rational_detail {

template <class T>
concept IntegerLike =
    std::signed_integral<T> ||
    (!std::integral<T> && std::copyable<T> && requires(T first, T second) {
        T(0);
        T(1);
        { -first } -> std::same_as<T>;
        { first + second } -> std::same_as<T>;
        { first - second } -> std::same_as<T>;
        { first * second } -> std::same_as<T>;
        { first / second } -> std::same_as<T>;
        { first % second } -> std::same_as<T>;
        { first += second } -> std::same_as<T&>;
        { first -= second } -> std::same_as<T&>;
        { first /= second } -> std::same_as<T&>;
        { first == second } -> std::convertible_to<bool>;
        { first < second } -> std::convertible_to<bool>;
    });

}  // namespace rational_detail

template <rational_detail::IntegerLike T = long long>
struct Rational {
    static_assert(!std::signed_integral<T> || sizeof(T) <= sizeof(long long));

   private:
    static constexpr bool BUILTIN_INTEGER = std::signed_integral<T>;
    using Wide = std::conditional_t<BUILTIN_INTEGER, __int128_t, T>;
    using Magnitude = std::conditional_t<BUILTIN_INTEGER, __uint128_t, T>;

    T _numerator;
    T _denominator;

    static constexpr Magnitude magnitude(Wide value) {
        if constexpr (BUILTIN_INTEGER) {
            if (value < 0) {
                return static_cast<Magnitude>(-(value + 1)) + 1;
            }
            return static_cast<Magnitude>(value);
        } else {
            return value < 0 ? -value : value;
        }
    }

    static constexpr Magnitude gcd(Magnitude first, Magnitude second) {
        while (second != 0) {
            Magnitude remainder = first % second;
            first = second;
            second = remainder;
        }
        return first;
    }

    static constexpr T narrow(Wide value) {
        if constexpr (BUILTIN_INTEGER) {
            assert(Wide(std::numeric_limits<T>::min()) <= value);
            assert(value <= Wide(std::numeric_limits<T>::max()));
            return static_cast<T>(value);
        } else {
            return value;
        }
    }

    constexpr void assign_normalized(Wide numerator, Wide denominator) {
        assert(denominator != 0);
        if (numerator == 0) {
            _numerator = 0;
            _denominator = 1;
            return;
        }

        Magnitude divisor = gcd(magnitude(numerator), magnitude(denominator));
        numerator /= static_cast<Wide>(divisor);
        denominator /= static_cast<Wide>(divisor);
        if (denominator < 0) {
            numerator = -numerator;
            denominator = -denominator;
        }
        _numerator = narrow(numerator);
        _denominator = narrow(denominator);
    }

    static constexpr Rational from_wide(Wide numerator, Wide denominator) {
        Rational result;
        result.assign_normalized(numerator, denominator);
        return result;
    }

    static std::pair<long double, long long> decimal_scientific(const T& value) {
        std::ostringstream output;
        output << value;
        const std::string text = output.str();
        std::size_t begin = 0;
        int sign = 1;
        if (!text.empty() && (text[0] == '-' || text[0] == '+')) {
            if (text[0] == '-') sign = -1;
            begin = 1;
        }
        while (begin < text.size() && text[begin] == '0') ++begin;
        if (begin == text.size()) return std::make_pair(0.0L, 0LL);

        constexpr int DIGITS = std::numeric_limits<long double>::digits10 + 1;
        const std::size_t used = std::min<std::size_t>(DIGITS, text.size() - begin);
        long double significand = 0;
        for (std::size_t i = 0; i < used; ++i) {
            assert('0' <= text[begin + i] && text[begin + i] <= '9');
            significand = significand * 10 + (text[begin + i] - '0');
        }
        for (std::size_t i = 1; i < used; ++i) significand /= 10;
        const long long exponent = static_cast<long long>(text.size() - begin - 1);
        return std::make_pair(sign * significand, exponent);
    }

   public:
    constexpr Rational() : _numerator(0), _denominator(1) {}

    constexpr Rational(T integer) : _numerator(integer), _denominator(1) {}

    template <std::integral U>
        requires std::constructible_from<T, U> &&
                 (!std::same_as<std::remove_cv_t<U>, T>)
    constexpr Rational(U integer) : Rational(T(integer)) {}

    constexpr Rational(T numerator, T denominator) {
        assign_normalized(Wide(numerator), Wide(denominator));
    }

    constexpr T numerator() const {
        return _numerator;
    }

    constexpr T denominator() const {
        return _denominator;
    }

    constexpr bool is_integer() const {
        return _denominator == 1;
    }

    constexpr int sign() const {
        return (_numerator > 0) - (_numerator < 0);
    }

    constexpr Rational reciprocal() const {
        assert(_numerator != 0);
        return from_wide(Wide(_denominator), Wide(_numerator));
    }

    constexpr Rational abs() const {
        return _numerator < 0 ? -*this : *this;
    }

    constexpr long double to_long_double() const
        requires requires(const T& value) { static_cast<long double>(value); }
    {
        return static_cast<long double>(_numerator) / static_cast<long double>(_denominator);
    }

    long double to_long_double() const
        requires(!requires(const T& value) { static_cast<long double>(value); })
    {
        const auto [numerator, numerator_exponent] = decimal_scientific(_numerator);
        const auto [denominator, denominator_exponent] = decimal_scientific(_denominator);
        return numerator / denominator *
               std::pow(10.0L, numerator_exponent - denominator_exponent);
    }

    template <std::floating_point F>
    explicit constexpr operator F() const
        requires requires(const T& value) { static_cast<long double>(value); }
    {
        return static_cast<F>(to_long_double());
    }

    template <std::floating_point F>
    explicit operator F() const
        requires(!requires(const T& value) { static_cast<long double>(value); })
    {
        return static_cast<F>(to_long_double());
    }

    constexpr T trunc() const {
        return _numerator / _denominator;
    }

    constexpr T floor() const {
        T quotient = _numerator / _denominator;
        if (_numerator < 0 && _numerator % _denominator != 0) quotient -= T(1);
        return quotient;
    }

    constexpr T ceil() const {
        T quotient = _numerator / _denominator;
        if (0 < _numerator && _numerator % _denominator != 0) quotient += T(1);
        return quotient;
    }

    constexpr Rational operator+() const {
        return *this;
    }

    constexpr Rational operator-() const {
        return from_wide(-Wide(_numerator), Wide(_denominator));
    }

    constexpr Rational& operator+=(const Rational& other) {
        Magnitude common =
            gcd(static_cast<Magnitude>(_denominator), static_cast<Magnitude>(other._denominator));
        Wide left_scale = Wide(other._denominator) / static_cast<Wide>(common);
        Wide right_scale = Wide(_denominator) / static_cast<Wide>(common);
        Wide numerator =
            Wide(_numerator) * left_scale + Wide(other._numerator) * right_scale;

        // With both operands already reduced, every factor shared by the new
        // numerator and denominator must divide `common`.  Restricting the
        // second gcd to that value avoids a full-size gcd against the product
        // of both denominators, which is especially important for BigInt.
        Magnitude reduction = common == Magnitude(1)
                                  ? Magnitude(1)
                                  : gcd(magnitude(numerator), common);
        if (reduction != Magnitude(1)) {
            numerator /= static_cast<Wide>(reduction);
        }
        Wide remaining_denominator = Wide(other._denominator);
        if (reduction != Magnitude(1)) {
            remaining_denominator /= static_cast<Wide>(reduction);
        }
        _numerator = narrow(numerator);
        _denominator = narrow(right_scale * remaining_denominator);
        return *this;
    }

    constexpr Rational& operator-=(const Rational& other) {
        return *this += -other;
    }

    constexpr Rational& operator*=(const Rational& other) {
        Magnitude first_gcd = gcd(magnitude(Wide(_numerator)), static_cast<Magnitude>(other._denominator));
        Magnitude second_gcd = gcd(magnitude(Wide(other._numerator)), static_cast<Magnitude>(_denominator));
        assign_normalized((Wide(_numerator) / static_cast<Wide>(first_gcd)) *
                              (Wide(other._numerator) / static_cast<Wide>(second_gcd)),
                          (Wide(_denominator) / static_cast<Wide>(second_gcd)) *
                              (Wide(other._denominator) / static_cast<Wide>(first_gcd)));
        return *this;
    }

    constexpr Rational& operator/=(const Rational& other) {
        return *this *= other.reciprocal();
    }

    friend constexpr Rational operator+(Rational left, const Rational& right) {
        return left += right;
    }

    friend constexpr Rational operator-(Rational left, const Rational& right) {
        return left -= right;
    }

    friend constexpr Rational operator*(Rational left, const Rational& right) {
        return left *= right;
    }

    friend constexpr Rational operator/(Rational left, const Rational& right) {
        return left /= right;
    }

    friend constexpr bool operator==(const Rational& left, const Rational& right) {
        return left._numerator == right._numerator && left._denominator == right._denominator;
    }

    friend constexpr std::strong_ordering operator<=>(const Rational& left, const Rational& right) {
        Wide first = Wide(left._numerator) * Wide(right._denominator);
        Wide second = Wide(right._numerator) * Wide(left._denominator);
        if (first < second) return std::strong_ordering::less;
        if (second < first) return std::strong_ordering::greater;
        return std::strong_ordering::equal;
    }

    friend std::ostream& operator<<(std::ostream& output, const Rational& value) {
        output << value._numerator;
        if (value._denominator != 1) {
            output << '/' << value._denominator;
        }
        return output;
    }

    friend std::istream& operator>>(std::istream& input, Rational& value) {
        std::string token;
        if (!(input >> token)) return input;

        std::size_t slash = token.find('/');
        if (slash != std::string::npos && token.find('/', slash + 1) != std::string::npos) {
            input.setstate(std::ios::failbit);
            return input;
        }

        T numerator = 0;
        T denominator = 1;
        std::istringstream numerator_input(token.substr(0, slash));
        if (!(numerator_input >> numerator) || numerator_input.peek() != std::char_traits<char>::eof()) {
            input.setstate(std::ios::failbit);
            return input;
        }
        if (slash != std::string::npos) {
            std::istringstream denominator_input(token.substr(slash + 1));
            if (!(denominator_input >> denominator) ||
                denominator_input.peek() != std::char_traits<char>::eof()) {
                input.setstate(std::ios::failbit);
                return input;
            }
        }
        value = Rational(numerator, denominator);
        return input;
    }
};

template <rational_detail::IntegerLike T>
constexpr Rational<T> abs(const Rational<T>& value) {
    return value.abs();
}

}  // namespace math
}  // namespace m1une

namespace std {

// Integer/rational common types already follow from implicit integer
// construction. Mixing a floating scalar explicitly chooses approximation.
template <m1une::math::rational_detail::IntegerLike T, floating_point F>
struct common_type<m1une::math::Rational<T>, F> {
    using type = long double;
};

template <floating_point F, m1une::math::rational_detail::IntegerLike T>
struct common_type<F, m1une::math::Rational<T>> {
    using type = long double;
};

}  // namespace std


#line 12 "math/stern_brocot_tree.hpp"

namespace m1une {
namespace math {

enum class SternBrocotDirection {
    Left,
    Right,
};

struct SternBrocotRun {
    SternBrocotDirection direction;
    uint64_t count;

    friend bool operator==(const SternBrocotRun&, const SternBrocotRun&) = default;
};

struct SternBrocotPath {
    std::vector<SternBrocotRun> runs;

    bool empty() const {
        return runs.empty();
    }

    uint64_t depth() const {
        uint64_t result = 0;
        for (const SternBrocotRun& run : runs) {
            assert(run.count <= std::numeric_limits<uint64_t>::max() - result);
            result += run.count;
        }
        return result;
    }

    void push(SternBrocotDirection direction, uint64_t count = 1) {
        if (count == 0) return;
        if (!runs.empty() && runs.back().direction == direction) {
            assert(count <= std::numeric_limits<uint64_t>::max() - runs.back().count);
            runs.back().count += count;
        } else {
            runs.push_back(SternBrocotRun{direction, count});
        }
    }

    bool move_up(uint64_t count = 1) {
        if (depth() < count) return false;
        while (count > 0) {
            SternBrocotRun& run = runs.back();
            uint64_t removed = run.count < count ? run.count : count;
            run.count -= removed;
            count -= removed;
            if (run.count == 0) runs.pop_back();
        }
        return true;
    }

    SternBrocotPath ancestor(uint64_t count) const {
        SternBrocotPath result = *this;
        [[maybe_unused]] bool valid = result.move_up(count);
        assert(valid);
        return result;
    }

    friend bool operator==(const SternBrocotPath&, const SternBrocotPath&) = default;
};

template <std::signed_integral T = long long>
struct SternBrocotBounds {
    std::pair<T, T> left;
    std::pair<T, T> right;
};

template <std::signed_integral T>
SternBrocotPath stern_brocot_path(T numerator, T denominator) {
    assert(0 < numerator);
    assert(0 < denominator);
    SternBrocotPath result;
    while (numerator != denominator) {
        if (numerator < denominator) {
            T count = (denominator - 1) / numerator;
            result.push(SternBrocotDirection::Left, uint64_t(count));
            denominator -= count * numerator;
        } else {
            T count = (numerator - 1) / denominator;
            result.push(SternBrocotDirection::Right, uint64_t(count));
            numerator -= count * denominator;
        }
    }
    assert(numerator == 1);
    return result;
}

template <std::signed_integral T = long long>
Rational<T> stern_brocot_decode(const SternBrocotPath& path) {
    using Wide = __int128_t;
    Wide left_numerator = 0;
    Wide left_denominator = 1;
    Wide right_numerator = 1;
    Wide right_denominator = 0;
    Wide numerator = 1;
    Wide denominator = 1;
    [[maybe_unused]] const Wide maximum = std::numeric_limits<T>::max();
    for (const SternBrocotRun& run : path.runs) {
        assert(run.count > 0);
        assert(Wide(run.count) <= maximum);
        Wide count = run.count;
        if (run.direction == SternBrocotDirection::Left) {
            right_numerator = numerator + (count - 1) * left_numerator;
            right_denominator = denominator + (count - 1) * left_denominator;
            numerator += count * left_numerator;
            denominator += count * left_denominator;
        } else {
            left_numerator = numerator + (count - 1) * right_numerator;
            left_denominator = denominator + (count - 1) * right_denominator;
            numerator += count * right_numerator;
            denominator += count * right_denominator;
        }
        assert(numerator <= maximum);
        assert(denominator <= maximum);
        assert(left_numerator <= maximum);
        assert(left_denominator <= maximum);
        assert(right_numerator <= maximum);
        assert(right_denominator <= maximum);
    }
    return Rational<T>(T(numerator), T(denominator));
}

template <std::signed_integral T>
uint64_t stern_brocot_depth(T numerator, T denominator) {
    return stern_brocot_path(numerator, denominator).depth();
}

inline SternBrocotPath stern_brocot_lca_path(
    const SternBrocotPath& first,
    const SternBrocotPath& second
) {
    SternBrocotPath result;
    int limit = int(first.runs.size() < second.runs.size() ? first.runs.size() : second.runs.size());
    for (int i = 0; i < limit; i++) {
        if (first.runs[i].direction != second.runs[i].direction) break;
        uint64_t common =
            first.runs[i].count < second.runs[i].count
                ? first.runs[i].count
                : second.runs[i].count;
        result.push(first.runs[i].direction, common);
        if (first.runs[i].count != second.runs[i].count) break;
    }
    return result;
}

template <std::signed_integral T>
Rational<T> stern_brocot_lca(
    T first_numerator,
    T first_denominator,
    T second_numerator,
    T second_denominator
) {
    SternBrocotPath path = stern_brocot_lca_path(
        stern_brocot_path(first_numerator, first_denominator),
        stern_brocot_path(second_numerator, second_denominator)
    );
    return stern_brocot_decode<T>(path);
}

template <std::signed_integral T>
Rational<T> stern_brocot_ancestor(T numerator, T denominator, uint64_t up) {
    SternBrocotPath path = stern_brocot_path(numerator, denominator);
    [[maybe_unused]] bool valid = path.move_up(up);
    assert(valid);
    return stern_brocot_decode<T>(path);
}

template <std::signed_integral T>
Rational<T> stern_brocot_parent(T numerator, T denominator) {
    return stern_brocot_ancestor(numerator, denominator, 1);
}

template <std::signed_integral T>
Rational<T> stern_brocot_move(
    T numerator,
    T denominator,
    SternBrocotDirection direction,
    uint64_t count = 1
) {
    SternBrocotPath path = stern_brocot_path(numerator, denominator);
    path.push(direction, count);
    return stern_brocot_decode<T>(path);
}

template <std::signed_integral T = long long>
SternBrocotBounds<T> stern_brocot_bounds(const SternBrocotPath& path) {
    using Wide = __int128_t;
    Wide left_numerator = 0;
    Wide left_denominator = 1;
    Wide right_numerator = 1;
    Wide right_denominator = 0;
    Wide numerator = 1;
    Wide denominator = 1;
    [[maybe_unused]] const Wide maximum = std::numeric_limits<T>::max();

    for (const SternBrocotRun& run : path.runs) {
        assert(run.count > 0);
        assert(Wide(run.count) <= maximum);
        Wide count = run.count;
        if (run.direction == SternBrocotDirection::Left) {
            right_numerator = numerator + (count - 1) * left_numerator;
            right_denominator = denominator + (count - 1) * left_denominator;
            numerator += count * left_numerator;
            denominator += count * left_denominator;
        } else {
            left_numerator = numerator + (count - 1) * right_numerator;
            left_denominator = denominator + (count - 1) * right_denominator;
            numerator += count * right_numerator;
            denominator += count * right_denominator;
        }
        assert(numerator <= maximum);
        assert(denominator <= maximum);
        assert(left_numerator <= maximum);
        assert(left_denominator <= maximum);
        assert(right_numerator <= maximum);
        assert(right_denominator <= maximum);
    }
    SternBrocotBounds<T> result;
    result.left = {T(left_numerator), T(left_denominator)};
    result.right = {T(right_numerator), T(right_denominator)};
    return result;
}

}  // namespace math
}  // namespace m1une

m1une's library

This documentation is automatically generated by online-judge-tools/verification-helper

Stern-Brocot Tree
(math/stern_brocot_tree.hpp)

Overview

Model

Path Types

Functions

Contracts

Example

Notes

Depends on

Required by

Verified with

Code

Stern-Brocot Tree (math/stern_brocot_tree.hpp)

Overview

Model

Path Types

Functions

Contracts

Example

Notes

Depends on

Required by

Verified with

Code

Stern-Brocot Tree
(math/stern_brocot_tree.hpp)