rational_approximation.hpp¶
Find the closest reduced fractions on both sides of x/y whose numerators and denominators are at most limit. 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.
Verified by rational_approximation.
\[
\displaystyle \frac{L}{R}<\frac{x}{y}<\frac{U}{V}
\]
Implementation¶
#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 numerator = 0;
std::int64_t denominator = 1;
};
namespace rational_approximation_detail {
inline bool equal(bounded_rational first, bounded_rational second) {
return static_cast<__int128>(first.numerator) * second.denominator ==
static_cast<__int128>(first.denominator) * second.numerator;
}
inline bool less_equal(bounded_rational first, bounded_rational second) {
return static_cast<__int128>(first.numerator) * second.denominator <=
static_cast<__int128>(first.denominator) * second.numerator;
}
inline bool bounded(bounded_rational value, std::int64_t limit) {
return value.numerator <= limit && value.denominator <= limit;
}
inline bounded_rational multiply_add(bounded_rational first,
std::int64_t multiplier,
bounded_rational second) {
return {first.numerator * multiplier + second.numerator,
first.denominator * multiplier + second.denominator};
}
template <class Predicate> std::int64_t maximum_true(Predicate predicate) {
std::int64_t accepted = 0;
std::int64_t rejected = 1;
while (predicate(rejected)) {
rejected *= 2;
}
while (rejected - accepted > 1) {
std::int64_t middle = (accepted + rejected) / 2;
if (predicate(middle)) {
accepted = middle;
} else {
rejected = middle;
}
}
return accepted;
}
} // namespace rational_approximation_detail
/// @brief Find the closest reduced fractions on both sides of `x/y` whose
/// numerators and denominators are at most `limit`. 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 limit, std::int64_t x, std::int64_t y) {
assert(limit >= 1 && x >= 1 && y >= 1);
bounded_rational target{x, y};
bounded_rational lower{0, 1};
bounded_rational upper{1, 0};
while (true) {
std::int64_t lower_steps =
rational_approximation_detail::maximum_true([&](std::int64_t count) {
bounded_rational candidate =
rational_approximation_detail::multiply_add(upper, count,
lower);
return rational_approximation_detail::bounded(candidate, limit) &&
rational_approximation_detail::less_equal(candidate, target);
});
lower = rational_approximation_detail::multiply_add(upper, lower_steps,
lower);
if (rational_approximation_detail::equal(lower, target)) {
upper = lower;
break;
}
std::int64_t upper_steps =
rational_approximation_detail::maximum_true([&](std::int64_t count) {
bounded_rational candidate =
rational_approximation_detail::multiply_add(lower, count,
upper);
return rational_approximation_detail::bounded(candidate, limit) &&
rational_approximation_detail::less_equal(target, candidate);
});
upper = rational_approximation_detail::multiply_add(lower, upper_steps,
upper);
if (rational_approximation_detail::equal(upper, target)) {
lower = upper;
break;
}
if (lower_steps == 0 && upper_steps == 0) {
break;
}
}
return {lower, upper};
}
} // 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 numerator = 0;
std::int64_t denominator = 1;
};
namespace rational_approximation_detail {
inline bool equal(bounded_rational first, bounded_rational second) {
return static_cast<__int128>(first.numerator) * second.denominator ==
static_cast<__int128>(first.denominator) * second.numerator;
}
inline bool less_equal(bounded_rational first, bounded_rational second) {
return static_cast<__int128>(first.numerator) * second.denominator <=
static_cast<__int128>(first.denominator) * second.numerator;
}
inline bool bounded(bounded_rational value, std::int64_t limit) {
return value.numerator <= limit && value.denominator <= limit;
}
inline bounded_rational multiply_add(bounded_rational first,
std::int64_t multiplier,
bounded_rational second) {
return {first.numerator * multiplier + second.numerator,
first.denominator * multiplier + second.denominator};
}
template <class Predicate> std::int64_t maximum_true(Predicate predicate) {
std::int64_t accepted = 0;
std::int64_t rejected = 1;
while (predicate(rejected)) {
rejected *= 2;
}
while (rejected - accepted > 1) {
std::int64_t middle = (accepted + rejected) / 2;
if (predicate(middle)) {
accepted = middle;
} else {
rejected = middle;
}
}
return accepted;
}
} // namespace rational_approximation_detail
/// @brief Find the closest reduced fractions on both sides of `x/y` whose
/// numerators and denominators are at most `limit`. 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 limit, std::int64_t x, std::int64_t y) {
assert(limit >= 1 && x >= 1 && y >= 1);
bounded_rational target{x, y};
bounded_rational lower{0, 1};
bounded_rational upper{1, 0};
while (true) {
std::int64_t lower_steps =
rational_approximation_detail::maximum_true([&](std::int64_t count) {
bounded_rational candidate =
rational_approximation_detail::multiply_add(upper, count,
lower);
return rational_approximation_detail::bounded(candidate, limit) &&
rational_approximation_detail::less_equal(candidate, target);
});
lower = rational_approximation_detail::multiply_add(upper, lower_steps,
lower);
if (rational_approximation_detail::equal(lower, target)) {
upper = lower;
break;
}
std::int64_t upper_steps =
rational_approximation_detail::maximum_true([&](std::int64_t count) {
bounded_rational candidate =
rational_approximation_detail::multiply_add(lower, count,
upper);
return rational_approximation_detail::bounded(candidate, limit) &&
rational_approximation_detail::less_equal(target, candidate);
});
upper = rational_approximation_detail::multiply_add(lower, upper_steps,
upper);
if (rational_approximation_detail::equal(upper, target)) {
lower = upper;
break;
}
if (lower_steps == 0 && upper_steps == 0) {
break;
}
}
return {lower, upper};
}
} // namespace noya