// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// Rational.cpp
// Implementation of the arbitrary-precision rational number class
//
// Improvements:
//   - IEEE 754 decomposition via std::bit_cast (C++23)
//   - Fixed exponent computation bug for subnormal numbers
//   - Added denominator==0 check
//   - Use std::move in the move constructor
//   - Fixed rounding bug in toDecimal

#include <math/core/mp/Rational.hpp>
#include <math/core/mp/Int/IntOps.hpp>
#include <math/core/mp/Float.hpp>
#include <math/core/mp/Int/IntGCD.hpp>
#include <math/core/mp/Int/IntModular.hpp>
#include <bit>
#include <numeric>      // std::gcd
#include <cmath>        // std::abs, std::max
#include <math/core/mp/Int/compiler_intrinsics.hpp>

namespace {
    // Helper for the double-word fast path
    // Check whether both numerator and denominator fit in 1 limb
    inline bool isSmallRational(const sangi::Int& num, const sangi::Int& den) {
        return num.size() <= 1 && den.size() <= 1 && num.isNormal() && den.isNormal();
    }

    // Check whether both numerator and denominator are at most 2 limbs (and not 1 limb)
    inline bool isMediumRationalPair(const sangi::Rational&, const sangi::Rational&) {
        return false;  // 2-limb path disabled (regression investigation)
    }

    // Get an unsigned 64-bit value from an Int (0 if the value is 0)
    inline uint64_t toU64(const sangi::Int& x) {
        return x.isZero() ? 0 : x.word(0);
    }

    // =========================================================================
    // 128-bit unsigned arithmetic helpers (for the 2-limb fast path)
    // =========================================================================
    struct u128 {
        uint64_t lo, hi;
        bool isZero() const { return lo == 0 && hi == 0; }
    };

    inline u128 u128_from_int(const sangi::Int& x) {
        if (x.isZero()) return {0, 0};
        uint64_t lo = x.word(0);
        uint64_t hi = (x.size() >= 2) ? x.word(1) : 0;
        return {lo, hi};
    }

    inline u128 u128_mul(u128 a, u128 b) {
        // 128x128 → only the lower 128 bits of the lower 256 bits (enough for gcd)
        // However Rational multiplication needs 256 bits, so use a separate function
        uint64_t hi;
        uint64_t lo = _umul128(a.lo, b.lo, &hi);
        hi += a.lo * b.hi + a.hi * b.lo;
        return {lo, hi};
    }

    // 128x128 → 256-bit (4 limb)
    struct u256 {
        uint64_t w[4]; // w[0]=LSW
        bool isZero() const { return w[0]==0 && w[1]==0 && w[2]==0 && w[3]==0; }
        size_t wordCount() const {
            if (w[3]) return 4;
            if (w[2]) return 3;
            if (w[1]) return 2;
            if (w[0]) return 1;
            return 0;
        }
    };

    inline u256 u256_mul128(u128 a, u128 b) {
        // Full 256-bit product of (a.hi:a.lo) × (b.hi:b.lo)
        u256 r = {};
        uint64_t c;
        // a.lo * b.lo → r[1]:r[0]
        r.w[0] = _umul128(a.lo, b.lo, &r.w[1]);
        // a.lo * b.hi → accumulate into r[2]:r[1]
        uint64_t t_lo = _umul128(a.lo, b.hi, &c);
        uint64_t carry = 0;
        r.w[1] += t_lo; if (r.w[1] < t_lo) carry = 1;
        r.w[2] = c + carry;
        // a.hi * b.lo → accumulate into r[2]:r[1]
        t_lo = _umul128(a.hi, b.lo, &c);
        carry = 0;
        r.w[1] += t_lo; if (r.w[1] < t_lo) carry = 1;
        r.w[2] += c + carry; if (r.w[2] < c + carry) r.w[3] = 1;
        // a.hi * b.hi → accumulate into r[3]:r[2]
        t_lo = _umul128(a.hi, b.hi, &c);
        carry = 0;
        r.w[2] += t_lo; if (r.w[2] < t_lo) carry = 1;
        r.w[3] += c + carry;
        return r;
    }

    // 128-bit mod 64-bit → 64-bit remainder (safe version)
    inline uint64_t u128_mod64(u128 a, uint64_t d) {
        // (a.hi:a.lo) mod d
        // _udiv128 precondition: hi < d
        uint64_t r = a.hi % d;
        _udiv128(r, a.lo, d, &r);
        return r;
    }

    // 128-bit GCD (Euclidean algorithm)
    inline u128 u128_gcd(u128 a, u128 b) {
        // Fast path: when both fit in 1 limb
        if (a.hi == 0 && b.hi == 0) {
            return {std::gcd(a.lo, b.lo), 0};
        }
        while (!b.isZero()) {
            u128 r;
            if (a.hi == 0 && b.hi == 0) {
                r = {a.lo % b.lo, 0};
            } else if (b.hi == 0) {
                r = {u128_mod64(a, b.lo), 0};
            } else if (a.hi == 0) {
                // a < b, so a mod b = a → swap
                r = a;
            } else {
                // Both 128-bit: compute a mod b by repeated subtraction
                // (the GCD recursion depth is O(log(min(a,b))), so both being
                //  128-bit is rare, and within the first 1-2 iterations one side
                //  drops to 64-bit)
                r = a;
                while (r.hi > b.hi || (r.hi == b.hi && r.lo >= b.lo)) {
                    uint64_t borrow = (r.lo < b.lo) ? 1 : 0;
                    r.lo -= b.lo;
                    r.hi -= b.hi + borrow;
                }
            }
            a = b;
            b = r;
        }
        return a;
    }

    // 128-bit / 64-bit division → 128-bit quotient
    inline u128 u128_div64(u128 a, uint64_t d) {
        u128 q;
        q.hi = a.hi / d;
        uint64_t r_hi = a.hi % d;
        q.lo = _udiv128(r_hi, a.lo, d, &r_hi);
        return q;
    }

    // 256-bit / 128-bit → no estimate needed; for GCD where only the remainder is
    //   needed, handled separately
    // mod from 256-bit into u128: for reduction
    inline u128 u256_mod128(u256 a, u128 m) {
        // Take mod from the high words downward
        // a = (w3:w2:w1:w0), m = 128-bit
        // (w3:w2) mod m → r2, (r2:w1) mod m → r1, (r1:w0) mod m → result
        // Each step is 128-bit mod 128-bit
        u128 cur = {a.w[2], a.w[3]};
        // cur mod m
        u128 r = cur;
        if (r.hi > m.hi || (r.hi == m.hi && r.lo >= m.lo)) {
            r = u128_gcd(r, m); // ← this is gcd, not mod. Implement properly
            // Strictly mod is needed, but Rational reduction uses gcd(num, g) where g
            // is the original gcd, so g always fits in 64-bit. 128-bit mod 128-bit is
            // unnecessary.
        }
        // In practice, the gcd needed for Rational reduction is gcd(numerator, g) with
        // g ≤ 64-bit, so u256_mod128 is not used. Instead u256_mod64 is implemented.
        return {0, 0}; // placeholder
    }

    // 256-bit mod 64-bit → 64-bit remainder
    inline uint64_t u256_mod64(u256 a, uint64_t d) {
        // Compute (w3:w2:w1:w0) mod d from the high words downward
        uint64_t r = a.w[3] % d;
        _udiv128(r, a.w[2], d, &r);
        _udiv128(r, a.w[1], d, &r);
        _udiv128(r, a.w[0], d, &r);
        return r;
    }

    // 256-bit / 64-bit → 256-bit quotient
    inline u256 u256_div64(u256 a, uint64_t d) {
        u256 q = {};
        uint64_t r = 0;
        q.w[3] = _udiv128(r, a.w[3], d, &r);
        q.w[2] = _udiv128(r, a.w[2], d, &r);
        q.w[1] = _udiv128(r, a.w[1], d, &r);
        q.w[0] = _udiv128(r, a.w[0], d, &r);
        return q;
    }

    // 256-bit addition
    inline u256 u256_add(u256 a, u256 b) {
        u256 r = {};
        uint64_t c = 0;
        r.w[0] = a.w[0] + b.w[0]; c = (r.w[0] < a.w[0]) ? 1 : 0;
        r.w[1] = a.w[1] + b.w[1] + c; c = (r.w[1] < a.w[1] || (c && r.w[1] == a.w[1])) ? 1 : 0;
        r.w[2] = a.w[2] + b.w[2] + c; c = (r.w[2] < a.w[2] || (c && r.w[2] == a.w[2])) ? 1 : 0;
        r.w[3] = a.w[3] + b.w[3] + c;
        return r;
    }

    // 256-bit subtraction (assumes a >= b)
    inline u256 u256_sub(u256 a, u256 b) {
        u256 r = {};
        uint64_t borrow = 0;
        r.w[0] = a.w[0] - b.w[0]; borrow = (a.w[0] < b.w[0]) ? 1 : 0;
        r.w[1] = a.w[1] - b.w[1] - borrow; borrow = (a.w[1] < b.w[1] + borrow) ? 1 : 0;
        r.w[2] = a.w[2] - b.w[2] - borrow; borrow = (a.w[2] < b.w[2] + borrow) ? 1 : 0;
        r.w[3] = a.w[3] - b.w[3] - borrow;
        return r;
    }

    // 256-bit comparison (a >= b)
    inline bool u256_ge(u256 a, u256 b) {
        if (a.w[3] != b.w[3]) return a.w[3] > b.w[3];
        if (a.w[2] != b.w[2]) return a.w[2] > b.w[2];
        if (a.w[1] != b.w[1]) return a.w[1] > b.w[1];
        return a.w[0] >= b.w[0];
    }

    // Build a Rational from u256 (up to 4 limbs)
    // Uses the public constructor Rational(Int, Int, bool reduce)
    inline sangi::Rational makeRationalFromU256(u256 num, int num_sign, u256 den) {
        if (num.isZero()) return sangi::Rational();
        size_t nn = num.wordCount(), dn = den.wordCount();
        sangi::Int n_i = sangi::Int::fromRawWordsPreNormalized(
            std::span<const uint64_t>(num.w, nn), num_sign);
        sangi::Int d_i = sangi::Int::fromRawWordsPreNormalized(
            std::span<const uint64_t>(den.w, dn), 1);
        return sangi::Rational(std::move(n_i), std::move(d_i), false);
    }

    // Build a signed Rational (from 1-limb numerator and denominator)
    inline sangi::Rational makeSmallRational(uint64_t num, int num_sign, uint64_t den) {
        if (num == 0) return sangi::Rational();
        sangi::Int n(num);
        if (num_sign < 0) n = -n;
        return sangi::Rational(std::move(n), sangi::Int(den), false);
    }
} // anonymous namespace

namespace sangi {

    //==========================================================================
    // Constructors
    //==========================================================================

    Rational::Rational()
        : numerator_(Int(0))
        , denominator_(Int(1))
    {
    }

    Rational::Rational(int num, int den)
        : numerator_(Int(num))
        , denominator_(Int(den))
    {
        if (den == 0) {
            numerator_ = Int::NaN();
            denominator_ = Int(1);
            return;
        }
        reduce();
    }

    Rational::Rational(const Int& num)
        : numerator_(num)
        , denominator_(Int(1))
    {
    }

