rational.hpp

Go to the documentation of this file.

#pragma once
 
#include "core/bigint.hpp"

 
#include <cmath>
#include <iostream>
#include <stdexcept>
#include <numeric>
#include <type_traits>
#include <compare>
#include <limits>
#include <functional>
#include <cstdint>
#include <utility>
#include <concepts>
#include <cassert>
#include <cstddef>
 
 
 
namespace pgl {
template <class T>
struct is_Rational : std::false_type {};
 
template <class T>
struct is_Rational<Rational<T>> : std::true_type {};
 
template <class T>
inline constexpr bool is_Rational_v = is_Rational<T>::value;
 
template <class T>
concept RationalConcept = is_Rational_v<T>;

template <class T>
struct rational_int {
    using type = T;
};
 
template <class Int>
struct rational_int<Rational<Int>> {
    using type = Int;
};
 
template <class T>
using rational_int_t = typename rational_int<T>::type;

template <class T>
using grid_number_t = std::conditional_t<std::signed_integral<rational_int_t<T>>,
                                         rational_int_t<T>, std::int64_t>;
 
namespace detail {
inline pgl::BigInt pow2(int k) {
    pgl::BigInt result(1), base(2);
    while (k > 0) {
        if (k & 1) result *= base;
        base *= base;
        k >>= 1;
    }
    return result;
}
} // namespace detail
 

template <class Int = int64_t>
class Rational {
private:
    Int num;                 
    Int den;                 
    bool normalized_ = true; 

    static constexpr std::size_t heapReductionLimbs = 8;

    static constexpr bool reductionUrgent(const Int& n, const Int& d) {
        if constexpr (requires { n.fitsInt64(); }) {
            // BigInt, and the answer turns on what the gcd itself will cost.
            //
            // While both parts are inline (a single int128), BigInt divides on
            // its native fast path, so the gcd is a handful of machine
            // divisions. Reduce at half that width, exactly as the fixed-width
            // rule below does, and the impending multiply or cross-addition
            // stays inline too.
            //
            // Once a part has spilled onto the heap that bargain is off: every
            // step of the gcd becomes a multi-limb long division, which costs
            // far more than carrying the extra limbs through the one multiply
            // the reduction was meant to shrink — and @ref operandParts throws
            // the reduced form away with its locals, so the same value pays
            // again at every use. Past that point reduce only to stop unbounded
            // growth, not to save width.
            if (n.fitsInt128() && d.fitsInt128()) {
                return !n.fitsInt64() || !d.fitsInt64();
            }
            return !n.fitsLimbs(heapReductionLimbs) || !d.fitsLimbs(heapReductionLimbs);
        } else {
            constexpr int half = pgl::detail::numeric_limits<Int>::digits / 2 - 1;
            const Int limit = Int(1) << half;
            return pgl::detail::abs(n) >= limit || d >= limit;
        }
    }

    static constexpr bool comparisonNeedsReduction(const Int& n, const Int& d) {
        if constexpr (requires { n.fitsInt64(); }) {
            return false;
        } else {
            return reductionUrgent(n, d);
        }
    }

    constexpr bool storedInteger() const {
        // An arbitrary-precision Int answers this off its inline store; building
        // a value to compare against would cost more than the shortcut saves.
        if constexpr (requires { den.isOne(); }) {
            return den.isOne();
        } else {
            return den == Int(1);
        }
    }

    constexpr void operandParts(Int& n, Int& d) const {
        n = num;
        d = den;
        if (!normalized_ && reductionUrgent(n, d)) {
            const Int g = pgl::detail::gcd(pgl::detail::abs(n), d);
            n /= g;
            d /= g;
        }
    }
 
public:
    constexpr Rational() : num(0), den(1) {}

    constexpr Rational(RationalConcept auto n) : num(0), den(1) {
        n.simplify();
        num = Int(n.numerator());
        den = Int(n.denominator());
    }

    constexpr Rational(pgl::detail::extended_integral auto n) : num(n), den(1) {}

