stern_brocot.hpp

SECTIONMath INCLUDEnoya/math/stern_brocot.hpp

在正有理数与 Stern–Brocot 树的左右路径间转换；用于分数区间、连分数式跳跃和有理逼近。

\[ \displaystyle \frac{a}{b}\oplus\frac{c}{d}=\frac{a+c}{b+d} \]

Complexity: Time: \(O(\log \max (\mathrm{numerator} , \mathrm{denominator}))\) for path and node queries; \(O(r)\) to decode a path of \(r\) runs. Space: \(O(r)\) for the returned path of \(r\) runs; \(O(1)\) for other queries.

AC 记录：stern_brocot_tree。

跳到代码 · GitHub ↗

Implementation

ReferenceSourceExpanded

当前头文件，省略 include guard；依赖见 #include。

/// @complexity Time: O(log max(numerator,denominator)) for path and node queries;
/// O(r) to decode a path of r runs.
/// Space: O(r) for the returned path of r runs; O(1) for other queries.

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <limits>
#include <numeric>
#include <optional>
#include <utility>
#include <vector>

namespace noya {

using stern_brocot_step = std::pair<char, std::uint64_t>;

namespace stern_brocot_detail {
// Consume one Euclidean run; the caller stops at the root num == den.
inline stern_brocot_step next(std::uint64_t &num, std::uint64_t &den) {
  if (num < den) {
    std::uint64_t cnt = (den - 1) / num;
    den -= cnt * num;
    return {'L', cnt};
  }
  std::uint64_t cnt = (num - 1) / den;
  num -= cnt * den;
  return {'R', cnt};
}

struct bounds {
  using u128 = unsigned __int128;
  // Farey neighbors; their mediant is the current node.
  u128 ln = 0, ld = 1, rn = 1, rd = 0;
  void add(char dir, std::uint64_t cnt) {
    if (dir == 'L') {
      rn += u128(cnt) * ln;
      rd += u128(cnt) * ld;
    } else {
      ln += u128(cnt) * rn;
      ld += u128(cnt) * rd;
    }
  }
  std::pair<std::uint64_t, std::uint64_t> value() const {
    constexpr u128 mx = std::numeric_limits<std::uint64_t>::max();
    assert(ln + rn <= mx && ld + rd <= mx);
    return {std::uint64_t(ln + rn), std::uint64_t(ld + rd)};
  }
};
} // namespace stern_brocot_detail

/// @brief Run-length encoded L/R path from 1/1 to a positive reduced fraction.
inline std::vector<stern_brocot_step> stern_brocot_path(std::uint64_t num,
                                                        std::uint64_t den) {
  assert(num >= 1 && den >= 1);
  assert(std::gcd(num, den) == 1);
  std::vector<stern_brocot_step> pat;
  while (num != den) pat.push_back(stern_brocot_detail::next(num, den));
  return pat;
}

/// @brief Decode a run-length encoded Stern-Brocot path to its reduced
/// fraction; asserts when a 64-bit numerator or denominator would overflow.
inline std::pair<std::uint64_t, std::uint64_t>
stern_brocot_fraction(const std::vector<stern_brocot_step> &pat) {
  using u128 = unsigned __int128;
  stern_brocot_detail::bounds cur;
  constexpr u128 mx = std::numeric_limits<std::uint64_t>::max();
  for (auto [dir, cnt] : pat) {
    assert((dir == 'L' || dir == 'R') && cnt >= 1);
    cur.add(dir, cnt);
    assert(cur.ln + cur.rn <= mx && cur.ld + cur.rd <= mx);
  }
  return cur.value();
}

/// @brief Return the Stern-Brocot node at the given depth on the root-to-node
/// path, or nullopt when the requested depth is greater than the node depth.
inline std::optional<std::pair<std::uint64_t, std::uint64_t>>
stern_brocot_ancestor(std::uint64_t dep, std::uint64_t num, std::uint64_t den) {
  assert(num >= 1 && den >= 1);
  assert(std::gcd(num, den) == 1);
  stern_brocot_detail::bounds cur;
  while (dep > 0 && num != den) {
    auto [dir, cnt] = stern_brocot_detail::next(num, den);
    cnt = std::min(dep, cnt);
    cur.add(dir, cnt);
    dep -= cnt;
  }
  if (dep > 0) return std::nullopt;
  return cur.value();
}

/// @brief Return the lowest common ancestor of positive reduced fractions a/b and c/d.
inline std::pair<std::uint64_t, std::uint64_t>
stern_brocot_lca(std::uint64_t a, std::uint64_t b, std::uint64_t c,
                 std::uint64_t d) {
  assert(a >= 1 && b >= 1 && c >= 1 && d >= 1);
  assert(std::gcd(a, b) == 1 && std::gcd(c, d) == 1);
  stern_brocot_detail::bounds cur;
  while (a != b && c != d) {
    auto [lhs, ls] = stern_brocot_detail::next(a, b);
    auto [rhs, rs] = stern_brocot_detail::next(c, d);
    if (lhs != rhs) break;
    cur.add(lhs, std::min(ls, rs));
    if (ls != rs) break;
  }
  return cur.value();
}

/// @brief Return the open rational interval containing precisely the
/// descendants of a positive reduced fraction. The boundaries may be 0/1 or
/// 1/0.
inline std::pair<std::pair<std::uint64_t, std::uint64_t>,
                 std::pair<std::uint64_t, std::uint64_t>>
stern_brocot_subtree_range(std::uint64_t num, std::uint64_t den) {
  using u128 = unsigned __int128;
  assert(num >= 1 && den >= 1);
  assert(std::gcd(num, den) == 1);
  stern_brocot_detail::bounds cur;
  while (num != den) {
    auto [dir, cnt] = stern_brocot_detail::next(num, den);
    cur.add(dir, cnt);
  }
  constexpr u128 mx = std::numeric_limits<std::uint64_t>::max();
  assert(cur.ln <= mx && cur.ld <= mx);
  assert(cur.rn <= mx && cur.rd <= mx);
  return {{std::uint64_t(cur.ln), std::uint64_t(cur.ld)},
          {std::uint64_t(cur.rn), std::uint64_t(cur.rd)}};
}

} // namespace noya

 // NOYA_STERN_BROCOT_HPP