    Rational::Rational(const Int& num, const Int& den, bool doReduce)
        : numerator_(num)
        , denominator_(den)
    {
        if (denominator_.isZero()) {
            numerator_ = Int::NaN();
            denominator_ = Int(1);
            return;
        }
        if (doReduce) reduce();
    }

    Rational::Rational(Int&& num, Int&& den, bool doReduce)
        : numerator_(std::move(num))
        , denominator_(std::move(den))
    {
        if (denominator_.isZero()) {
            numerator_ = Int::NaN();
            denominator_ = Int(1);
            return;
        }
        if (doReduce) reduce();
    }

    Rational::Rational(const Rational& other)
        : numerator_(other.numerator_)
        , denominator_(other.denominator_)
    {
    }

    Rational::Rational(Rational&& other) noexcept
        : numerator_(std::move(other.numerator_))
        , denominator_(std::move(other.denominator_))
    {
    }

    // String constructor
    // Supported formats: "3/7", "42", "-1/3", "1.5" (decimal point only in base 10)
    Rational::Rational(std::string_view str, int base)
        : numerator_(Int(0))
        , denominator_(Int(1))
    {
        if (str.empty()) {
            numerator_ = Int::NaN();
            return;
        }

        // Look for "/" → fraction form
        auto slashPos = str.find('/');
        if (slashPos != std::string_view::npos) {
            auto numStr = str.substr(0, slashPos);
            auto denStr = str.substr(slashPos + 1);
            if (numStr.empty() || denStr.empty()) {
                numerator_ = Int::NaN();
                return;
            }
            numerator_ = Int(numStr, base);
            denominator_ = Int(denStr, base);
            if (denominator_.isZero()) {
                numerator_ = Int::NaN();
                denominator_ = Int(1);
                return;
            }
            reduce();
            return;
        }

        // Look for "." → decimal form (base 10 only)
        auto dotPos = str.find('.');
        if (dotPos != std::string_view::npos && base == 10) {
            auto intPart = str.substr(0, dotPos);
            auto decPart = str.substr(dotPos + 1);

            bool negative = false;
            if (!intPart.empty() && intPart[0] == '-') {
                negative = true;
                intPart = intPart.substr(1);
            }

            int n = static_cast<int>(decPart.length());
            Int intVal = intPart.empty() ? Int(0) : Int(intPart, 10);
            Int decVal = decPart.empty() ? Int(0) : Int(decPart, 10);
            Int den = sangi::pow(Int(10), static_cast<unsigned int>(n));

            numerator_ = intVal * den + decVal;
            denominator_ = den;
            if (negative) numerator_ = -numerator_;
            reduce();
            return;
        }

        // Integer form
        numerator_ = Int(str, base);
        denominator_ = Int(1);
    }

    Rational::Rational(double value)
        : numerator_(Int(0))
        , denominator_(Int(1))
    {
        *this = value;
    }

    Rational::Rational(float value)
        : numerator_(Int(0))
        , denominator_(Int(1))
    {
        *this = value;
    }

    // Arbitrary-precision float → rational (exact conversion)
    // Float = mantissa * 2^exponent (mantissa is positive, sign is held in is_negative_)
    Rational::Rational(const Float& value)
        : numerator_(Int(0))
        , denominator_(Int(1))
    {
        // Special states
        if (value.isNaN() || value.isInfinity()) {
            numerator_ = Int::NaN();
            return;
        }
        if (value.mantissa().isZero()) {
            return; // 0/1
        }

        Int m = value.mantissa(); // always positive
        int64_t e = value.exponent();

        if (e >= 0) {
            // mantissa * 2^e → integer
            numerator_ = m << static_cast<int>(e);
            denominator_ = Int(1);
        } else {
            // mantissa / 2^(-e) → rational
            numerator_ = m;
            denominator_ = Int(1) << static_cast<int>(-e);
            reduce();
        }

        if (value.isNegative()) {
            numerator_ = -numerator_;
        }
    }

    //==========================================================================
    // Assignment operators
    //==========================================================================

    Rational& Rational::operator=(const Rational& other) {
        numerator_ = other.numerator_;
        denominator_ = other.denominator_;
        return *this;
    }

    Rational& Rational::operator=(Rational&& other) noexcept {
        numerator_ = std::move(other.numerator_);
        denominator_ = std::move(other.denominator_);
        return *this;
    }

    Rational& Rational::operator=(int value) {
        numerator_ = Int(value);
        denominator_ = Int(1);
        return *this;
    }

    Rational& Rational::operator=(const Int& value) {
        numerator_ = value;
        denominator_ = Int(1);
        return *this;
    }

    // IEEE 754 double-precision float → rational
    // [Source] https://en.wikipedia.org/wiki/Double-precision_floating-point_format
    // [Improvement] Use std::bit_cast (C++23), fixed subnormal bug
    Rational& Rational::operator=(double value) {
        uint64_t bits = std::bit_cast<uint64_t>(value);

        constexpr uint64_t mantissaMask = (1ULL << 52) - 1;
        constexpr uint64_t exponentMask = (1ULL << 11) - 1;

        uint64_t mantissa = bits & mantissaMask;
        bits >>= 52;
        int exponent = static_cast<int>(bits & exponentMask);
        bits >>= 11;
        int sign = static_cast<int>(bits);

        // Infinity or NaN
        if (exponent == 2047) {
            numerator_ = Int::NaN();
            denominator_ = Int(1);
            return *this;
        }

        int e;
        if (exponent != 0) {
            // Normalized number: add the hidden bit
            mantissa |= 1ULL << 52;
            e = exponent - 1023 - 52;
        } else {
            // Subnormal number: effective exponent is 1 - 1023 (fixes sangi bug)
            e = 1 - 1023 - 52;  // = -1074
        }

        if (mantissa == 0) {
            numerator_ = Int(0);
            denominator_ = Int(1);
            return *this;
        }

        if (e > 0) {
            numerator_ = Int(mantissa) << e;
            denominator_ = Int(1);
        } else if (e < 0) {
            // The GCD of mantissa and 2^(-e) is a power of 2 only, so compute directly
            int tz = std::countr_zero(mantissa);
            int common = std::min(tz, -e);
            numerator_ = Int(mantissa >> common);
            denominator_ = Int(1) << (-e - common);
        } else {
            numerator_ = Int(mantissa);
            denominator_ = Int(1);
        }

        if (sign) {
            numerator_ = -numerator_;
        }

        return *this;
    }

    // IEEE 754 single-precision float → rational
    // [Source] https://en.wikipedia.org/wiki/Single-precision_floating-point_format
    Rational& Rational::operator=(float value) {
        uint32_t bits = std::bit_cast<uint32_t>(value);

        constexpr uint32_t mantissaMask = (1U << 23) - 1;
        constexpr uint32_t exponentMask = (1U << 8) - 1;

        uint32_t mantissa = bits & mantissaMask;
        bits >>= 23;
        int exponent = static_cast<int>(bits & exponentMask);
        bits >>= 8;
        int sign = static_cast<int>(bits);

        // Infinity or NaN
        if (exponent == 255) {
            numerator_ = Int::NaN();
            denominator_ = Int(1);
            return *this;
        }

        int e;
        if (exponent != 0) {
            mantissa |= 1U << 23;
            e = exponent - 127 - 23;
        } else {
            e = 1 - 127 - 23;  // = -149
        }

        if (mantissa == 0) {
            numerator_ = Int(0);
            denominator_ = Int(1);
            return *this;
        }

        if (e > 0) {
            numerator_ = Int(static_cast<uint64_t>(mantissa)) << e;
            denominator_ = Int(1);
        } else if (e < 0) {
            // The GCD of mantissa and 2^(-e) is a power of 2 only, so compute directly
            int tz = std::countr_zero(mantissa);
            int common = std::min(tz, -e);
            numerator_ = Int(static_cast<uint64_t>(mantissa >> common));
            denominator_ = Int(1) << (-e - common);
        } else {
            numerator_ = Int(static_cast<uint64_t>(mantissa));
            denominator_ = Int(1);
        }

        if (sign) {
            numerator_ = -numerator_;
        }

        return *this;
    }

    //==========================================================================
    // Unary operators
    //==========================================================================

    Rational Rational::operator-() const& {
        Rational result;
        result.numerator_ = -numerator_;
        result.denominator_ = denominator_;
        return result;
    }

    Rational Rational::operator-() && {
        numerator_ = -numerator_;
        return std::move(*this);
    }

    //==========================================================================
    // Reduction
    //==========================================================================

    bool Rational::reduce() {
        // Special-state check
        if (numerator_.isSpecialState() || denominator_.isSpecialState()) {
            return false;
        }

        // denominator==0 → NaN
        if (denominator_.isZero()) {
            numerator_ = Int::NaN();
            denominator_ = Int(1);
            return true;
        }

        // Adjust signs so the denominator is positive
        if (denominator_.isNegative()) {
            numerator_ = -numerator_;
            denominator_ = -denominator_;
        }

        // numerator==0 → normalize to 0/1
        if (numerator_.isZero()) {
            numerator_.setSign(0);  // normalize negative zero
            denominator_ = Int(1);
            return false;
        }

        // Reduce by GCD (gcd() handles absolute values internally, so abs() is unneeded)
        Int g = gcd(numerator_, denominator_);
        if (!g.isOne()) {
            // 1-limb fast path: divExactWord (Hensel) is ~2x vs divmod_1
            if (g.size() == 1) {
                uint64_t gw = g.word(0);
                IntOps::divExactWord(numerator_, gw);
                IntOps::divExactWord(denominator_, gw);
            } else {
                // multi-limb: mpn::divexact (Hensel lifting) — lighter than ordinary division
                IntOps::divExactInPlace(numerator_, g);
                IntOps::divExactInPlace(denominator_, g);
            }
            return true;
        }
        return false;
    }

    //==========================================================================
    // Conversion
    //==========================================================================

    double Rational::toDouble() const {
        return numerator_.toDouble() / denominator_.toDouble();
    }

    float Rational::toFloat() const {
        return static_cast<float>(toDouble());
    }

    Rational::operator int() const {
        return static_cast<int>(toDouble());
    }

    Rational::operator float() const {
        return toFloat();
    }

    Rational::operator double() const {
        return toDouble();
    }

    //==========================================================================
    // String conversion
    //==========================================================================

    std::string Rational::toString(int base) const {
        if (isNaN()) return "NaN";

        if (denominator_.isOne()) {
            return numerator_.toString(base);
        }
        return numerator_.toString(base) + "/" + denominator_.toString(base);
    }

    // Convert to a string containing a decimal point
    // [Ported from] sangi Rational::ToNumeric
    // [Fix] rounding carry bug (t[nn] → s[nn])
    std::string Rational::toDecimal(int digits) const {
        if (isNaN()) return "NaN";
        if (isZero()) return "0";
        if (digits <= 0) digits = 0;

        // Integer part
        Int num = abs(numerator_);
        Int q = num / denominator_;
        Int r = num % denominator_;

        std::string sign = numerator_.isNegative() ? "-" : "";
        std::string intPart = q.toString();

        if (r.isZero() || digits == 0) {
            return sign + intPart;
        }

        // Compute the fractional part one digit at a time
        std::string decPart;
        Int ten(10);
        for (int i = 0; i <= digits; i++) {  // compute one extra digit (for rounding)
            r = r * ten;
            Int digit = r / denominator_;
            r = r % denominator_;
            decPart += digit.toString();
        }

        // Rounding
        char roundDigit = decPart[digits];
        decPart.resize(digits);

        if (roundDigit >= '5') {
            bool carried = true;
            for (int i = digits - 1; i >= 0 && carried; i--) {
                if (decPart[i] == '9') {
                    decPart[i] = '0';
                } else {
                    decPart[i]++;
                    carried = false;
                }
            }
            if (carried) {
                // Carry into the integer part
                q += 1;
                intPart = q.toString();
            }
        }

        return sign + intPart + "." + decPart;
    }