    template <std::same_as<Int> T>
        requires(!pgl::detail::extended_integral<Int> && !std::floating_point<Int>
                 && !RationalConcept<Int>)
    constexpr Rational(T n) : num(std::move(n)), den(1) {}

    constexpr Rational(Int n, Int d, bool normalized = false)
        : num(n), den(d), normalized_(normalized) {
        if (den < 0) {
            num = -num;
            den = -den;
        }
 
        assert(den > 0);
        // The reduction is deferred: it happens lazily at a read, a comparison,
        // or when an arithmetic step would otherwise risk overflow. Integers are
        // already in lowest terms, so flag them normalized to skip that work.
        if (den == 1) normalized_ = true;
    }

    template<std::floating_point Float>
    explicit constexpr Rational(Float f,
                       int digits = pgl::detail::numeric_limits<Int>::digits > 0
                           ? std::min(pgl::detail::numeric_limits<Float>::digits,
                                      pgl::detail::numeric_limits<Int>::digits / 2 - 4)
                           : pgl::detail::numeric_limits<Float>::digits) {
        if (!std::isfinite(f)) {
            throw std::domain_error("pgl::Rational: cannot construct from a non-finite floating-point value");
        }
        if (f == 0) {
            num = 0;
            den = 1;
            return;
        }
 
        // |f| lies in [2^(exponent - 1), 2^exponent), so exponent is how many
        // of the requested bits the integer part takes for itself.
        int exponent = 0;
        std::frexp(f, &exponent);
        const int fractionBits = exponent > 0 ? digits - exponent : digits;
        if (fractionBits <= 0) {
            // No fraction bit fits beside the integer part: the value is that
            // part, exactly, whenever the storage type holds it.
            if constexpr (pgl::detail::numeric_limits<Int>::is_bounded) {
                if (exponent > pgl::detail::numeric_limits<Int>::digits) {
                    throw std::overflow_error("pgl::Rational: floating-point value does not fit the integer type");
                }
            }
            num = Int(std::trunc(f));
            den = Int(1);
            normalized_ = true;
            return;
        }
 
        if constexpr (requires(Int x, Float g) { x * g; }) {
            bool negative = false;
            if (f < 0) {
                negative = true;
                f = -f;
            }
            den = Int(1) << fractionBits;
            // Wrap in Int(): for the Boost int128 fallback `den * f` yields a
            // double (see the double-interop shim in numeric.hpp), so convert it
            // back explicitly. For native __int128 / built-in ints this is the
            // same truncating conversion that the implicit assignment did.
            num = Int(negative ? -den * f : den * f);
 
            normalized_ = false;
            simplify();
        } else {
            // Int has no floating-point multiply (e.g. pgl::BigInt). A float's
            // significand always fits in a 128-bit integer, so build the value
            // with the requested number of fractional bits in Rational<int128>
            // and widen it. This makes float -> Rational<BigInt> work, which the
            // int128 -> BigInt type promotion relies on for mixed predicates.
            const Rational<pgl::int128> r(f, digits);
            num = Int(r.numerator());
            den = Int(r.denominator());
            normalized_ = true;
        }
    }

    constexpr Int numerator() const noexcept {
        if (normalized_) return num;
        return num / pgl::detail::gcd(pgl::detail::abs(num), den);
    }

    constexpr Int denominator() const noexcept {
        if (normalized_) return den;
        return den / pgl::detail::gcd(pgl::detail::abs(num), den);
    }

    constexpr void simplify() {
        if (normalized_) {
            return;
        }
        normalized_ = true;
        assert(den != 0);
 
        if (num == 0) {
            den = 1;
            return;
        }
 
        if (den == 1 || num == 1 || num == -1) {
            return;   // a unit numerator or denominator is already in lowest terms
        }
 
        Int g = pgl::detail::gcd(pgl::detail::abs(num), den);
        num /= g;
        den /= g;
    }