#ifndef NOYA_STERN_BROCOT_HPP
#define NOYA_STERN_BROCOT_HPP 1

/// @complexity Time: O(log max(numerator,denominator)) for path and node queries;
/// O(r) to decode a path of r runs.
/// Space: O(r) for the returned path of r runs; O(1) for other queries.

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <limits>
#include <numeric>
#include <optional>
#include <utility>
#include <vector>

namespace noya {

using stern_brocot_step = std::pair<char, std::uint64_t>;

namespace stern_brocot_detail {
// Consume one Euclidean run; the caller stops at the root num == den.
inline stern_brocot_step next(std::uint64_t &num, std::uint64_t &den) {
  if (num < den) {
    std::uint64_t cnt = (den - 1) / num;
    den -= cnt * num;
    return {'L', cnt};
  }
  std::uint64_t cnt = (num - 1) / den;
  num -= cnt * den;
  return {'R', cnt};
}

struct bounds {
  using u128 = unsigned __int128;
  // Farey neighbors; their mediant is the current node.
  u128 ln = 0, ld = 1, rn = 1, rd = 0;
  void add(char dir, std::uint64_t cnt) {
    if (dir == 'L') {
      rn += u128(cnt) * ln;
      rd += u128(cnt) * ld;
    } else {
      ln += u128(cnt) * rn;
      ld += u128(cnt) * rd;
    }
  }
  std::pair<std::uint64_t, std::uint64_t> value() const {
    constexpr u128 mx = std::numeric_limits<std::uint64_t>::max();
    assert(ln + rn <= mx && ld + rd <= mx);
    return {std::uint64_t(ln + rn), std::uint64_t(ld + rd)};
  }
};
} // namespace stern_brocot_detail

/// @brief Run-length encoded L/R path from 1/1 to a positive reduced fraction.
inline std::vector<stern_brocot_step> stern_brocot_path(std::uint64_t num,
                                                        std::uint64_t den) {
  assert(num >= 1 && den >= 1);
  assert(std::gcd(num, den) == 1);
  std::vector<stern_brocot_step> pat;
  while (num != den) pat.push_back(stern_brocot_detail::next(num, den));
  return pat;
}

/// @brief Decode a run-length encoded Stern-Brocot path to its reduced
/// fraction; asserts when a 64-bit numerator or denominator would overflow.
inline std::pair<std::uint64_t, std::uint64_t>
stern_brocot_fraction(const std::vector<stern_brocot_step> &pat) {
  using u128 = unsigned __int128;
  stern_brocot_detail::bounds cur;
  constexpr u128 mx = std::numeric_limits<std::uint64_t>::max();
  for (auto [dir, cnt] : pat) {
    assert((dir == 'L' || dir == 'R') && cnt >= 1);
    cur.add(dir, cnt);
    assert(cur.ln + cur.rn <= mx && cur.ld + cur.rd <= mx);
  }
  return cur.value();
}

/// @brief Return the Stern-Brocot node at the given depth on the root-to-node
/// path, or nullopt when the requested depth is greater than the node depth.
inline std::optional<std::pair<std::uint64_t, std::uint64_t>>
stern_brocot_ancestor(std::uint64_t dep, std::uint64_t num, std::uint64_t den) {
  assert(num >= 1 && den >= 1);
  assert(std::gcd(num, den) == 1);
  stern_brocot_detail::bounds cur;
  while (dep > 0 && num != den) {
    auto [dir, cnt] = stern_brocot_detail::next(num, den);
    cnt = std::min(dep, cnt);
    cur.add(dir, cnt);
    dep -= cnt;
  }
  if (dep > 0) return std::nullopt;
  return cur.value();
}

/// @brief Return the lowest common ancestor of positive reduced fractions a/b and c/d.
inline std::pair<std::uint64_t, std::uint64_t>
stern_brocot_lca(std::uint64_t a, std::uint64_t b, std::uint64_t c,
                 std::uint64_t d) {
  assert(a >= 1 && b >= 1 && c >= 1 && d >= 1);
  assert(std::gcd(a, b) == 1 && std::gcd(c, d) == 1);
  stern_brocot_detail::bounds cur;
  while (a != b && c != d) {
    auto [lhs, ls] = stern_brocot_detail::next(a, b);
    auto [rhs, rs] = stern_brocot_detail::next(c, d);
    if (lhs != rhs) break;
    cur.add(lhs, std::min(ls, rs));
    if (ls != rs) break;
  }
  return cur.value();
}

/// @brief Return the open rational interval containing precisely the
/// descendants of a positive reduced fraction. The boundaries may be 0/1 or
/// 1/0.
inline std::pair<std::pair<std::uint64_t, std::uint64_t>,
                 std::pair<std::uint64_t, std::uint64_t>>
stern_brocot_subtree_range(std::uint64_t num, std::uint64_t den) {
  using u128 = unsigned __int128;
  assert(num >= 1 && den >= 1);
  assert(std::gcd(num, den) == 1);
  stern_brocot_detail::bounds cur;
  while (num != den) {
    auto [dir, cnt] = stern_brocot_detail::next(num, den);
    cur.add(dir, cnt);
  }
  constexpr u128 mx = std::numeric_limits<std::uint64_t>::max();
  assert(cur.ln <= mx && cur.ld <= mx);
  assert(cur.rn <= mx && cur.rd <= mx);
  return {{std::uint64_t(cur.ln), std::uint64_t(cur.ld)},
          {std::uint64_t(cur.rn), std::uint64_t(cur.rd)}};
}

} // namespace noya