    //==========================================================================
    // Stream output
    //==========================================================================

    std::ostream& operator<<(std::ostream& os, const Rational& r) {
        return os << r.toString();
    }

    std::istream& operator>>(std::istream& is, Rational& r) {
        std::string str;
        is >> str;
        if (!str.empty()) {
            r = Rational(str);
        }
        return is;
    }

    //==========================================================================
    // Comparison operators
    //==========================================================================

    std::partial_ordering operator<=>(const Rational& lhs, const Rational& rhs) {
        if (lhs.isNaN() || rhs.isNaN()) {
            return std::partial_ordering::unordered;
        }

        // Early decision by sign
        int lsign = lhs.numerator_.getSign();
        int rsign = rhs.numerator_.getSign();
        if (lsign != rsign) {
            if (lsign > rsign) return std::partial_ordering::greater;
            return std::partial_ordering::less;
        }
        // Both zero
        if (lsign == 0) return std::partial_ordering::equivalent;

        // Same sign: early decision by size (same technique as GMP mpq_cmp)
        // a/b vs c/d → a*d vs c*b
        // Approximate the size of the products by limb count
        size_t ad_size = lhs.numerator_.size() + rhs.denominator_.size();
        size_t cb_size = rhs.numerator_.size() + lhs.denominator_.size();
        if (ad_size != cb_size) {
            // Positive case: ad_size > cb_size → a*d > c*b → lhs > rhs
            // Negative case: reversed
            if (lsign > 0)
                return (ad_size > cb_size) ? std::partial_ordering::greater
                                           : std::partial_ordering::less;
            else
                return (ad_size > cb_size) ? std::partial_ordering::less
                                           : std::partial_ordering::greater;
        }

        // 1-limb fast path: avoid Int construction with a raw 128-bit comparison
        if (ad_size == 2) {
            uint64_t a = lhs.numerator_.word(0), d = rhs.denominator_.word(0);
            uint64_t c = rhs.numerator_.word(0), b = lhs.denominator_.word(0);
            // Compare a*d vs c*b in 128-bit
            uint64_t ad_hi, ad_lo, cb_hi, cb_lo;
            ad_lo = _umul128(a, d, &ad_hi);
            cb_lo = _umul128(c, b, &cb_hi);
            if (ad_hi != cb_hi) {
                bool gt = (ad_hi > cb_hi);
                if (lsign < 0) gt = !gt;
                return gt ? std::partial_ordering::greater : std::partial_ordering::less;
            }
            if (ad_lo == cb_lo) return std::partial_ordering::equivalent;
            bool gt = (ad_lo > cb_lo);
            if (lsign < 0) gt = !gt;
            return gt ? std::partial_ordering::greater : std::partial_ordering::less;
        }

        // 2-limb fast path: 128x128 → 256-bit comparison
        if (ad_size <= 4 && lhs.numerator_.size() <= 2 && lhs.denominator_.size() <= 2
            && rhs.numerator_.size() <= 2 && rhs.denominator_.size() <= 2) {
            u128 a128 = u128_from_int(lhs.numerator_);
            u128 d128 = u128_from_int(rhs.denominator_);
            u128 c128 = u128_from_int(rhs.numerator_);
            u128 b128 = u128_from_int(lhs.denominator_);
            u256 ad256 = u256_mul128(a128, d128);
            u256 cb256 = u256_mul128(c128, b128);
            if (u256_ge(ad256, cb256) && !u256_ge(cb256, ad256))  // ad > cb
                return (lsign > 0) ? std::partial_ordering::greater : std::partial_ordering::less;
            if (u256_ge(cb256, ad256) && !u256_ge(ad256, cb256))  // cb > ad
                return (lsign > 0) ? std::partial_ordering::less : std::partial_ordering::greater;
            return std::partial_ordering::equivalent;  // ad == cb
        }

        // ── P3-3: precise early decision by bitLength ──
        // Compare by bit length rather than limb count: decided if the difference is ≥ 2
        {
            int64_t ad_bits = static_cast<int64_t>(lhs.numerator_.bitLength())
                            + static_cast<int64_t>(rhs.denominator_.bitLength());
            int64_t cb_bits = static_cast<int64_t>(rhs.numerator_.bitLength())
                            + static_cast<int64_t>(lhs.denominator_.bitLength());
            int64_t bit_diff = ad_bits - cb_bits;
            if (bit_diff >= 2) {
                // a*d is definitely larger than c*b
                return (lsign > 0) ? std::partial_ordering::greater
                                   : std::partial_ordering::less;
            }
            if (bit_diff <= -2) {
                return (lsign > 0) ? std::partial_ordering::less
                                   : std::partial_ordering::greater;
            }
        }

        // ── P3-3: high-limb partial-product comparison ──
        // Without computing the full product, compare the products of the top 2 limbs for an early decision
        {
            const uint64_t* a_data = lhs.numerator_.data();
            size_t a_n = lhs.numerator_.size();
            const uint64_t* b_data = lhs.denominator_.data();
            size_t b_n = lhs.denominator_.size();
            const uint64_t* c_data = rhs.numerator_.data();
            size_t c_n = rhs.numerator_.size();
            const uint64_t* d_data = rhs.denominator_.data();
            size_t d_n = rhs.denominator_.size();

            // Compare by the product of the top 1 limb each (at most 128-bit × 2)
            uint64_t ad_top_hi, ad_top_lo, cb_top_hi, cb_top_lo;
            ad_top_lo = _umul128(a_data[a_n - 1], d_data[d_n - 1], &ad_top_hi);
            cb_top_lo = _umul128(c_data[c_n - 1], b_data[b_n - 1], &cb_top_hi);

            // Decided if the high products differ (ad_size == cb_size, so positions match)
            if (ad_top_hi != cb_top_hi) {
                bool gt = (ad_top_hi > cb_top_hi);
                if (lsign < 0) gt = !gt;
                return gt ? std::partial_ordering::greater
                          : std::partial_ordering::less;
            }
            // hi equal → compare lo (if the difference is ≥ 2, a low-word carry cannot flip it)
            if (ad_top_lo > cb_top_lo + 1) {
                return (lsign > 0) ? std::partial_ordering::greater
                                   : std::partial_ordering::less;
            }
            if (cb_top_lo > ad_top_lo + 1) {
                return (lsign > 0) ? std::partial_ordering::less
                                   : std::partial_ordering::greater;
            }
        }

        // Early decision by floating-point approximation
        // Compute a double approximation from the top few limbs and avoid the multiplication if far enough apart
        {
            double la = lhs.numerator_.toDouble();
            double lb = lhs.denominator_.toDouble();
            double ra = rhs.numerator_.toDouble();
            double rb = rhs.denominator_.toDouble();
            if (lb > 0.0 && rb > 0.0) {
                double lv = la / lb;
                double rv = ra / rb;
                // Decided if the relative difference is large enough (double precision: ~15 digits)
                double diff = lv - rv;
                double scale = std::max(std::abs(lv), std::abs(rv));
                if (scale > 0.0 && std::abs(diff) > scale * 1e-10) {
                    bool gt = (diff > 0.0);
                    return gt ? std::partial_ordering::greater
                              : std::partial_ordering::less;
                }
            }
        }

        // General case: compare exactly via cross multiplication (3-argument form reuses buffers)
        Int ad, cb;
        IntOps::mulUnchecked(lhs.numerator_, rhs.denominator_, ad);
        IntOps::mulUnchecked(lhs.denominator_, rhs.numerator_, cb);
        if (ad == cb) return std::partial_ordering::equivalent;
        return (ad > cb) ? std::partial_ordering::greater
                         : std::partial_ordering::less;
    }

    bool operator==(const Rational& lhs, const Rational& rhs) {
        if (lhs.isNaN() || rhs.isNaN()) return false;
        // Comparison of reduced fractions: check numerator and denominator each for equality
        return lhs.numerator_ == rhs.numerator_
            && lhs.denominator_ == rhs.denominator_;
    }

    // --- Rational vs Int ---
    std::partial_ordering operator<=>(const Rational& lhs, const Int& rhs) {
        if (lhs.isNaN() || rhs.isNaN()) {
            return std::partial_ordering::unordered;
        }
        // a/b vs n → a vs n*b
        // Early decision by sign
        int lsign = lhs.numerator_.getSign();
        int rsign = rhs.getSign();
        if (lsign != rsign) {
            if (lsign > rsign) return std::partial_ordering::greater;
            return std::partial_ordering::less;
        }
        if (lsign == 0) return std::partial_ordering::equivalent;

        // Early decision by size
        size_t a_size = lhs.numerator_.size();
        size_t nb_size = rhs.size() + lhs.denominator_.size();
        if (a_size != nb_size) {
            if (lsign > 0)
                return (a_size > nb_size) ? std::partial_ordering::greater
                                          : std::partial_ordering::less;
            else
                return (a_size > nb_size) ? std::partial_ordering::less
                                          : std::partial_ordering::greater;
        }

        // 1-limb fast path
        if (a_size == 1 && lhs.denominator_.size() == 1 && rhs.size() == 1) {
            uint64_t a = lhs.numerator_.word(0);
            uint64_t nb_hi, nb_lo;
            nb_lo = _umul128(rhs.word(0), lhs.denominator_.word(0), &nb_hi);
            if (nb_hi != 0) {
                // n*b > 64bit → a < n*b (a is 1 limb)
                return (lsign > 0) ? std::partial_ordering::less
                                   : std::partial_ordering::greater;
            }
            if (a == nb_lo) return std::partial_ordering::equivalent;
            bool gt = (a > nb_lo);
            if (lsign < 0) gt = !gt;
            return gt ? std::partial_ordering::greater : std::partial_ordering::less;
        }

        // General case
        Int nb = rhs * lhs.denominator_;
        if (lhs.numerator_ == nb) return std::partial_ordering::equivalent;
        return (lhs.numerator_ > nb) ? std::partial_ordering::greater
                                     : std::partial_ordering::less;
    }

    bool operator==(const Rational& lhs, const Int& rhs) {
        if (lhs.isNaN() || rhs.isNaN()) return false;
        // Reduced fraction: denominator is 1 and numerators are equal
        return lhs.denominator_.isOne() && lhs.numerator_ == rhs;
    }

    // --- Rational vs int ---
    std::partial_ordering operator<=>(const Rational& lhs, int rhs) {
        return lhs <=> Int(rhs);
    }

    bool operator==(const Rational& lhs, int rhs) {
        if (lhs.isNaN()) return false;
        return lhs.denominator_.isOne() && lhs.numerator_ == Int(rhs);
    }