    [[nodiscard]] constexpr Rational simplified() const {
        if (normalized_) {
            return *this;
        }
        if (num == 0) {
            return Rational(Int(0), Int(1), true);
        }
        if (den == 1 || num == 1 || num == -1) {
            return Rational(num, den, true);
        }
        const Int g = pgl::detail::gcd(pgl::detail::abs(num), den);
        return Rational(num / g, den / g, true);
    }

    constexpr void simplifyIfLarge() {
        if (!normalized_ && reductionUrgent(num, den)) {
            simplify();
        }
    }

    [[nodiscard]] constexpr Rational simplifiedIfLarge() const {
        if (!normalized_ && reductionUrgent(num, den)) {
            return simplified();
        }
        return *this;
    }

    [[nodiscard]] constexpr bool isInteger() const {
        return normalized_ ? den == Int(1) : num % den == Int(0);
    }

    explicit constexpr operator float() const {
        if (den == Int(1)) {
            return static_cast<float>(num);
        }
        return static_cast<float>(num) / static_cast<float>(den);
    }

    explicit constexpr operator double() const {
        if (den == Int(1)) {
            return static_cast<double>(num);
        }
        return static_cast<double>(num) / static_cast<double>(den);
    }

    explicit constexpr operator long double() const {
        if (den == Int(1)) {
            return static_cast<long double>(num);
        }
        return static_cast<long double>(num) / static_cast<long double>(den);
    }

    explicit constexpr operator int() const {
        if (den == Int(1)) {
            return static_cast<int>(num);
        }
        return static_cast<int>(num / den);
    }

    explicit constexpr operator int64_t() const {
        if (den == Int(1)) {
            return static_cast<int64_t>(num);
        }
        return static_cast<int64_t>(num / den);
    }

    explicit constexpr operator Int() const
        requires(!std::same_as<Int, int> && !std::same_as<Int, int64_t>) {
        if (den == Int(1)) {
            return num;
        }
        return num / den;
    }

    template <class OtherInt>
    explicit constexpr operator Rational<OtherInt>() const {
        return Rational<OtherInt>(num, den, normalized_);
    }

    template <std::floating_point Float = double>
    constexpr Float lowerBound() const {
        using Double = pgl::detail::promoted_number_t<Float>;
 
        Double result = static_cast<Double>(num) / static_cast<Double>(den);
        Float flt_result = static_cast<Float>(result);
 
#ifndef NDEBUG
        int i = 0;
#endif
        while (static_cast<Double>(flt_result) > result) {
            flt_result = std::nextafter(flt_result, -pgl::detail::numeric_limits<Float>::infinity());
            assert(i++ < 10); // Normally one iteration should be enough
        }
        // The quotient is itself a rounding of the exact value, so when the
        // narrowing landed exactly on it nothing above tells which side of the
        // value the candidate is on: settle that exactly. One step then
        // suffices, since a float strictly below the quotient is below every
        // value the quotient can have rounded from.
        if (static_cast<Double>(flt_result) == result && (*this <=> flt_result) < 0) {
            flt_result = std::nextafter(flt_result, -pgl::detail::numeric_limits<Float>::infinity());
        }
 
        return flt_result;
    }

    template <std::floating_point Float = double>
    constexpr Float upperBound() const {
        using Double = pgl::detail::promoted_number_t<Float>;
 
        Double result = static_cast<Double>(num) / static_cast<Double>(den);
        Float flt_result = static_cast<Float>(result);
 
#ifndef NDEBUG
        int i = 0;
#endif
        while (static_cast<Double>(flt_result) < result) {
            flt_result = std::nextafter(flt_result, pgl::detail::numeric_limits<Float>::infinity());
            assert(i++ < 10); // Normally one iteration should be enough
        }
        // See lowerBound: a candidate equal to the rounded quotient may still
        // be below the exact value, and one exact step settles it.
        if (static_cast<Double>(flt_result) == result && (*this <=> flt_result) > 0) {
            flt_result = std::nextafter(flt_result, pgl::detail::numeric_limits<Float>::infinity());
        }
 
        return flt_result;
    }
 