#endif // NOYA_STERN_BROCOT_HPP

/// @complexity Time: O(log max(numerator,denominator)) for path and node queries;
/// O(r) to decode a path of r runs.
/// Space: O(r) for the returned path of r runs; O(1) for other queries.

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <limits>
#include <numeric>
#include <optional>
#include <utility>
#include <vector>

namespace noya {

using stern_brocot_step = std::pair<char, std::uint64_t>;

namespace stern_brocot_detail {
// Consume one Euclidean run; the caller stops at the root num == den.
inline stern_brocot_step next(std::uint64_t &num, std::uint64_t &den) {
  if (num < den) {
    std::uint64_t cnt = (den - 1) / num;
    den -= cnt * num;
    return {'L', cnt};
  }
  std::uint64_t cnt = (num - 1) / den;
  num -= cnt * den;
  return {'R', cnt};
}

struct bounds {
  using u128 = unsigned __int128;
  // Farey neighbors; their mediant is the current node.
  u128 ln = 0, ld = 1, rn = 1, rd = 0;
  void add(char dir, std::uint64_t cnt) {
    if (dir == 'L') {
      rn += u128(cnt) * ln;
      rd += u128(cnt) * ld;
    } else {
      ln += u128(cnt) * rn;
      ld += u128(cnt) * rd;
    }
  }
  std::pair<std::uint64_t, std::uint64_t> value() const {
    constexpr u128 mx = std::numeric_limits<std::uint64_t>::max();
    assert(ln + rn <= mx && ld + rd <= mx);
    return {std::uint64_t(ln + rn), std::uint64_t(ld + rd)};
  }
};
} // namespace stern_brocot_detail

/// @brief Run-length encoded L/R path from 1/1 to a positive reduced fraction.
inline std::vector<stern_brocot_step> stern_brocot_path(std::uint64_t num,
                                                        std::uint64_t den) {
  assert(num >= 1 && den >= 1);
  assert(std::gcd(num, den) == 1);
  std::vector<stern_brocot_step> pat;
  while (num != den) pat.push_back(stern_brocot_detail::next(num, den));
  return pat;
}

/// @brief Decode a run-length encoded Stern-Brocot path to its reduced
/// fraction; asserts when a 64-bit numerator or denominator would overflow.
inline std::pair<std::uint64_t, std::uint64_t>
stern_brocot_fraction(const std::vector<stern_brocot_step> &pat) {
  using u128 = unsigned __int128;
  stern_brocot_detail::bounds cur;
  constexpr u128 mx = std::numeric_limits<std::uint64_t>::max();
  for (auto [dir, cnt] : pat) {
    assert((dir == 'L' || dir == 'R') && cnt >= 1);
    cur.add(dir, cnt);
    assert(cur.ln + cur.rn <= mx && cur.ld + cur.rd <= mx);
  }
  return cur.value();
}

/// @brief Return the Stern-Brocot node at the given depth on the root-to-node
/// path, or nullopt when the requested depth is greater than the node depth.
inline std::optional<std::pair<std::uint64_t, std::uint64_t>>
stern_brocot_ancestor(std::uint64_t dep, std::uint64_t num, std::uint64_t den) {
  assert(num >= 1 && den >= 1);
  assert(std::gcd(num, den) == 1);
  stern_brocot_detail::bounds cur;
  while (dep > 0 && num != den) {
    auto [dir, cnt] = stern_brocot_detail::next(num, den);
    cnt = std::min(dep, cnt);
    cur.add(dir, cnt);
    dep -= cnt;
  }
  if (dep > 0) return std::nullopt;
  return cur.value();
}

/// @brief Return the lowest common ancestor of positive reduced fractions a/b and c/d.
inline std::pair<std::uint64_t, std::uint64_t>
stern_brocot_lca(std::uint64_t a, std::uint64_t b, std::uint64_t c,
                 std::uint64_t d) {
  assert(a >= 1 && b >= 1 && c >= 1 && d >= 1);
  assert(std::gcd(a, b) == 1 && std::gcd(c, d) == 1);
  stern_brocot_detail::bounds cur;
  while (a != b && c != d) {
    auto [lhs, ls] = stern_brocot_detail::next(a, b);
    auto [rhs, rs] = stern_brocot_detail::next(c, d);
    if (lhs != rhs) break;
    cur.add(lhs, std::min(ls, rs));
    if (ls != rs) break;
  }
  return cur.value();
}

/// @brief Return the open rational interval containing precisely the
/// descendants of a positive reduced fraction. The boundaries may be 0/1 or
/// 1/0.
inline std::pair<std::pair<std::uint64_t, std::uint64_t>,
                 std::pair<std::uint64_t, std::uint64_t>>
stern_brocot_subtree_range(std::uint64_t num, std::uint64_t den) {
  using u128 = unsigned __int128;
  assert(num >= 1 && den >= 1);
  assert(std::gcd(num, den) == 1);
  stern_brocot_detail::bounds cur;
  while (num != den) {
    auto [dir, cnt] = stern_brocot_detail::next(num, den);
    cur.add(dir, cnt);
  }
  constexpr u128 mx = std::numeric_limits<std::uint64_t>::max();
  assert(cur.ln <= mx && cur.ld <= mx);
  assert(cur.rn <= mx && cur.rd <= mx);
  return {{std::uint64_t(cur.ln), std::uint64_t(cur.ld)},
          {std::uint64_t(cur.rn), std::uint64_t(cur.rd)}};
}

} // namespace noya

 // NOYA_STERN_BROCOT_HPP