    //==========================================================================
    // Arithmetic (Rational <-> Rational)
    // Addition: the Osawagawa method (reduces intermediate digit count with GCD)
    // [Source] http://www.oishi.info.waseda.ac.jp/~samukawa/LongintRational.pdf p47-48
    //==========================================================================

    Rational operator+(const Rational& lhs, const Rational& rhs) {
        // NaN propagation
        if (lhs.isNaN() || rhs.isNaN()) [[unlikely]] {
            Rational r; r.numerator_ = Int::NaN(); r.denominator_ = Int(1); return r;
        }
        // Double-word fast path: only when numerator and denominator are all 1 limb and
        // the intermediates (a*d_g, c*b_g) also fit in 1 limb
        if (isSmallRational(lhs.numerator_, lhs.denominator_) &&
            isSmallRational(rhs.numerator_, rhs.denominator_)) {
            uint64_t abs_a = toU64(lhs.numerator_), b = toU64(lhs.denominator_);
            uint64_t abs_c = toU64(rhs.numerator_), d = toU64(rhs.denominator_);
            if (abs_a == 0) return rhs;
            if (abs_c == 0) return lhs;
            int sign_a = lhs.numerator_.getSign();
            int sign_c = rhs.numerator_.getSign();

            uint64_t g = std::gcd(b, d);
            uint64_t d_g = d / g;  // d/g
            uint64_t b_g = b / g;  // b/g

            // Compute a * (d/g) and c * (b/g) in 128-bit
            uint64_t ad_hi, cb_hi;
            uint64_t ad_lo = _umul128(abs_a, d_g, &ad_hi);
            uint64_t cb_lo = _umul128(abs_c, b_g, &cb_hi);

            // Compute the numerator via 128-bit add/subtract
            uint64_t num_lo, num_hi;
            int num_sign;
            if (sign_a == sign_c) {
                // ad + cb (128-bit)
                num_lo = ad_lo + cb_lo;
                num_hi = ad_hi + cb_hi + (num_lo < ad_lo ? 1 : 0);
                num_sign = sign_a;
                // 129-bit overflow → general path
                // (with 1-limb inputs, num ≤ 2 * (2^64-1)^2 / g < 2^129, so practically never happens)
                // the top bit of num_hi may be set, but it is representable in a 2-limb Int
            } else {
                // |ad - cb| (128-bit)
                if (ad_hi > cb_hi || (ad_hi == cb_hi && ad_lo >= cb_lo)) {
                    uint64_t borrow = (ad_lo < cb_lo) ? 1 : 0;
                    num_lo = ad_lo - cb_lo;
                    num_hi = ad_hi - cb_hi - borrow;
                    num_sign = sign_a;
                } else {
                    uint64_t borrow = (cb_lo < ad_lo) ? 1 : 0;
                    num_lo = cb_lo - ad_lo;
                    num_hi = cb_hi - ad_hi - borrow;
                    num_sign = sign_c;
                }
            }
            if (num_lo == 0 && num_hi == 0) return Rational();

            // denominator = (b/g) * d (128-bit)
            uint64_t den_hi;
            uint64_t den_lo = _umul128(b_g, d, &den_hi);

            // Reduction: f = gcd(numerator, g)
            // g is 64-bit, so compute numerator mod g in 128-bit, then a 64-bit gcd
            if (g != 1) {
                uint64_t num_mod_g;
                if (num_hi == 0) {
                    num_mod_g = num_lo % g;
                } else {
                    // 128-bit mod: (num_hi:num_lo) mod g
                    // num_hi % g < g, so the _udiv128 precondition is satisfied
                    uint64_t hi_rem = num_hi % g;
                    _udiv128(hi_rem, num_lo, g, &num_mod_g);
                }
                uint64_t f = std::gcd(num_mod_g, g);
                if (f != 1) {
                    // num /= f, den = b_g * (d / f)
                    // 128-bit / 64-bit division
                    if (num_hi == 0) {
                        num_lo /= f;
                    } else {
                        // (num_hi:num_lo) / f
                        uint64_t q_hi = num_hi / f;
                        uint64_t r_hi = num_hi % f;
                        uint64_t rem;
                        // _udiv128(hi, lo, div, &rem) → returns the quotient (assumes hi < div)
                        num_lo = _udiv128(r_hi, num_lo, f, &rem);
                        num_hi = q_hi;
                    }
                    // den = b_g * (d / f) — d/f is 64-bit, b_g is 64-bit
                    den_lo = _umul128(b_g, d / f, &den_hi);
                }
            }

            // Build the result as a Rational (handles up to 2 limbs)
            if (num_hi == 0 && den_hi == 0) {
                return makeSmallRational(num_lo, num_sign, den_lo);
            }
            // 2 limb: build the Int directly
            uint64_t nw[2] = { num_lo, num_hi };
            uint64_t dw[2] = { den_lo, den_hi };
            size_t nn = (num_hi != 0) ? 2 : 1;
            size_t dn = (den_hi != 0) ? 2 : 1;
            Int num_i = Int::fromRawWordsPreNormalized(std::span<const uint64_t>(nw, nn), num_sign);
            Int den_i = Int::fromRawWordsPreNormalized(std::span<const uint64_t>(dw, dn), 1);
            Rational result;
            result.numerator_ = std::move(num_i);
            result.denominator_ = std::move(den_i);
            return result;
        }
        // Quad-word fast path: numerator and denominator at most 2 limbs (20-38 digits)
        // Compute GCD/multiplication directly with 128-bit arithmetic, avoiding Int construction overhead
        if (isMediumRationalPair(lhs, rhs)) {
            u128 abs_a = u128_from_int(lhs.numerator_);
            u128 b = u128_from_int(lhs.denominator_);
            u128 abs_c = u128_from_int(rhs.numerator_);
            u128 d = u128_from_int(rhs.denominator_);
            if (abs_a.isZero()) return rhs;
            if (abs_c.isZero()) return lhs;
            int sign_a = lhs.numerator_.getSign();
            int sign_c = rhs.numerator_.getSign();

            // g = gcd(b, d), 128-bit GCD
            u128 g = u128_gcd(b, d);
            // d_g = d / g, b_g = b / g (128-bit / 128-bit → simplified: g is usually small)
            u128 d_g, b_g;
            if (g.hi == 0 && g.lo == 1) {
                d_g = d; b_g = b;
            } else if (g.hi == 0) {
                d_g = u128_div64(d, g.lo);
                b_g = u128_div64(b, g.lo);
            } else {
                // when g is 128-bit, d/g and b/g are small → fall back to the general path
                goto general_add;
            }

            // Compute a * d_g and c * b_g in 256-bit
            u256 ad = u256_mul128(abs_a, d_g);
            u256 cb = u256_mul128(abs_c, b_g);

            // Compute the numerator via 256-bit add/subtract
            u256 num256;
            int num_sign;
            if (sign_a == sign_c) {
                num256 = u256_add(ad, cb);
                num_sign = sign_a;
            } else {
                if (u256_ge(ad, cb)) {
                    num256 = u256_sub(ad, cb);
                    num_sign = sign_a;
                } else {
                    num256 = u256_sub(cb, ad);
                    num_sign = sign_c;
                }
            }
            if (num256.isZero()) return Rational();

            // denominator = b_g * d (256-bit)
            u256 den256 = u256_mul128(b_g, d);

            // Reduction: f = gcd(num256 mod g, g) — g is 64-bit, so 256-bit mod 64-bit
            if (g.hi == 0 && g.lo != 1) {
                uint64_t num_mod_g = u256_mod64(num256, g.lo);
                uint64_t f = std::gcd(num_mod_g, g.lo);
                if (f != 1) {
                    num256 = u256_div64(num256, f);
                    den256 = u256_mul128(b_g, u128_div64(d, f));
                }
            }

            return makeRationalFromU256(num256, num_sign, den256);
        }
    general_add:
        if (lhs.denominator_ == rhs.denominator_) {
            // Same denominator: (a + c) / b
            Int num;
            IntOps::addUnchecked(lhs.numerator_, rhs.numerator_, num);
            return Rational(std::move(num), lhs.denominator_, true);
        }

        Int e = gcd(lhs.denominator_, rhs.denominator_);
        if (!e.isOne()) {
            // The Osawagawa method (3-argument form to reduce temporary Ints)
            Int d, t, t2;
            IntOps::divUnchecked(lhs.denominator_, e, d);
            IntOps::divUnchecked(rhs.denominator_, e, t);
            IntOps::mulUnchecked(lhs.numerator_, t, t);       // t = a * (d_rhs / e)
            IntOps::mulUnchecked(rhs.numerator_, d, t2);       // t2 = c * d
            IntOps::addUnchecked(t, t2, t);                    // t = a*(d_rhs/e) + c*d
            Int f = gcd(t, e);
            if (!f.isOne()) {
                Int num, den;
                IntOps::divUnchecked(t, f, num);
                IntOps::divUnchecked(rhs.denominator_, f, t2);
                IntOps::mulUnchecked(d, t2, den);
                return Rational(std::move(num), std::move(den), false);
            } else {
                Int den;
                IntOps::mulUnchecked(d, rhs.denominator_, den);
                return Rational(std::move(t), std::move(den), false);
            }
        } else {
            // gcd(b,d) == 1: the result is already reduced
            Int t1, t2, num, den;
            IntOps::mulUnchecked(lhs.numerator_, rhs.denominator_, t1);
            IntOps::mulUnchecked(lhs.denominator_, rhs.numerator_, t2);
            IntOps::addUnchecked(t1, t2, num);
            IntOps::mulUnchecked(lhs.denominator_, rhs.denominator_, den);
            return Rational(std::move(num), std::move(den), false);
        }
    }