    // Arithmetic results are left unreduced (the gcd is deferred); each operand
    // is only reduced first when its magnitude would otherwise risk overflow.
    // An integer result (denominator 1) is trivially in lowest terms.

    constexpr bool safeRaw() const {
        return normalized_ || !reductionUrgent(num, den);
    }
 
    constexpr Rational operator+(const Rational& r) const {
        // Two integers add and multiply as integers, over the 1 they already
        // share: the general form below would cross-multiply by it three times.
        if (storedInteger() && r.storedInteger()) {
            return Rational(num + r.num, den, true);
        }
        if (safeRaw() && r.safeRaw()) {
            const Int rd = den * r.den;
            return Rational(num * r.den + r.num * den, rd, rd == 1);
        }
        Int an, ad, bn, bd;
        operandParts(an, ad);
        r.operandParts(bn, bd);
        const Int rd = ad * bd;
        return Rational(an * bd + bn * ad, rd, rd == 1);
    }
 
    constexpr Rational operator-(const Rational& r) const {
        // Two integers add and multiply as integers, over the 1 they already
        // share: the general form below would cross-multiply by it three times.
        if (storedInteger() && r.storedInteger()) {
            return Rational(num - r.num, den, true);
        }
        if (safeRaw() && r.safeRaw()) {
            const Int rd = den * r.den;
            return Rational(num * r.den - r.num * den, rd, rd == 1);
        }
        Int an, ad, bn, bd;
        operandParts(an, ad);
        r.operandParts(bn, bd);
        const Int rd = ad * bd;
        return Rational(an * bd - bn * ad, rd, rd == 1);
    }
 
    constexpr Rational<Int> operator*(const Rational<Int>& r) const {
        if (storedInteger() && r.storedInteger()) {
            return Rational(num * r.num, den, true);
        }
        if (safeRaw() && r.safeRaw()) {
            const Int rd = den * r.den;
            return Rational(num * r.num, rd, rd == 1);
        }
        Int an, ad, bn, bd;
        operandParts(an, ad);
        r.operandParts(bn, bd);
        const Int rd = ad * bd;
        return Rational(an * bn, rd, rd == 1);
    }
 
    constexpr Rational reciprocal() const {
        assert(num != 0);
        // Swapping numerator and denominator preserves gcd(|num|, den), so the
        // normalization state carries over unchanged.
        return num < 0 ? Rational(-den, -num, normalized_)
                       : Rational(den, num, normalized_);
    }
 
    constexpr Rational operator/(const Rational& r) const {
        return *this * r.reciprocal();
    }
 
    constexpr Rational operator-() const {
        // Negation preserves gcd(|num|, den), so the normalization state holds.
        return Rational(-num, den, normalized_);
    }
 
    // Compound
    constexpr Rational& operator+=(const Rational& r) { return *this = *this + r; }
    constexpr Rational& operator-=(const Rational& r) { return *this = *this - r; }
    constexpr Rational& operator*=(const Rational& r) { return *this = *this * r; }
    constexpr Rational& operator/=(const Rational& r) { return *this = *this / r; }
 
    // Mixed integer operators. Templating on the (deduced) argument type makes
    // an integer operand bind without a conversion, so these win cleanly over
    // operator OP(const Rational&) instead of tying with it. The tie — and the
    // resulting ambiguity — only arose when `int -> Int` is itself a
    // user-defined conversion, i.e. for a class-type Int such as BigInt.
    template <class I>
        requires (pgl::detail::extended_integral<I> || std::same_as<I, Int>)
    constexpr Rational operator+(const I& x) const { return *this + Rational(Int(x), true); }
 
    template <class I>
        requires (pgl::detail::extended_integral<I> || std::same_as<I, Int>)
    constexpr Rational operator-(const I& x) const { return *this - Rational(Int(x), true); }
 
