Rational Approximation
(math/rational_approximation.hpp)

• View this file on GitHub
• Last update: 2026-07-16 21:30:39+09:00
• Include: #include "math/rational_approximation.hpp"

Overview

rational_approximation brackets a positive rational number by the closest reduced fractions whose numerators and denominators are both bounded.

For positive maximum, numerator, and denominator, it returns:

• the largest reduced $a/b\leq numerator/denominator$ with $1\leq a,b\leq maximum$, or 0/1 if none exists;
• the smallest reduced $c/d\geq numerator/denominator$ with $1\leq c,d\leq maximum$, or 1/0 if none exists.

The algorithm walks the Stern–Brocot interval with maximal run lengths. It uses only exact unsigned 128-bit intermediate arithmetic and does not use floating point.

Interface

template <std::integral T>
struct RationalApproximationResult {
    using fraction_type = std::pair<T, T>;

    fraction_type lower;
    fraction_type upper;
};

template <std::integral T>
RationalApproximationResult<T> rational_approximation(
    T maximum,
    T numerator,
    T denominator
);

Member or function	Description	Complexity
result.lower	Pair (numerator, denominator) for the lower bound.	$O(1)$ access
result.upper	Pair (numerator, denominator) for the upper bound.	$O(1)$ access
rational_approximation(maximum, numerator, denominator)	Computes both closest bounded fractions.	$O(\log \max(maximum,numerator,denominator))$ time and $O(1)$ memory

T may be a signed or unsigned integral type other than bool, with at most 64 bits. All three arguments must be positive. Returned ordinary fractions are reduced; 0/1 and 1/0 are boundary sentinels and are not members of the positive bounded set.

Example

#include "math/rational_approximation.hpp"

#include <cassert>
#include <utility>

int main() {
    auto result = m1une::math::rational_approximation(5LL, 7LL, 10LL);

    assert(result.lower == std::pair<long long, long long>(2, 3));
    assert(result.upper == std::pair<long long, long long>(3, 4));
}

Required by

• Math All (math/all.hpp)

Verified with

• verify/math/math_algorithms.test.cpp
• verify/math/rational_approximation.test.cpp

Code

#ifndef M1UNE_MATH_RATIONAL_APPROXIMATION_HPP
#define M1UNE_MATH_RATIONAL_APPROXIMATION_HPP 1

#include <algorithm>
#include <cassert>
#include <concepts>
#include <cstdint>
#include <type_traits>
#include <utility>

namespace m1une {
namespace math {

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
struct RationalApproximationResult {
    using fraction_type = std::pair<T, T>;

    fraction_type lower;
    fraction_type upper;
};

namespace rational_approximation_detail {

using Wide = __uint128_t;

struct Fraction {
    Wide numerator;
    Wide denominator;
};

inline bool equal(
    const Fraction& fraction,
    Wide numerator,
    Wide denominator
) {
    return fraction.numerator * denominator ==
           numerator * fraction.denominator;
}

inline Wide coordinate_bound(
    Wide maximum,
    Wide offset,
    Wide step
) {
    return step == 0 ? maximum : (maximum - offset) / step;
}

inline Fraction multiply_add(
    const Fraction& first,
    Wide multiplier,
    const Fraction& second
) {
    return {
        first.numerator * multiplier + second.numerator,
        first.denominator * multiplier + second.denominator
    };
}

}  // namespace rational_approximation_detail

// Returns the closest reduced fractions below and above numerator/denominator
// whose positive numerator and denominator are both at most maximum.
template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
RationalApproximationResult<T> rational_approximation(
    T maximum,
    T numerator,
    T denominator
) {
    assert(maximum > 0);
    assert(numerator > 0);
    assert(denominator > 0);
    static_assert(sizeof(T) <= sizeof(uint64_t));

    using rational_approximation_detail::Fraction;
    using rational_approximation_detail::Wide;
    using rational_approximation_detail::coordinate_bound;
    using rational_approximation_detail::equal;
    using rational_approximation_detail::multiply_add;

    const Wide limit = static_cast<std::make_unsigned_t<T>>(maximum);
    const Wide target_numerator =
        static_cast<std::make_unsigned_t<T>>(numerator);
    const Wide target_denominator =
        static_cast<std::make_unsigned_t<T>>(denominator);
    Fraction lower{0, 1};
    Fraction upper{1, 0};

    while (true) {
        Wide lower_multiplier = limit;
        lower_multiplier = std::min(
            lower_multiplier,
            coordinate_bound(limit, lower.numerator, upper.numerator)
        );
        lower_multiplier = std::min(
            lower_multiplier,
            coordinate_bound(limit, lower.denominator, upper.denominator)
        );
        const Wide lower_slack =
            target_numerator * lower.denominator -
            lower.numerator * target_denominator;
        const Wide lower_step =
            upper.numerator * target_denominator -
            target_numerator * upper.denominator;
        assert(lower_step != 0);
        lower_multiplier = std::min(
            lower_multiplier, lower_slack / lower_step
        );
        lower = multiply_add(upper, lower_multiplier, lower);
        if (equal(lower, target_numerator, target_denominator)) {
            upper = lower;
            break;
        }

        Wide upper_multiplier = limit;
        upper_multiplier = std::min(
            upper_multiplier,
            coordinate_bound(limit, upper.numerator, lower.numerator)
        );
        upper_multiplier = std::min(
            upper_multiplier,
            coordinate_bound(limit, upper.denominator, lower.denominator)
        );
        const Wide upper_slack =
            upper.numerator * target_denominator -
            target_numerator * upper.denominator;
        const Wide upper_step =
            target_numerator * lower.denominator -
            lower.numerator * target_denominator;
        assert(upper_step != 0);
        upper_multiplier = std::min(
            upper_multiplier, upper_slack / upper_step
        );
        upper = multiply_add(lower, upper_multiplier, upper);
        if (equal(upper, target_numerator, target_denominator)) {
            lower = upper;
            break;
        }

        if (lower_multiplier == 0 && upper_multiplier == 0) break;
    }

    RationalApproximationResult<T> result;
    result.lower = {
        static_cast<T>(lower.numerator),
        static_cast<T>(lower.denominator)
    };
    result.upper = {
        static_cast<T>(upper.numerator),
        static_cast<T>(upper.denominator)
    };
    return result;
}

}  // namespace math
}  // namespace m1une

