rational_approximation.hpp

SECTIONMath INCLUDEnoya/rational_approximation.hpp

在分子、分母均有上界时，求给定有理数两侧最接近的最简分数。

\[ \displaystyle \frac{L}{R}<\frac{x}{y}<\frac{U}{V} \]

Complexity: Time: O(log N log max(x,y)). Space: O(1).

AC 记录：rational_approximation。

跳到代码 · GitHub ↗

Implementation

ReferenceSourceExpanded

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

/// @complexity Time: O(log N log max(x,y)).  Space: O(1).

#include <cassert>
#include <cstdint>
#include <utility>

namespace noya {

struct bounded_rational {
  std::int64_t num = 0;
  std::int64_t den = 1;
};

namespace rational_approximation_detail {

inline bool equal(bounded_rational a, bounded_rational b) {
  return static_cast<__int128>(a.num) * b.den ==
         static_cast<__int128>(a.den) * b.num;
}

inline bool less_equal(bounded_rational a, bounded_rational b) {
  return static_cast<__int128>(a.num) * b.den <=
         static_cast<__int128>(a.den) * b.num;
}

inline bool bounded(bounded_rational val, std::int64_t lim) {
  return val.num <= lim && val.den <= lim;
}

inline bounded_rational multiply_add(bounded_rational a, std::int64_t mul,
                                     bounded_rational b) {
  return {a.num * mul + b.num, a.den * mul + b.den};
}

template <class Predicate> std::int64_t maximum_true(Predicate pre) {
  std::int64_t acc = 0;
  std::int64_t rej = 1;
  while (pre(rej)) {
    rej *= 2;
  }
  while (rej - acc > 1) {
    std::int64_t mid = (acc + rej) / 2;
    if (pre(mid)) {
      acc = mid;
    } else {
      rej = mid;
    }
  }
  return acc;
}

} // namespace rational_approximation_detail

/// @brief Find the closest reduced fractions on both sides of `x/y` whose
/// numerators and denominators are at most `lim`.  Alternating jumps between
/// the two Stern--Brocot boundaries consumes one continued-fraction run at a
/// time; each run length is maximized by integer binary search.  Missing lower
/// and upper answers are represented by `0/1` and `1/0`.
inline std::pair<bounded_rational, bounded_rational>
rational_approximation(std::int64_t lim, std::int64_t x, std::int64_t y) {
  assert(lim >= 1 && x >= 1 && y >= 1);
  bounded_rational tar{x, y};
  bounded_rational lo{0, 1};
  bounded_rational hi{1, 0};
  while (true) {
    std::int64_t ls =
        rational_approximation_detail::maximum_true([&](std::int64_t cnt) {
          bounded_rational can =
              rational_approximation_detail::multiply_add(hi, cnt, lo);
          return rational_approximation_detail::bounded(can, lim) &&
                 rational_approximation_detail::less_equal(can, tar);
        });
    lo = rational_approximation_detail::multiply_add(hi, ls, lo);
    if (rational_approximation_detail::equal(lo, tar)) {
      hi = lo;
      break;
    }

    std::int64_t rs =
        rational_approximation_detail::maximum_true([&](std::int64_t cnt) {
          bounded_rational can =
              rational_approximation_detail::multiply_add(lo, cnt, hi);
          return rational_approximation_detail::bounded(can, lim) &&
                 rational_approximation_detail::less_equal(tar, can);
        });
    hi = rational_approximation_detail::multiply_add(lo, rs, hi);
    if (rational_approximation_detail::equal(hi, tar)) {
      lo = hi;
      break;
    }
    if (ls == 0 && rs == 0) {
      break;
    }
  }
  return {lo, hi};
}

} // namespace noya

#ifndef NOYA_RATIONAL_APPROXIMATION_HPP
#define NOYA_RATIONAL_APPROXIMATION_HPP 1

/// @complexity Time: O(log N log max(x,y)).  Space: O(1).

#include <cassert>
#include <cstdint>
#include <utility>