    template <class I>
        requires (pgl::detail::extended_integral<I> || std::same_as<I, Int>)
    constexpr Rational operator*(const I& x) const {
        Int n, d;
        operandParts(n, d);
        return Rational(n * Int(x), d);
    }
 
    template <class I>
        requires (pgl::detail::extended_integral<I> || std::same_as<I, Int>)
    constexpr Rational operator/(const I& x) const {
        Int n, d;
        operandParts(n, d);
        return Rational(n, d * Int(x));
    }
 
    friend constexpr Rational operator+(Int x, const Rational& r) { return r + x; }
    friend constexpr Rational operator-(Int x, const Rational& r) { return r + (-x); }
    friend constexpr Rational operator*(Int x, const Rational& r) { return r * x; }
    friend constexpr Rational operator/(Int x, const Rational& r) {
        return x == 1 ? r.reciprocal() : Rational(x) / r;
    }

    template <class A, class B>
    static constexpr std::strong_ordering compareValues(const A& a, const B& b) {
        if constexpr (requires {
                          { a <=> b } -> std::same_as<std::strong_ordering>;
                      }) {
            return a <=> b;
        } else {
            if (a < b)
                return std::strong_ordering::less;
            if (b < a)
                return std::strong_ordering::greater;
            return std::strong_ordering::equal;
        }
    }

    constexpr std::strong_ordering operator<=>(const Rational& r) const {
        using Wide = pgl::detail::promoted_number_t<Int>;
        // Deferred fractions compare exactly by cross-multiplying in the wider
        // type (both denominators are positive, so the sign is preserved, and
        // compareValues yields the equal case). Reduce first only when an operand
        // is large enough that the widened product could itself overflow.
        // Two integers compare as integers: both cross products below would be
        // a multiplication by one, and over a BigInt that is a real one.
        if (storedInteger() && r.storedInteger()) {
            return compareValues(num, r.num);
        }
        if ((!normalized_ && comparisonNeedsReduction(num, den)) ||
            (!r.normalized_ && comparisonNeedsReduction(r.num, r.den))) {
            Int an, ad, bn, bd;
            operandParts(an, ad);
            r.operandParts(bn, bd);
            return compareValues(static_cast<Wide>(an) * bd, static_cast<Wide>(bn) * ad);
        }
        return compareValues(static_cast<Wide>(num) * r.den,
                             static_cast<Wide>(r.num) * den);
    }

    template <class U>
        requires (!std::same_as<U, Int>)
    constexpr std::strong_ordering operator<=>(const Rational<U>& r) const {
        using Common = std::common_type_t<Int, U>;
        using Wide = pgl::detail::promoted_number_t<Common>;
        // Both parts of the other value are wanted, and a deferred fraction
        // reduces itself anew on each read, so reduce it once and read twice.
        const Rational<U> other = r.simplified();
        const Wide lhs = static_cast<Wide>(num) * static_cast<Wide>(other.denominator());
        const Wide rhs = static_cast<Wide>(other.numerator()) * static_cast<Wide>(den);
        return compareValues(lhs, rhs);
    }

    constexpr bool operator==(const Rational& r) const {
        using Wide = pgl::detail::promoted_number_t<Int>;
        if (storedInteger() && r.storedInteger()) {
            return num == r.num;
        }
        // Deferred fractions are equal iff their cross products match; reduce
        // first only when an operand could overflow the widened product.
        if ((!normalized_ && comparisonNeedsReduction(num, den)) ||
            (!r.normalized_ && comparisonNeedsReduction(r.num, r.den))) {
            Int an, ad, bn, bd;
            operandParts(an, ad);
            r.operandParts(bn, bd);
            return static_cast<Wide>(an) * bd == static_cast<Wide>(bn) * ad;
        }
        return static_cast<Wide>(num) * r.den == static_cast<Wide>(r.num) * den;
    }
    constexpr bool operator!=(const Rational& r) const {
        return !(*this == r);
    }