    Rational operator-(const Rational& lhs, const Rational& rhs) {
        // NaN propagation
        if (lhs.isNaN() || rhs.isNaN()) [[unlikely]] {
            Rational r; r.numerator_ = Int::NaN(); r.denominator_ = Int(1); return r;
        }
        // Double-word fast path (same structure as addition, with rhs sign flipped)
        if (isSmallRational(lhs.numerator_, lhs.denominator_) &&
            isSmallRational(rhs.numerator_, rhs.denominator_)) {
            uint64_t abs_a = toU64(lhs.numerator_), b = toU64(lhs.denominator_);
            uint64_t abs_c = toU64(rhs.numerator_), d = toU64(rhs.denominator_);
            if (abs_c == 0) return lhs;
            if (abs_a == 0) return -rhs;
            int sign_a = lhs.numerator_.getSign();
            int sign_c = -(rhs.numerator_.getSign());  // sign flip

            uint64_t g = std::gcd(b, d);
            uint64_t d_g = d / g;
            uint64_t b_g = b / g;

            uint64_t ad_hi, cb_hi;
            uint64_t ad_lo = _umul128(abs_a, d_g, &ad_hi);
            uint64_t cb_lo = _umul128(abs_c, b_g, &cb_hi);

            // Compute the numerator via 128-bit add/subtract (same logic as operator+)
            uint64_t num_lo, num_hi;
            int num_sign;
            if (sign_a == sign_c) {
                num_lo = ad_lo + cb_lo;
                num_hi = ad_hi + cb_hi + (num_lo < ad_lo ? 1 : 0);
                num_sign = sign_a;
            } else {
                if (ad_hi > cb_hi || (ad_hi == cb_hi && ad_lo >= cb_lo)) {
                    uint64_t borrow = (ad_lo < cb_lo) ? 1 : 0;
                    num_lo = ad_lo - cb_lo;
                    num_hi = ad_hi - cb_hi - borrow;
                    num_sign = sign_a;
                } else {
                    uint64_t borrow = (cb_lo < ad_lo) ? 1 : 0;
                    num_lo = cb_lo - ad_lo;
                    num_hi = cb_hi - ad_hi - borrow;
                    num_sign = sign_c;
                }
            }
            if (num_lo == 0 && num_hi == 0) return Rational();

            uint64_t den_hi;
            uint64_t den_lo = _umul128(b_g, d, &den_hi);

            if (g != 1) {
                uint64_t num_mod_g;
                if (num_hi == 0) {
                    num_mod_g = num_lo % g;
                } else {
                    uint64_t hi_rem = num_hi % g;
                    _udiv128(hi_rem, num_lo, g, &num_mod_g);
                }
                uint64_t f = std::gcd(num_mod_g, g);
                if (f != 1) {
                    if (num_hi == 0) {
                        num_lo /= f;
                    } else {
                        uint64_t q_hi = num_hi / f;
                        uint64_t r_hi = num_hi % f;
                        uint64_t rem;
                        num_lo = _udiv128(r_hi, num_lo, f, &rem);
                        num_hi = q_hi;
                    }
                    den_lo = _umul128(b_g, d / f, &den_hi);
                }
            }

            if (num_hi == 0 && den_hi == 0) {
                return makeSmallRational(num_lo, num_sign, den_lo);
            }
            uint64_t nw[2] = { num_lo, num_hi };
            uint64_t dw[2] = { den_lo, den_hi };
            size_t nn = (num_hi != 0) ? 2 : 1;
            size_t dn = (den_hi != 0) ? 2 : 1;
            Int num_i = Int::fromRawWordsPreNormalized(std::span<const uint64_t>(nw, nn), num_sign);
            Int den_i = Int::fromRawWordsPreNormalized(std::span<const uint64_t>(dw, dn), 1);
            Rational result;
            result.numerator_ = std::move(num_i);
            result.denominator_ = std::move(den_i);
            return result;
        }
        // Quad-word fast path (same structure as operator+, with rhs sign flipped)
        if (isMediumRationalPair(lhs, rhs)) {
            u128 abs_a = u128_from_int(lhs.numerator_);
            u128 b = u128_from_int(lhs.denominator_);
            u128 abs_c = u128_from_int(rhs.numerator_);
            u128 d = u128_from_int(rhs.denominator_);
            if (abs_c.isZero()) return lhs;
            if (abs_a.isZero()) return -rhs;
            int sign_a = lhs.numerator_.getSign();
            int sign_c = -(rhs.numerator_.getSign());  // sign flip

            u128 g = u128_gcd(b, d);
            u128 d_g, b_g;
            if (g.hi == 0 && g.lo == 1) {
                d_g = d; b_g = b;
            } else if (g.hi == 0) {
                d_g = u128_div64(d, g.lo);
                b_g = u128_div64(b, g.lo);
            } else {
                goto general_sub;  // fall back to the general path
            }

            u256 ad = u256_mul128(abs_a, d_g);
            u256 cb = u256_mul128(abs_c, b_g);

            u256 num256;
            int num_sign;
            if (sign_a == sign_c) {
                num256 = u256_add(ad, cb);
                num_sign = sign_a;
            } else {
                if (u256_ge(ad, cb)) {
                    num256 = u256_sub(ad, cb);
                    num_sign = sign_a;
                } else {
                    num256 = u256_sub(cb, ad);
                    num_sign = sign_c;
                }
            }
            if (num256.isZero()) return Rational();

            u256 den256 = u256_mul128(b_g, d);

            if (g.hi == 0 && g.lo != 1) {
                uint64_t num_mod_g = u256_mod64(num256, g.lo);
                uint64_t f = std::gcd(num_mod_g, g.lo);
                if (f != 1) {
                    num256 = u256_div64(num256, f);
                    den256 = u256_mul128(b_g, u128_div64(d, f));
                }
            }

            return makeRationalFromU256(num256, num_sign, den256);
        }
    general_sub:
        if (lhs.denominator_ == rhs.denominator_) {
            return Rational(lhs.numerator_ - rhs.numerator_, lhs.denominator_, true);
        }

        Int e = gcd(lhs.denominator_, rhs.denominator_);
        if (!e.isOne()) {
            // The Osawagawa method (3-argument form to reduce temporary Ints) — same structure as operator+
            Int d, t, t2;
            IntOps::divUnchecked(lhs.denominator_, e, d);
            IntOps::divUnchecked(rhs.denominator_, e, t);
            IntOps::mulUnchecked(lhs.numerator_, t, t);       // t = a * (d_rhs / e)
            IntOps::mulUnchecked(rhs.numerator_, d, t2);       // t2 = c * d
            IntOps::subUnchecked(t, t2, t);                    // t = a*(d_rhs/e) - c*d
            Int f = gcd(t, e);
            if (!f.isOne()) {
                Int num, den;
                IntOps::divUnchecked(t, f, num);
                IntOps::divUnchecked(rhs.denominator_, f, t2);
                IntOps::mulUnchecked(d, t2, den);
                return Rational(std::move(num), std::move(den), false);
            } else {
                Int den;
                IntOps::mulUnchecked(d, rhs.denominator_, den);
                return Rational(std::move(t), std::move(den), false);
            }
        } else {
            // gcd(b,d) == 1: the result is already reduced
            Int t1, t2, num, den;
            IntOps::mulUnchecked(lhs.numerator_, rhs.denominator_, t1);
            IntOps::mulUnchecked(lhs.denominator_, rhs.numerator_, t2);
            IntOps::subUnchecked(t1, t2, num);
            IntOps::mulUnchecked(lhs.denominator_, rhs.denominator_, den);
            return Rational(std::move(num), std::move(den), false);
        }
    }

    Rational operator*(const Rational& lhs, const Rational& rhs) {
        // NaN propagation
        if (lhs.isNaN() || rhs.isNaN()) [[unlikely]] {
            Rational r; r.numerator_ = Int::NaN(); r.denominator_ = Int(1); return r;
        }
        // Double-word fast path: if numerator and denominator are all 1 limb, compute directly in 64/128-bit
        if (isSmallRational(lhs.numerator_, lhs.denominator_) &&
            isSmallRational(rhs.numerator_, rhs.denominator_)) {
            uint64_t a = toU64(lhs.numerator_), b = toU64(lhs.denominator_);
            uint64_t c = toU64(rhs.numerator_), d = toU64(rhs.denominator_);
            if (a == 0 || c == 0) return Rational();
            int sign = lhs.numerator_.getSign() * rhs.numerator_.getSign();
            // Pre-reduce via cross-GCD
            uint64_t g1 = std::gcd(a, d);
            uint64_t g2 = std::gcd(b, c);
            a /= g1; d /= g1;
            c /= g2; b /= g2;
            // 128-bit multiplication
            uint64_t num_hi, den_hi;
            uint64_t num_lo = _umul128(a, c, &num_hi);
            uint64_t den_lo = _umul128(b, d, &den_hi);
            if (num_hi == 0 && den_hi == 0) {
                return makeSmallRational(num_lo, sign, den_lo);
            }
            // When it fits in 2 limbs: build the Int directly
            uint64_t nw[2] = { num_lo, num_hi };
            uint64_t dw[2] = { den_lo, den_hi };
            std::span<const uint64_t> ns(nw, num_hi ? 2 : 1);
            std::span<const uint64_t> ds(dw, den_hi ? 2 : 1);
            Int num_i = Int::fromRawWordsPreNormalized(ns, sign);
            Int den_i = Int::fromRawWordsPreNormalized(ds, 1);
            Rational result;
            result.numerator_ = std::move(num_i);
            result.denominator_ = std::move(den_i);
            return result;
        }
        // Quad-word fast path: cross-GCD + multiplication for ≤ 2 limbs
        if (isMediumRationalPair(lhs, rhs)) {
            {
                u128 a = u128_from_int(lhs.numerator_);
                u128 b = u128_from_int(lhs.denominator_);
                u128 c = u128_from_int(rhs.numerator_);
                u128 d = u128_from_int(rhs.denominator_);
                if (a.isZero() || c.isZero()) return Rational();
                int sign = lhs.numerator_.getSign() * rhs.numerator_.getSign();

                // Cross-GCD: g1 = gcd(a,d), g2 = gcd(b,c)
                u128 g1 = u128_gcd(a, d);
                u128 g2 = u128_gcd(b, c);
                // Fast path only when g1, g2 fit in 64-bit
                if (g1.hi == 0 && g2.hi == 0) {
                    if (g1.lo != 1) { a = u128_div64(a, g1.lo); d = u128_div64(d, g1.lo); }
                    if (g2.lo != 1) { b = u128_div64(b, g2.lo); c = u128_div64(c, g2.lo); }
                    // result: num = a*c, den = b*d (at most 256-bit)
                    u256 num256 = u256_mul128(a, c);
                    u256 den256 = u256_mul128(b, d);
                    return makeRationalFromU256(num256, sign, den256);
                }
            }
        }

        // GCD optimization: cross-reduce, then multiply
        // gcd() handles absolute values internally, so abs() is unneeded
        Int g1 = gcd(lhs.numerator_, rhs.denominator_);
        Int g2 = gcd(lhs.denominator_, rhs.numerator_);

        Int num, den;
        if (g1.isOne() && g2.isOne()) {
            IntOps::mulUnchecked(lhs.numerator_, rhs.numerator_, num);
            IntOps::mulUnchecked(lhs.denominator_, rhs.denominator_, den);
        } else if (g1.isOne()) {
            Int t; IntOps::divUnchecked(rhs.numerator_, g2, t);
            IntOps::mulUnchecked(lhs.numerator_, t, num);
            IntOps::divUnchecked(lhs.denominator_, g2, t);
            IntOps::mulUnchecked(t, rhs.denominator_, den);
        } else if (g2.isOne()) {
            Int t; IntOps::divUnchecked(lhs.numerator_, g1, t);
            IntOps::mulUnchecked(t, rhs.numerator_, num);
            IntOps::divUnchecked(rhs.denominator_, g1, t);
            IntOps::mulUnchecked(lhs.denominator_, t, den);
        } else {
            Int t1, t2;
            IntOps::divUnchecked(lhs.numerator_, g1, t1);
            IntOps::divUnchecked(rhs.numerator_, g2, t2);
            IntOps::mulUnchecked(t1, t2, num);
            IntOps::divUnchecked(lhs.denominator_, g2, t1);
            IntOps::divUnchecked(rhs.denominator_, g1, t2);
            IntOps::mulUnchecked(t1, t2, den);
        }

        // Adjust the sign (when the denominator is negative)
        if (den.isNegative()) {
            num = -num;
            den = -den;
        }
        Rational result;
        result.numerator_ = std::move(num);
        result.denominator_ = std::move(den);
        return result;
    }