namespace noya {

struct bounded_rational {
  std::int64_t num = 0;
  std::int64_t den = 1;
};

namespace rational_approximation_detail {

inline bool equal(bounded_rational a, bounded_rational b) {
  return static_cast<__int128>(a.num) * b.den ==
         static_cast<__int128>(a.den) * b.num;
}

inline bool less_equal(bounded_rational a, bounded_rational b) {
  return static_cast<__int128>(a.num) * b.den <=
         static_cast<__int128>(a.den) * b.num;
}

inline bool bounded(bounded_rational val, std::int64_t lim) {
  return val.num <= lim && val.den <= lim;
}

inline bounded_rational multiply_add(bounded_rational a, std::int64_t mul,
                                     bounded_rational b) {
  return {a.num * mul + b.num, a.den * mul + b.den};
}

template <class Predicate> std::int64_t maximum_true(Predicate pre) {
  std::int64_t acc = 0;
  std::int64_t rej = 1;
  while (pre(rej)) {
    rej *= 2;
  }
  while (rej - acc > 1) {
    std::int64_t mid = (acc + rej) / 2;
    if (pre(mid)) {
      acc = mid;
    } else {
      rej = mid;
    }
  }
  return acc;
}

} // namespace rational_approximation_detail

/// @brief Find the closest reduced fractions on both sides of `x/y` whose
/// numerators and denominators are at most `lim`.  Alternating jumps between
/// the two Stern--Brocot boundaries consumes one continued-fraction run at a
/// time; each run length is maximized by integer binary search.  Missing lower
/// and upper answers are represented by `0/1` and `1/0`.
inline std::pair<bounded_rational, bounded_rational>
rational_approximation(std::int64_t lim, std::int64_t x, std::int64_t y) {
  assert(lim >= 1 && x >= 1 && y >= 1);
  bounded_rational tar{x, y};
  bounded_rational lo{0, 1};
  bounded_rational hi{1, 0};
  while (true) {
    std::int64_t ls =
        rational_approximation_detail::maximum_true([&](std::int64_t cnt) {
          bounded_rational can =
              rational_approximation_detail::multiply_add(hi, cnt, lo);
          return rational_approximation_detail::bounded(can, lim) &&
                 rational_approximation_detail::less_equal(can, tar);
        });
    lo = rational_approximation_detail::multiply_add(hi, ls, lo);
    if (rational_approximation_detail::equal(lo, tar)) {
      hi = lo;
      break;
    }

    std::int64_t rs =
        rational_approximation_detail::maximum_true([&](std::int64_t cnt) {
          bounded_rational can =
              rational_approximation_detail::multiply_add(lo, cnt, hi);
          return rational_approximation_detail::bounded(can, lim) &&
                 rational_approximation_detail::less_equal(tar, can);
        });
    hi = rational_approximation_detail::multiply_add(lo, rs, hi);
    if (rational_approximation_detail::equal(hi, tar)) {
      lo = hi;
      break;
    }
    if (ls == 0 && rs == 0) {
      break;
    }
  }
  return {lo, hi};
}

} // namespace noya

#endif // NOYA_RATIONAL_APPROXIMATION_HPP

#include <cassert>
#include <cstdint>
#include <utility>

/// @complexity Time: O(log N log max(x,y)).  Space: O(1).

namespace noya {

struct bounded_rational {
  std::int64_t num = 0;
  std::int64_t den = 1;
};

namespace rational_approximation_detail {

inline bool equal(bounded_rational a, bounded_rational b) {
  return static_cast<__int128>(a.num) * b.den ==
         static_cast<__int128>(a.den) * b.num;
}

inline bool less_equal(bounded_rational a, bounded_rational b) {
  return static_cast<__int128>(a.num) * b.den <=
         static_cast<__int128>(a.den) * b.num;
}

inline bool bounded(bounded_rational val, std::int64_t lim) {
  return val.num <= lim && val.den <= lim;
}

inline bounded_rational multiply_add(bounded_rational a, std::int64_t mul,
                                     bounded_rational b) {
  return {a.num * mul + b.num, a.den * mul + b.den};
}

template <class Predicate> std::int64_t maximum_true(Predicate pre) {
  std::int64_t acc = 0;
  std::int64_t rej = 1;
  while (pre(rej)) {
    rej *= 2;
  }
  while (rej - acc > 1) {
    std::int64_t mid = (acc + rej) / 2;
    if (pre(mid)) {
      acc = mid;
    } else {
      rej = mid;
    }
  }
  return acc;
}

} // namespace rational_approximation_detail

/// @brief Find the closest reduced fractions on both sides of `x/y` whose
/// numerators and denominators are at most `lim`.  Alternating jumps between
/// the two Stern--Brocot boundaries consumes one continued-fraction run at a
/// time; each run length is maximized by integer binary search.  Missing lower
/// and upper answers are represented by `0/1` and `1/0`.
inline std::pair<bounded_rational, bounded_rational>
rational_approximation(std::int64_t lim, std::int64_t x, std::int64_t y) {
  assert(lim >= 1 && x >= 1 && y >= 1);
  bounded_rational tar{x, y};
  bounded_rational lo{0, 1};
  bounded_rational hi{1, 0};
  while (true) {
    std::int64_t ls =
        rational_approximation_detail::maximum_true([&](std::int64_t cnt) {
          bounded_rational can =
              rational_approximation_detail::multiply_add(hi, cnt, lo);
          return rational_approximation_detail::bounded(can, lim) &&
                 rational_approximation_detail::less_equal(can, tar);
        });
    lo = rational_approximation_detail::multiply_add(hi, ls, lo);
    if (rational_approximation_detail::equal(lo, tar)) {
      hi = lo;
      break;
    }

    std::int64_t rs =
        rational_approximation_detail::maximum_true([&](std::int64_t cnt) {
          bounded_rational can =
              rational_approximation_detail::multiply_add(lo, cnt, hi);
          return rational_approximation_detail::bounded(can, lim) &&
                 rational_approximation_detail::less_equal(tar, can);
        });
    hi = rational_approximation_detail::multiply_add(lo, rs, hi);
    if (rational_approximation_detail::equal(hi, tar)) {
      lo = hi;
      break;
    }
    if (ls == 0 && rs == 0) {
      break;
    }
  }
  return {lo, hi};
}

} // namespace noya