#endif  // M1UNE_MATH_RATIONAL_APPROXIMATION_HPP

#line 1 "math/rational_approximation.hpp"



#include <algorithm>
#include <cassert>
#include <concepts>
#include <cstdint>
#include <type_traits>
#include <utility>

namespace m1une {
namespace math {

template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
struct RationalApproximationResult {
    using fraction_type = std::pair<T, T>;

    fraction_type lower;
    fraction_type upper;
};

namespace rational_approximation_detail {

using Wide = __uint128_t;

struct Fraction {
    Wide numerator;
    Wide denominator;
};

inline bool equal(
    const Fraction& fraction,
    Wide numerator,
    Wide denominator
) {
    return fraction.numerator * denominator ==
           numerator * fraction.denominator;
}

inline Wide coordinate_bound(
    Wide maximum,
    Wide offset,
    Wide step
) {
    return step == 0 ? maximum : (maximum - offset) / step;
}

inline Fraction multiply_add(
    const Fraction& first,
    Wide multiplier,
    const Fraction& second
) {
    return {
        first.numerator * multiplier + second.numerator,
        first.denominator * multiplier + second.denominator
    };
}

}  // namespace rational_approximation_detail

// Returns the closest reduced fractions below and above numerator/denominator
// whose positive numerator and denominator are both at most maximum.
template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
RationalApproximationResult<T> rational_approximation(
    T maximum,
    T numerator,
    T denominator
) {
    assert(maximum > 0);
    assert(numerator > 0);
    assert(denominator > 0);
    static_assert(sizeof(T) <= sizeof(uint64_t));

    using rational_approximation_detail::Fraction;
    using rational_approximation_detail::Wide;
    using rational_approximation_detail::coordinate_bound;
    using rational_approximation_detail::equal;
    using rational_approximation_detail::multiply_add;

    const Wide limit = static_cast<std::make_unsigned_t<T>>(maximum);
    const Wide target_numerator =
        static_cast<std::make_unsigned_t<T>>(numerator);
    const Wide target_denominator =
        static_cast<std::make_unsigned_t<T>>(denominator);
    Fraction lower{0, 1};
    Fraction upper{1, 0};

    while (true) {
        Wide lower_multiplier = limit;
        lower_multiplier = std::min(
            lower_multiplier,
            coordinate_bound(limit, lower.numerator, upper.numerator)
        );
        lower_multiplier = std::min(
            lower_multiplier,
            coordinate_bound(limit, lower.denominator, upper.denominator)
        );
        const Wide lower_slack =
            target_numerator * lower.denominator -
            lower.numerator * target_denominator;
        const Wide lower_step =
            upper.numerator * target_denominator -
            target_numerator * upper.denominator;
        assert(lower_step != 0);
        lower_multiplier = std::min(
            lower_multiplier, lower_slack / lower_step
        );
        lower = multiply_add(upper, lower_multiplier, lower);
        if (equal(lower, target_numerator, target_denominator)) {
            upper = lower;
            break;
        }

        Wide upper_multiplier = limit;
        upper_multiplier = std::min(
            upper_multiplier,
            coordinate_bound(limit, upper.numerator, lower.numerator)
        );
        upper_multiplier = std::min(
            upper_multiplier,
            coordinate_bound(limit, upper.denominator, lower.denominator)
        );
        const Wide upper_slack =
            upper.numerator * target_denominator -
            target_numerator * upper.denominator;
        const Wide upper_step =
            target_numerator * lower.denominator -
            lower.numerator * target_denominator;
        assert(upper_step != 0);
        upper_multiplier = std::min(
            upper_multiplier, upper_slack / upper_step
        );
        upper = multiply_add(lower, upper_multiplier, upper);
        if (equal(upper, target_numerator, target_denominator)) {
            lower = upper;
            break;
        }

        if (lower_multiplier == 0 && upper_multiplier == 0) break;
    }

    RationalApproximationResult<T> result;
    result.lower = {
        static_cast<T>(lower.numerator),
        static_cast<T>(lower.denominator)
    };
    result.upper = {
        static_cast<T>(upper.numerator),
        static_cast<T>(upper.denominator)
    };
    return result;
}

}  // namespace math
}  // namespace m1une

Back to top page