    Rational operator/(const Rational& lhs, const Rational& rhs) {
        // NaN propagation
        if (lhs.isNaN() || rhs.isNaN()) [[unlikely]] {
            Rational r; r.numerator_ = Int::NaN(); r.denominator_ = Int(1); return r;
        }
        // a/b ÷ c/d = a/b × d/c = (a*d)/(b*c)
        if (rhs.numerator_.isZero()) {
            Rational result;
            result.numerator_ = Int::NaN();
            result.denominator_ = Int(1);
            return result;
        }

        // Double-word fast path
        if (isSmallRational(lhs.numerator_, lhs.denominator_) &&
            isSmallRational(rhs.numerator_, rhs.denominator_)) {
            uint64_t a = toU64(lhs.numerator_), b = toU64(lhs.denominator_);
            uint64_t c = toU64(rhs.numerator_), d = toU64(rhs.denominator_);
            if (a == 0) return Rational();
            // a/b ÷ c/d = (a*d)/(b*c) — Cross-GCD: gcd(a,c) and gcd(b,d)
            int sign = lhs.numerator_.getSign() * rhs.numerator_.getSign();
            uint64_t g1 = std::gcd(a, c);
            uint64_t g2 = std::gcd(b, d);
            a /= g1; c /= g1;
            d /= g2; b /= g2;
            uint64_t num_hi, den_hi;
            uint64_t num_lo = _umul128(a, d, &num_hi);
            uint64_t den_lo = _umul128(b, c, &den_hi);
            if (num_hi == 0 && den_hi == 0) {
                return makeSmallRational(num_lo, sign, den_lo);
            }
            uint64_t nw[2] = { num_lo, num_hi };
            uint64_t dw[2] = { den_lo, den_hi };
            std::span<const uint64_t> ns(nw, num_hi ? 2 : 1);
            std::span<const uint64_t> ds(dw, den_hi ? 2 : 1);
            Int num_i = Int::fromRawWordsPreNormalized(ns, sign);
            Int den_i = Int::fromRawWordsPreNormalized(ds, 1);
            Rational result;
            result.numerator_ = std::move(num_i);
            result.denominator_ = std::move(den_i);
            return result;
        }
        // Quad-word fast path: ≤ 2 limbs
        if (isMediumRationalPair(lhs, rhs)) {
            {
                u128 a = u128_from_int(lhs.numerator_);
                u128 b = u128_from_int(lhs.denominator_);
                u128 c = u128_from_int(rhs.numerator_);
                u128 d = u128_from_int(rhs.denominator_);
                if (a.isZero()) return Rational();
                int sign = lhs.numerator_.getSign() * rhs.numerator_.getSign();

                // Cross-GCD: g1 = gcd(a,c), g2 = gcd(b,d)
                u128 g1 = u128_gcd(a, c);
                u128 g2 = u128_gcd(b, d);
                if (g1.hi == 0 && g2.hi == 0) {
                    if (g1.lo != 1) { a = u128_div64(a, g1.lo); c = u128_div64(c, g1.lo); }
                    if (g2.lo != 1) { b = u128_div64(b, g2.lo); d = u128_div64(d, g2.lo); }
                    // num = a*d, den = b*c
                    u256 num256 = u256_mul128(a, d);
                    u256 den256 = u256_mul128(b, c);
                    return makeRationalFromU256(num256, sign, den256);
                }
            }
        }

        // Multiply by the reciprocal (cross-reduction optimization)
        // gcd() handles absolute values internally, so abs() is unneeded
        Int g1 = gcd(lhs.numerator_, rhs.numerator_);
        Int g2 = gcd(lhs.denominator_, rhs.denominator_);

        Int num, den;
        if (g1.isOne() && g2.isOne()) {
            IntOps::mulUnchecked(lhs.numerator_, rhs.denominator_, num);
            IntOps::mulUnchecked(lhs.denominator_, rhs.numerator_, den);
        } else {
            Int a_div_g1, d_div_g2, b_div_g2, c_div_g1;
            IntOps::divUnchecked(lhs.numerator_, g1, a_div_g1);
            IntOps::divUnchecked(rhs.denominator_, g2, d_div_g2);
            IntOps::mulUnchecked(a_div_g1, d_div_g2, num);
            IntOps::divUnchecked(lhs.denominator_, g2, b_div_g2);
            IntOps::divUnchecked(rhs.numerator_, g1, c_div_g1);
            IntOps::mulUnchecked(b_div_g2, c_div_g1, den);
        }

        if (den.isNegative()) {
            num.negate();
            den.negate();
        }
        Rational result;
        result.numerator_ = std::move(num);
        result.denominator_ = std::move(den);
        return result;
    }

    Rational operator%(const Rational& lhs, const Rational& rhs) {
        // a % b = a - floor(a/b) * b
        Int m = floor(lhs / rhs);
        return lhs - Rational(m) * rhs;
    }

    //==========================================================================
    // Arithmetic (Rational <-> Int)
    //==========================================================================

    Rational operator+(const Rational& lhs, const Int& rhs) {
        // a/b + n = (a + n*b) / b  (no reduction needed: gcd(a,b)==1 and gcd(n*b,b)==b)
        Int t, num;
        IntOps::mulUnchecked(rhs, lhs.denominator_, t);
        IntOps::addUnchecked(lhs.numerator_, t, num);
        return Rational(std::move(num), lhs.denominator_, false);
    }

    Rational operator+(const Int& lhs, const Rational& rhs) {
        return rhs + lhs;
    }

    Rational operator-(const Rational& lhs, const Int& rhs) {
        // a/b - n = (a - n*b) / b
        Int t, num;
        IntOps::mulUnchecked(rhs, lhs.denominator_, t);
        IntOps::subUnchecked(lhs.numerator_, t, num);
        return Rational(std::move(num), lhs.denominator_, false);
    }

    Rational operator-(const Int& lhs, const Rational& rhs) {
        // n - a/b = (n*b - a) / b
        Int t, num;
        IntOps::mulUnchecked(lhs, rhs.denominator_, t);
        IntOps::subUnchecked(t, rhs.numerator_, num);
        return Rational(std::move(num), rhs.denominator_, false);
    }

    Rational operator*(const Rational& lhs, const Int& rhs) {
        // Power of 2 → left-shift the numerator or right-shift the denominator (avoid GCD computation)
        if (!rhs.isZero() && !rhs.isNegative() && rhs.bitLength() <= 64) {
            uint64_t w = rhs.data()[0];
            if ((w & (w - 1)) == 0) {
                unsigned long shift = 0;
                _BitScanForward64(&shift, w);
                // Generalize the same logic as doubled(): if the denominator's low bits are 0, right-shift
                size_t den_tz = lhs.denominator_.countTrailingZeros();
                if (den_tz >= shift) {
                    return Rational(lhs.numerator_, lhs.denominator_ >> shift, false);
                } else {
                    // Right-shift the denominator by den_tz, left-shift the remainder into the numerator
                    int num_shift = static_cast<int>(shift - den_tz);
                    Int new_num = lhs.numerator_ << num_shift;
                    Int new_den = (den_tz > 0) ? (lhs.denominator_ >> static_cast<int>(den_tz)) : lhs.denominator_;
                    return Rational(std::move(new_num), std::move(new_den), false);
                }
            }
        }
        // (a/b) * n: pre-reduce via gcd(n, b)
        Int g = gcd(rhs, lhs.denominator_);
        if (g.isOne()) {
            Int num;
            IntOps::mulUnchecked(lhs.numerator_, rhs, num);
            return Rational(std::move(num), lhs.denominator_, false);
        } else {
            Int n_div_g, num, den;
            IntOps::divUnchecked(rhs, g, n_div_g);
            IntOps::mulUnchecked(lhs.numerator_, n_div_g, num);
            IntOps::divUnchecked(lhs.denominator_, g, den);
            return Rational(std::move(num), std::move(den), false);
        }
    }

    Rational operator*(const Int& lhs, const Rational& rhs) {
        return rhs * lhs;
    }

    Rational operator/(const Rational& lhs, const Int& rhs) {
        if (rhs.isZero()) {
            Rational result;
            result.numerator_ = Int::NaN();
            result.denominator_ = Int(1);
            return result;
        }
        // Power of 2 → right-shift the numerator or left-shift the denominator (avoid GCD computation)
        if (!rhs.isNegative() && rhs.bitLength() <= 64) {
            uint64_t w = rhs.data()[0];
            if ((w & (w - 1)) == 0) {
                unsigned long shift = 0;
                _BitScanForward64(&shift, w);
                // Generalize the same logic as halved()
                size_t num_tz = lhs.numerator_.countTrailingZeros();
                if (num_tz >= shift) {
                    return Rational(lhs.numerator_ >> shift, lhs.denominator_, false);
                } else {
                    int den_shift = static_cast<int>(shift - num_tz);
                    Int new_num = (num_tz > 0) ? (lhs.numerator_ >> static_cast<int>(num_tz)) : lhs.numerator_;
                    Int new_den = lhs.denominator_ << den_shift;
                    return Rational(std::move(new_num), std::move(new_den), false);
                }
            }
        }
        // (a/b) / n: pre-reduce via gcd(a, n)
        Int g = gcd(lhs.numerator_, rhs);
        Rational result;
        if (g.isOne()) {
            result.numerator_ = lhs.numerator_;
            IntOps::mulUnchecked(lhs.denominator_, rhs, result.denominator_);
        } else {
            IntOps::divUnchecked(lhs.numerator_, g, result.numerator_);
            Int n_div_g;
            IntOps::divUnchecked(rhs, g, n_div_g);
            IntOps::mulUnchecked(lhs.denominator_, n_div_g, result.denominator_);
        }
        if (result.denominator_.isNegative()) {
            result.numerator_.negate();
            result.denominator_.negate();
        }
        return result;
    }

    Rational operator/(const Int& lhs, const Rational& rhs) {
        if (rhs.numerator_.isZero()) {
            Rational result;
            result.numerator_ = Int::NaN();
            result.denominator_ = Int(1);
            return result;
        }
        // n / (a/b) = n*b / a: pre-reduce via gcd(n, a)
        Int g = gcd(lhs, rhs.numerator_);
        Rational result;
        if (g.isOne()) {
            IntOps::mulUnchecked(lhs, rhs.denominator_, result.numerator_);
            result.denominator_ = rhs.numerator_;
        } else {
            Int n_div_g;
            IntOps::divUnchecked(lhs, g, n_div_g);
            IntOps::mulUnchecked(n_div_g, rhs.denominator_, result.numerator_);
            IntOps::divUnchecked(rhs.numerator_, g, result.denominator_);
        }
        if (result.denominator_.isNegative()) {
            result.numerator_.negate();
            result.denominator_.negate();
        }
        return result;
    }

    //==========================================================================
    // Arithmetic (Rational <-> int)
    //==========================================================================

    Rational operator+(const Rational& lhs, int rhs) {
        return lhs + Int(rhs);
    }

    Rational operator+(int lhs, const Rational& rhs) {
        return Int(lhs) + rhs;
    }

    Rational operator-(const Rational& lhs, int rhs) {
        return lhs - Int(rhs);
    }

    Rational operator-(int lhs, const Rational& rhs) {
        return Int(lhs) - rhs;
    }

    Rational operator*(const Rational& lhs, int rhs) {
        return lhs * Int(rhs);
    }

    Rational operator*(int lhs, const Rational& rhs) {
        return Int(lhs) * rhs;
    }

    Rational operator/(const Rational& lhs, int rhs) {
        return lhs / Int(rhs);
    }

    Rational operator/(int lhs, const Rational& rhs) {
        return Int(lhs) / rhs;
    }