    template <class U>
        requires (!std::same_as<U, Int>)
    constexpr bool operator==(const Rational<U>& r) const {
        const Rational mine = simplified();
        const Rational<U> other = r.simplified();
        return mine.num == other.numerator() && mine.den == other.denominator();
    }
 
    // --- comparison against floating point ---
    //
    // Exact for every finite value and exponent. The float f is decomposed by
    // frexp into an integer significand times a power of two (f == sig * 2^E,
    // with |sig| < 2^digits so it fits an int128), and num/den <=> sig*2^E is
    // settled by cross-multiplying through the positive denominator in pgl::BigInt
    // so the binary shift can never overflow. No precision is lost regardless of
    // the magnitudes involved. Being members, these also cover the reversed forms
    // (`f == r`, `f < r`, ...) via C++20 rewriting. The float -> Rational
    // constructor is explicit and lossy, so an implicit conversion never silently
    // approximates f here; this overload is the only float comparison path.
    //
    // Not constexpr: pgl::BigInt arithmetic is not constexpr-capable.

    template <std::floating_point Float>
    std::partial_ordering operator<=>(Float f) const {
        if (std::isnan(f))
            return std::partial_ordering::unordered;
        if (std::isinf(f))
            return f > 0 ? std::partial_ordering::less
                         : std::partial_ordering::greater;
 
        // 0 and ±1 are settled in the integer domain, skipping the BigInt
        // decomposition below entirely: den > 0, so the answer is the numerator
        // against 0 or against ±den. (-0.0 compares equal to 0 and is handled
        // here as the zero it is.)
        if (f == 0) return compareValues(num, Int(0));
        if (f == 1) return compareValues(num, den);
        if (f == -1) return compareValues(num, -den);
 
        // f == sig * 2^E exactly.
        int exponent;
        const Float frac = std::frexp(f, &exponent);
        const int digits = pgl::detail::numeric_limits<Float>::digits;
        const pgl::int128 sig = static_cast<pgl::int128>(std::ldexp(frac, digits));
        const int E = exponent - digits;
 
        // Compare num/den against sig*2^E by multiplying through by den > 0:
        //   sign(num/den - f) == sign(num - den*sig*2^E).
        // The factor of 2^|E| is moved to whichever side keeps the exponent
        // non-negative, and all of it is done in BigInt so nothing overflows.
        pgl::BigInt lhs(num);
        pgl::BigInt rhs = pgl::BigInt(den) * pgl::BigInt(sig);
        if (E >= 0)
            rhs *= detail::pow2(E);
        else
            lhs *= detail::pow2(-E);
        return lhs <=> rhs;
    }

    template <std::floating_point Float>
    bool operator==(Float f) const {
        return std::isfinite(f) &&
               (*this <=> f) == std::partial_ordering::equivalent;
    }
 
    // --- comparison against integers ---
    //
    // An integer converts implicitly to Rational, so `r == i` and `i < r` already
    // work through the same-type operators. `i == r` does not, and neither does
    // `i != r`: C++20 builds the reversed candidate for `==` out of the *declared*
    // operator== functions, and the same-type one only matches after converting
    // the integer, which is not something a reversed candidate is formed through.
    // The relational forms escape this because they are rewritten via operator<=>,
    // which is why `i < r` compiles while `i == r` does not — an asymmetry that
    // surfaces wherever a predicate compares two coordinates of different types,
    // such as Rectangle<Point<Rational>>::boundaryContains(Point<int>).
    //
    // Declaring the integer comparison as its own member template fixes both
    // directions at once, exactly as the floating-point overloads above already
    // do: the template matches the integer directly, so the reversed candidate is
    // formed without any user-defined conversion. `Int` is admitted alongside the
    // built-in integers to cover a storage type that is not one of them, BigInt
    // being the case that matters here.

    template <class I>
    static constexpr bool isNegativeOne(const I& n) {
        if constexpr (std::is_unsigned_v<I>) {
            return false;
        } else {
            return n == -1;
        }
    }