    //==========================================================================
    // Compound assignment
    //==========================================================================

    Rational& Rational::operator+=(const Rational& rhs) {
        *this = *this + rhs;
        return *this;
    }

    Rational& Rational::operator-=(const Rational& rhs) {
        *this = *this - rhs;
        return *this;
    }

    Rational& Rational::operator*=(const Rational& rhs) {
        *this = *this * rhs;
        return *this;
    }

    Rational& Rational::operator/=(const Rational& rhs) {
        *this = *this / rhs;
        return *this;
    }

    Rational& Rational::operator+=(const Int& rhs) {
        *this = *this + rhs;
        return *this;
    }

    Rational& Rational::operator-=(const Int& rhs) {
        *this = *this - rhs;
        return *this;
    }

    Rational& Rational::operator*=(const Int& rhs) {
        // (a/b) * n: pre-reduce via gcd(n, b)
        Int g = gcd(rhs, denominator_);
        if (g.isOne()) {
            Int tmp;
            IntOps::mulUnchecked(numerator_, rhs, tmp);
            numerator_ = std::move(tmp);
        } else {
            Int n_div_g;
            IntOps::divUnchecked(rhs, g, n_div_g);
            Int tmp;
            IntOps::mulUnchecked(numerator_, n_div_g, tmp);
            numerator_ = std::move(tmp);
            IntOps::divUnchecked(denominator_, g, denominator_);
        }
        return *this;
    }

    Rational& Rational::operator/=(const Int& rhs) {
        if (rhs.isZero()) {
            numerator_ = Int::NaN();
            denominator_ = Int(1);
            return *this;
        }
        // (a/b) / n: pre-reduce via gcd(a, n)
        Int g = gcd(numerator_, rhs);
        if (g.isOne()) {
            Int tmp;
            IntOps::mulUnchecked(denominator_, rhs, tmp);
            denominator_ = std::move(tmp);
        } else {
            IntOps::divUnchecked(numerator_, g, numerator_);
            Int n_div_g, tmp;
            IntOps::divUnchecked(rhs, g, n_div_g);
            IntOps::mulUnchecked(denominator_, n_div_g, tmp);
            denominator_ = std::move(tmp);
        }
        if (denominator_.isNegative()) {
            numerator_.negate();
            denominator_.negate();
        }
        return *this;
    }

    Rational& Rational::operator+=(int rhs) {
        return *this += Int(rhs);
    }

    Rational& Rational::operator-=(int rhs) {
        return *this -= Int(rhs);
    }

    Rational& Rational::operator*=(int rhs) {
        return *this *= Int(rhs);
    }

    Rational& Rational::operator/=(int rhs) {
        return *this /= Int(rhs);
    }

    //==========================================================================
    // Increment / decrement
    //==========================================================================

    Rational& Rational::operator++() {
        // a/b + 1 = (a+b)/b (no reduction needed: if gcd(a,b)==1 then gcd(a+b,b)==1)
        if (numerator_.isZero()) {
            numerator_ = Int(1);
            denominator_ = Int(1);
        } else {
            numerator_ += denominator_;
        }
        return *this;
    }

    Rational Rational::operator++(int) {
        Rational old = *this;
        ++(*this);
        return old;
    }

    Rational& Rational::operator--() {
        if (numerator_.isZero()) {
            numerator_ = Int(-1);
            denominator_ = Int(1);
        } else {
            numerator_ -= denominator_;
        }
        return *this;
    }

    Rational Rational::operator--(int) {
        Rational old = *this;
        --(*this);
        return old;
    }

    //==========================================================================
    // Exponentiation
    //==========================================================================

    Rational Rational::pow(int exponent) const {
        Rational result;
        if (exponent > 0) {
            result.numerator_ = sangi::pow(numerator_, static_cast<unsigned int>(exponent));
            result.denominator_ = sangi::pow(denominator_, static_cast<unsigned int>(exponent));
        } else if (exponent < 0) {
            // A negative power of 0 is undefined → NaN
            if (isZero()) {
                result.numerator_ = Int::NaN();
                result.denominator_ = Int(1);
                return result;
            }
            unsigned int absExp = static_cast<unsigned int>(-exponent);
            result.numerator_ = sangi::pow(denominator_, absExp);
            result.denominator_ = sangi::pow(numerator_, absExp);
            // Adjust the sign if the denominator became negative
            if (result.denominator_.isNegative()) {
                result.numerator_ = -result.numerator_;
                result.denominator_ = -result.denominator_;
            }
        } else {
            // x^0 = 1
            result.numerator_ = Int(1);
            result.denominator_ = Int(1);
        }
        return result;
    }

    //==========================================================================
    // Shift (multiply/divide by powers of 2)
    //==========================================================================

    Rational Rational::operator<<(int n) const {
        if (n > 0) {
            return Rational(numerator_ << n, denominator_, true);
        } else if (n < 0) {
            return Rational(numerator_, denominator_ << (-n), true);
        }
        return *this;
    }

    Rational Rational::operator>>(int n) const {
        if (n > 0) {
            return Rational(numerator_, denominator_ << n, true);
        } else if (n < 0) {
            return Rational(numerator_ << (-n), denominator_, true);
        }
        return *this;
    }

    //==========================================================================
    // Phase 3: free functions
    //==========================================================================

    Rational abs(const Rational& x) {
        if (x.numerator().isNegative()) {
            return -x;
        }
        return x;
    }

    Int floor(const Rational& x) {
        Int q = x.numerator() / x.denominator();
        if (x.numerator().isNegative()) {
            Int r = x.numerator() % x.denominator();
            if (!r.isZero()) {
                q -= 1;
            }
        }
        return q;
    }

    Int ceil(const Rational& x) {
        Int q = x.numerator() / x.denominator();
        if (x.numerator().isPositive()) {
            Int r = x.numerator() % x.denominator();
            if (!r.isZero()) {
                q += 1;
            }
        }
        return q;
    }

    Rational square(const Rational& x) {
        return x * x;
    }

    Int height(const Rational& x) {
        Int n = abs(x.numerator());
        Int d = abs(x.denominator());
        return (n >= d) ? n : d;
    }

    Rational conj(const Rational& x) {
        return x;
    }

    // GCD of rationals (Wolfram definition)
    // gcd(r1, r2, ...) = the largest rational r such that every ri/r is an integer
    Rational gcd(const Rational& x, const Rational& y) {
        Int d = x.denominator() * y.denominator();
        Int n1 = x.numerator() * y.denominator();
        Int n2 = y.numerator() * x.denominator();
        return Rational(gcd(n1, n2), d);
    }

    //==========================================================================
    // swap
    //==========================================================================

    void swap(Rational& a, Rational& b) noexcept {
        using std::swap;
        swap(a.numerator_, b.numerator_);
        swap(a.denominator_, b.denominator_);
    }

    //==========================================================================
    // Utility functions: inv, trunc, sgn, frac, mediant
    //==========================================================================

    Rational inv(const Rational& x) {
        if (x.isNaN()) return x;
        // Inside the Rational(den, num) constructor:
        //   num==0 → NaN (reciprocal of zero)
        //   sign adjustment + reduction are handled by reduce()
        return Rational(x.denominator(), x.numerator());
    }

    Int trunc(const Rational& x) {
        // Truncation toward zero = integer division (C++ division is toward zero)
        return x.numerator() / x.denominator();
    }

    Int round(const Rational& x) {
        // Round to nearest integer (half-away-from-zero: 0.5 → 1, -0.5 → -1)
        // Take floor of |x| + 1/2 and apply the original sign
        Rational half(1, 2);
        if (x.isNegative()) {
            return -floor(-x + half);
        }
        return floor(x + half);
    }

    int sgn(const Rational& x) {
        if (x.isPositive()) return 1;
        if (x.isNegative()) return -1;
        return 0;
    }

    Rational frac(const Rational& x) {
        // Fractional part: x - floor(x) (always 0 <= frac(x) < 1)
        return x - Rational(floor(x));
    }

    Rational mediant(const Rational& a, const Rational& b) {
        // Mediant: (a_num + b_num) / (a_den + b_den)
        return Rational(a.numerator() + b.numerator(),
                        a.denominator() + b.denominator());
    }

    //==========================================================================
    // Convenience functions
    //==========================================================================

    Rational Rational::doubled() const {
        // x*2: if the denominator is even, halve it; otherwise double the numerator
        if (!denominator_.getBit(0)) {
            // even (least significant bit == 0)
            return Rational(numerator_, denominator_ >> 1, false);
        } else {
            return Rational(numerator_ << 1, denominator_, false);
        }
    }

    Rational Rational::halved() const {
        // x/2: if the numerator is even, halve it; otherwise double the denominator
        // getBit(0) does not depend on the sign, so it can be used directly on numerator_
        if (!numerator_.getBit(0)) {
            // even (least significant bit == 0)
            return Rational(numerator_ >> 1, denominator_, false);
        } else {
            return Rational(numerator_, denominator_ << 1, false);
        }
    }

    //==========================================================================
    // Rational → Float conversion (exact conversion)
    //==========================================================================

    Float Rational::toMpFloat(int precision) const {
        if (isNaN()) return Float::nan();

        // precision == 0 uses the default precision
        if (precision <= 0) {
            precision = Float::defaultPrecision();
        }

        if (isZero()) return Float::zero(precision);

        // Shift the numerator sufficiently before dividing to get a precision-bit mantissa
        Int num = abs(numerator_);
        Int den = denominator_;

        int numBits = static_cast<int>(num.bitLength());
        int denBits = static_cast<int>(den.bitLength());

        // Reserve bits for precision + guard bits
        int shift = precision + 10 - (numBits - denBits);
        if (shift < 0) shift = 0;

        Int shifted = num << shift;
        Int quotient = shifted / den;

        // Float(mantissa, exponent, is_negative)
        Float result(quotient, -static_cast<int64_t>(shift), isNegative());
        result.setPrecision(precision);
        return result;
    }

    //==========================================================================
    // Continued fractions
    //==========================================================================

    std::vector<Int> Rational::toContinuedFraction() const {
        std::vector<Int> result;
        if (isNaN()) return result;

        Int num = numerator_;   // may be negative
        Int den = denominator_; // always positive

        while (!den.isZero()) {
            // Floor division: q = floor(num / den)
            Int q = num / den;
            Int r = num - q * den;  // C++ truncation remainder

            // Correct toward floor: if r < 0 then q--, r += den
            if (r.isNegative()) {
                q -= 1;
                r = r + den;
            }

            result.push_back(q);
            num = den;
            den = r;
        }

        return result;
    }

    Rational Rational::fromContinuedFraction(const std::vector<Int>& cf) {
        if (cf.empty()) return Rational(0);

        // Compute from the back: [a0; a1, a2, ..., an]
        // = a0 + 1/(a1 + 1/(a2 + ... + 1/an))
        // Efficient method: compute the convergents by forward iteration
        // h_{-1} = 1, h_0 = a_0
        // k_{-1} = 0, k_0 = 1
        // h_i = a_i * h_{i-1} + h_{i-2}
        // k_i = a_i * k_{i-1} + k_{i-2}
        Int h_prev2(1), h_prev1(cf[0]);
        Int k_prev2(0), k_prev1(1);

        for (size_t i = 1; i < cf.size(); ++i) {
            Int h = cf[i] * h_prev1 + h_prev2;
            Int k = cf[i] * k_prev1 + k_prev2;
            h_prev2 = h_prev1;
            h_prev1 = h;
            k_prev2 = k_prev1;
            k_prev1 = k;
        }

        return Rational(h_prev1, k_prev1);
    }

    //==========================================================================
    // FLINT-23: Rational Reconstruction
    //==========================================================================
    // Reconstruct p/q from a mod m via the extended Euclidean algorithm
    // |p| ≤ N, 0 < q ≤ D, p/q ≡ a (mod m)

    Rational reconstructRational(const Int& a, const Int& m,
                                 const Int& N, const Int& D)
    {
        // Extended-GCD-style approach:
        // r_0 = m, r_1 = a mod m
        // Run the Euclidean algorithm until |r_i| ≤ N
        // At that point (t_i is q, r_i is p)
        Int r0 = m, r1 = a % m;
        if (r1.isNegative()) r1 += m;
        Int t0(0), t1(1);

        while (r1 > N) {
            Int q = r0 / r1;
            Int r2 = r0 - q * r1;
            Int t2 = t0 - q * t1;
            r0 = r1; r1 = r2;
            t0 = t1; t1 = t2;
        }

        // p = r1, q = t1
        Int p = r1;
        Int q = t1;
        if (q.isNegative()) { p = -p; q = -q; }

        // Check q ≤ D
        if (q > D) return Rational();  // reconstruction failed

        return Rational(p, q);
    }

    //==========================================================================
    // FLINT-24: Harmonic Numbers
    //==========================================================================
    // H_n = 1 + 1/2 + ... + 1/n
    // Divide and conquer: H(lo, hi) recurses to log(n) depth, reducing the number of multiplications

    static Rational harmonicRange(unsigned int lo, unsigned int hi) {
        if (lo > hi) return Rational(0);
        if (lo == hi) return Rational(1, (int)lo);
        unsigned int mid = lo + (hi - lo) / 2;
        return harmonicRange(lo, mid) + harmonicRange(mid + 1, hi);
    }

    Rational harmonicNumber(unsigned int n) {
        if (n == 0) return Rational(0);
        return harmonicRange(1, n);
    }

    //==========================================================================
    // FLINT-25: Dedekind Sum
    //==========================================================================
    // s(h,k) = Σ_{i=1}^{k-1} ((i/k))((hi/k))
    // ((x)) = x - floor(x) - 1/2 (x ∉ Z), ((x)) = 0 (x ∈ Z)
    //
    // Fast computation via the reciprocity law:
    // s(h,k) + s(k,h) = (h/k + k/h + 1/(hk)) / 12 - 1/4
    // Reduced to O(log(max(h,k))) via continued-fraction expansion

    Rational dedekindSum(const Int& h, const Int& k) {
        if (k <= 1) return Rational(0);
        Int hh = h % k;
        if (hh.isNegative()) hh += k;
        if (hh.isZero()) return Rational(0);

        // Recursion via the reciprocity law:
        // 12·k·s(h,k) = h(k-1)(2k-1) - k·Σ floor(ih/k) ... direct computation is O(k)
        // Reduced to O(log k) via continued-fraction expansion:
        // s(h,k) can be computed from the continued fraction [a_0; a_1, ..., a_n] of h/k
        //
        // Compute by the direct definition (used when k is small):
        if (k.bitLength() <= 20) {
            unsigned int kk = (unsigned int)k.toUInt64();
            unsigned int hmod = (unsigned int)hh.toUInt64();
            Rational sum(0);
            for (unsigned int i = 1; i < kk; i++) {
                // ((i/k)) = i/k - 1/2 (i/k is not an integer)
                // ((hi/k))
                unsigned int hi_mod_k = (unsigned int)((uint64_t)i * hmod % kk);
                if (hi_mod_k == 0) continue;  // ((integer)) = 0
                Rational saw_i = Rational(Int((int)i), Int((int)kk)) - Rational(1, 2);
                Rational saw_hi = Rational(Int((int)hi_mod_k), Int((int)kk)) - Rational(1, 2);
                sum += saw_i * saw_hi;
            }
            return sum;
        }

        // Large k: recursion via the reciprocity law
        // s(h,k) = (h²+k²+1)/(12hk) - 1/4 - s(k, h mod k) ... (when h,k are coprime)
        // First take the gcd
        Int g = IntGCD::gcd(hh, k);
        if (!g.isOne()) {
            // s(h,k) = s(h/g, k/g) (generalization for gcd > 1)
            return dedekindSum(hh / g, k / g);
        }

        // s(h,k) + s(k,h) = (h/k + k/h + 1/(hk))/12 - 1/4
        Rational lhs = (Rational(hh, k) + Rational(k, hh) +
                        Rational(Int(1), hh * k)) / Rational(12) -
                       Rational(1, 4);
        Rational sk_hmod = dedekindSum(k % hh, hh);
        return lhs - sk_hmod;
    }

    //==========================================================================
    // FLINT-26: Rational Enumeration
    //==========================================================================

    // Next positive rational in the Calkin-Wilf sequence
    // x = p/q → next is 1/(2·floor(p/q) + 1 - p/q) = q/(2·floor(p/q)·q - p + q)
    // Concisely: p/q → q/(2·floor(p/q)·q + q - p)
    Rational nextCalkinWilf(const Rational& x) {
        Int p = x.numerator(), q = x.denominator();
        if (p <= 0 || q <= 0) return Rational(1);
        Int fl = p / q;  // floor(p/q)
        return Rational(q, 2 * fl * q + q - p);
    }

    // Next positive rational in Stern-Brocot order (ascending by height)
    // Height = max(|p|, q). Within the same height, ascending order.
    Rational nextMinimal(const Rational& x) {
        Int p = x.numerator(), q = x.denominator();
        if (p <= 0 || q <= 0) return Rational(1);

        Int h = height(x);

        // Next of the same height: search for the next after p/q
        // Rationals of height h: max(p, q) = h
        // If q = h: increment p by 1 (p < h, gcd(p+1, h) == 1)
        // If p = h: increment q by 1 (q < h)

        if (q == h) {
            // denominator = h: increase the numerator
            for (Int np = p + 1; np < h; np += 1) {
                if (IntGCD::gcd(np, h).isOne())
                    return Rational(np, h);
            }
            // last of q = h → p = h, q = 1 (i.e. h/1 = h)
            // → next after p/q = (h-1)/h: move to the numerator=h series
            for (Int nq(1); nq <= h; nq += 1) {
                if (Rational(h, nq) > x && IntGCD::gcd(h, nq).isOne())
                    return Rational(h, nq);
            }
        } else if (p == h) {
            // numerator = h: increase the denominator → the value gets smaller
            for (Int nq = q + 1; nq < h; nq += 1) {
                if (IntGCD::gcd(h, nq).isOne())
                    return Rational(h, nq);
            }
        }

        // Move to the next height
        Int nh = h + 1;
        // Smallest rational of height nh: 1/nh
        return Rational(1, nh);
    }

    //==========================================================================
    // FLINT-27: Farey neighbors and minimal-height rationals
    //==========================================================================

    // Left and right neighbors of x = p/q in the Farey sequence F_Q
    // Left neighbor a/b: bp - aq = 1, b ≤ Q (choose the largest b)
    // Right neighbor c/d: cp - dq = 1 ... no, pc - qd = 1
    std::pair<Rational, Rational> fareyNeighbors(const Rational& x, const Int& Q) {
        Int p = x.numerator(), q = x.denominator();
        if (p.isNegative() || q <= 0 || q > Q) return {Rational(0), Rational(1)};

        // Left neighbor: a/b < p/q with bp - aq = 1, b ≤ Q
        // Find p·s + q·t = 1 by extended GCD → a = -t, b = s (adjustment needed)
        Int s, t;
        IntGCD::extendedGcd(p, q, s, t);
        // s·p + t·q = 1 → not (-t)·p - (-s)·q = ... but rather
        // a/b = (s·p - 1) / (s·q)... no, directly:
        // Left neighbor a/b: b·p - a·q = 1
        // b = s mod q (adjusted to positive), a = (b·p - 1) / q
        Int b_left = s;
        if (b_left.isNegative()) b_left += q;
        // Adjust so b ≤ Q: with k = floor((Q - b_left) / q), b_left += k*q
        if (b_left > 0) {
            Int k = (Q - b_left) / q;
            if (k > 0) b_left += k * q;
        }
        if (b_left <= 0) b_left = q - b_left;  // fallback
        // Largest b ≤ Q: maximum via q steps of b_left
        while (b_left + q <= Q) b_left += q;
        Int a_left = (b_left * p - 1) / q;
        Rational left(a_left, b_left);

        // Right neighbor: c/d > p/q with q·c - p·d = 1, d ≤ Q
        Int d_right = (-s);
        if (d_right.isNegative()) d_right += q;
        while (d_right + q <= Q) d_right += q;
        if (d_right <= 0) d_right += q;
        Int c_right = (p * d_right + 1) / q;
        Rational right(c_right, d_right);

        return {left, right};
    }

    // Minimal-height rational in the open interval (a, b) (Stern-Brocot tree search)
    Rational simplestBetween(const Rational& lo, const Rational& hi) {
        if (lo >= hi) return Rational();

        // Use the Stern-Brocot tree: start from 0/1 and 1/0
        // Midpoint = mediant. Move left or right by comparing with lo, hi.
        Int lp(0), lq(1);  // left boundary 0/1
        Int rp(1), rq(0);  // right boundary 1/0

        for (int iter = 0; iter < 10000; iter++) {
            Int mp = lp + rp, mq = lq + rq;
            Rational m(mp, mq);

            if (m <= lo) {
                lp = mp; lq = mq;
            } else if (m >= hi) {
                rp = mp; rq = mq;
            } else {
                return m;  // lo < m < hi
            }
        }
        return mediant(lo, hi);
    }

    //==========================================================================
    // FLINT-29: mod of a rational
    //==========================================================================
    // (p/q) mod m = p * q^{-1} mod m

    Int mod(const Rational& x, const Int& m) {
        Int p = x.numerator() % m;
        if (p.isNegative()) p += m;
        Int q = x.denominator();
        // q^{-1} mod m (returns 0 if q and m are not coprime)
        Int g = IntGCD::gcd(q, m);
        if (!g.isOne()) return Int(0);  // no inverse
        Int s, t;
        IntGCD::extendedGcd(q, m, s, t);
        Int q_inv = s % m;
        if (q_inv.isNegative()) q_inv += m;
        return (p * q_inv) % m;
    }

    //==========================================================================
    // FLINT-30: large-exponent power pow(Rational, Int)
    //==========================================================================

    Rational pow(const Rational& x, const Int& exponent) {
        if (exponent.isZero()) return Rational(1);

        bool neg_exp = exponent.isNegative();
        Int absExp = neg_exp ? -exponent : exponent;

        if (x.isZero()) {
            if (neg_exp) return Rational(Int::NaN(), Int(1));
            return Rational(0);
        }

        // Binary exponentiation
        Rational base = x;
        Rational result(1);
        while (!absExp.isZero()) {
            if (absExp.getBit(0))
                result *= base;
            base *= base;
            absExp >>= 1;
        }

        if (neg_exp) return inv(result);
        return result;
    }

} // namespace sangi