    template <class I>
        requires(pgl::detail::extended_integral<I> || std::same_as<I, Int>)
    constexpr std::strong_ordering operator<=>(const I& n) const {
        using Wide = pgl::detail::promoted_number_t<Int>;
 
        // The three constant tests come first: the integer is nearly always a
        // literal, and then they fold away at compile time along with every
        // branch not taken.
        if (n == 0) return compareValues(num, Int(0));
        if (n == 1) return compareValues(num, den);
        if (isNegativeOne(n)) return compareValues(num, -den);
        if (den == 1) return compareValues(num, static_cast<Wide>(n));
 
        // The general case multiplies through the positive denominator. Only the
        // integer is cast to the wide type: the denominator widens on its own
        // inside the product, and naming Wide for it too would copy an
        // arbitrary-precision denominator for nothing. A deferred fraction is
        // reduced first only when that product could overflow, exactly as the
        // Rational-vs-Rational comparison does.
        if (!normalized_ && comparisonNeedsReduction(num, den)) {
            Int an, ad;
            operandParts(an, ad);
            return compareValues(an, ad * static_cast<Wide>(n));
        }
        return compareValues(num, den * static_cast<Wide>(n));
    }

    template <class I>
        requires(pgl::detail::extended_integral<I> || std::same_as<I, Int>)
    constexpr bool operator==(const I& n) const {
        return (*this <=> n) == 0;
    }
 
    constexpr Rational& operator++() { return *this += 1; }
    constexpr Rational operator++(int) { Rational tmp = *this; ++(*this); return tmp; }
 
    constexpr Rational& operator--() { return *this -= 1; }
    constexpr Rational operator--(int) { Rational tmp = *this; --(*this); return tmp; }

    friend std::ostream& operator<<(std::ostream& os, const Rational& r) {
        // Reduced once, rather than once per part read.
        const Rational reduced = r.simplified();
        if (reduced.den == 1)
            return os << reduced.num;
        else
            return os << reduced.num << "/" << reduced.den;
    }

    friend std::istream& operator>>(std::istream& is, Rational& r) {
        Int n, d = 1;
        char sep = 0;
 
        if (!(is >> n)) return is;
 
        // A bare integer (no "/den" suffix) is a valid, successful parse, but
        // peek() at end-of-stream sets failbit even though nothing is wrong;
        // clear it so a plain integer at the end of the stream doesn't fail.
        if (!is.eof()) {
            if (is.peek() == '/') {
                is >> sep >> d;
            } else if (is.fail()) {
                is.clear(is.rdstate() & ~std::ios::failbit);
            }
        }
 
        r = Rational(n, d);
        return is;
    }
};

template <class T>
    requires(!pgl::detail::extended_integral<T> && !std::floating_point<T>
             && !RationalConcept<T>)
Rational(T) -> Rational<T>;
 
 
 
// Primary template: no valid type by default
template <int Bits>
struct select_int_ge
{
    using type = void;
};
 
// Specializations for supported sizes
template <> struct select_int_ge<8>  { using type = std::int8_t;  };
template <> struct select_int_ge<16> { using type = std::int16_t; };
template <> struct select_int_ge<32> { using type = std::int32_t; };
template <> struct select_int_ge<64> { using type = std::int64_t; };
template <> struct select_int_ge<128> { using type = pgl::int128; };
 
// Helper to round required bits up to the next supported size
constexpr int round_up_bits(int bits) {
    if (bits <= 8)   return 8;
    if (bits <= 16)  return 16;
    if (bits <= 32)  return 32;
    if (bits <= 64)  return 64;
    return 128;
}
 
template <typename T>
struct is_int128 : std::false_type {};
 
template <> struct is_int128<pgl::int128> : std::true_type {};
#ifdef __SIZEOF_INT128__
template <> struct is_int128<__uint128_t> : std::true_type {};
#endif
 
template <typename T>
concept NumericType =
    std::is_integral_v<T> ||
    std::is_floating_point_v<T> ||
    is_int128<T>::value;
 
template <NumericType T>
struct to_integer_with_digits {
private:
    static constexpr int bits = pgl::detail::numeric_limits<T>::digits;
    static constexpr int rounded = round_up_bits(bits);
 
public:
    using type = typename select_int_ge<rounded>::type;
};
 
template <typename T>
using to_integer_with_digits_t = typename to_integer_with_digits<T>::type;

using ERational = Rational<BigInt>;

template <class Number>
struct DivisionResult {
    using NumberType = std::remove_cvref_t<Number>;
    using type = std::conditional_t<
        std::floating_point<NumberType> || RationalConcept<NumberType>,
        NumberType,
        ERational>;
};

template <class Number>
using division_result_t = typename DivisionResult<Number>::type;
 
}// pgl
 
template<class U, class V>
struct std::common_type<pgl::Rational<U>, pgl::Rational<V>>{
    using type = pgl::Rational<std::common_type_t<U,V>>;
};
 
// Mixing a Rational with a floating-point type yields the floating-point type:
// the exact value is intentionally abandoned (the float→Rational conversion is
// explicit and lossy), so promotion collapses toward the inexact representation
// rather than silently approximating the float as a Rational.
template<class U, std::floating_point F>
struct std::common_type<pgl::Rational<U>, F>{
    using type = F;
};
template<std::floating_point F, class V>
struct std::common_type<F, pgl::Rational<V>>{
    using type = F;
};

template <class Int>
struct std::numeric_limits<pgl::Rational<Int>> {
    static constexpr bool is_specialized = true;
 
    static constexpr pgl::Rational<Int> min() noexcept {
        return pgl::Rational<Int>(pgl::detail::numeric_limits<Int>::min(), 1);
    }
 
    static constexpr pgl::Rational<Int> max() noexcept {
        return pgl::Rational<Int>(pgl::detail::numeric_limits<Int>::max(), 1);
    }
 
    static constexpr pgl::Rational<Int> lowest() noexcept {
        return min();
    }
 
    static constexpr int digits = pgl::detail::numeric_limits<Int>::digits;
    static constexpr int digits10 = pgl::detail::numeric_limits<Int>::digits10;
 
    static constexpr bool is_signed = pgl::detail::numeric_limits<Int>::is_signed;
    static constexpr bool is_integer = false;
    static constexpr bool is_exact = true;
 
    // Remaining members keep the specialization complete so generic trait
    // machinery (e.g. Boost.Multiprecision, used as the int128 fallback) can
    // classify Rational. With is_integer == false and max_exponent == 0 it is
    // treated as an "unknown" category, so Boost never tries to convert it.
    static constexpr int radix = 2;
    static constexpr int min_exponent = 0;
    static constexpr int min_exponent10 = 0;
    static constexpr int max_exponent = 0;
    static constexpr int max_exponent10 = 0;
    static constexpr bool is_bounded = pgl::detail::numeric_limits<Int>::is_bounded;
    static constexpr bool is_modulo = false;
    static constexpr bool is_iec559 = false;
    static constexpr bool has_infinity = false;
    static constexpr bool has_quiet_NaN = false;
    static constexpr bool has_signaling_NaN = false;
    static constexpr bool traps = false;
    static constexpr bool tinyness_before = false;
    static constexpr std::float_round_style round_style = std::round_toward_zero;
 
    static constexpr pgl::Rational<Int> epsilon() noexcept { return pgl::Rational<Int>(0, 1); }
    static constexpr pgl::Rational<Int> round_error() noexcept { return pgl::Rational<Int>(0, 1); }
    static constexpr pgl::Rational<Int> infinity() noexcept { return pgl::Rational<Int>(0, 1); }
    static constexpr pgl::Rational<Int> quiet_NaN() noexcept { return pgl::Rational<Int>(0, 1); }
    static constexpr pgl::Rational<Int> signaling_NaN() noexcept { return pgl::Rational<Int>(0, 1); }
    static constexpr pgl::Rational<Int> denorm_min() noexcept { return pgl::Rational<Int>(0, 1); }
};