// 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 #include #include #include #include #include #include // std::gcd #include // std::abs, std::max #include 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(num.w, nn), num_sign); sangi::Int d_i = sangi::Int::fromRawWordsPreNormalized( std::span(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(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(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(e); denominator_ = Int(1); } else { // mantissa / 2^(-e) → rational numerator_ = m; denominator_ = Int(1) << static_cast(-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(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(bits & exponentMask); bits >>= 11; int sign = static_cast(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(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(bits & exponentMask); bits >>= 8; int sign = static_cast(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(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(mantissa >> common)); denominator_ = Int(1) << (-e - common); } else { numerator_ = Int(static_cast(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(toDouble()); } Rational::operator int() const { return static_cast(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(lhs.numerator_.bitLength()) + static_cast(rhs.denominator_.bitLength()); int64_t cb_bits = static_cast(rhs.numerator_.bitLength()) + static_cast(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(nw, nn), num_sign); Int den_i = Int::fromRawWordsPreNormalized(std::span(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(nw, nn), num_sign); Int den_i = Int::fromRawWordsPreNormalized(std::span(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 ns(nw, num_hi ? 2 : 1); std::span 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 ns(nw, num_hi ? 2 : 1); std::span 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(shift - den_tz); Int new_num = lhs.numerator_ << num_shift; Int new_den = (den_tz > 0) ? (lhs.denominator_ >> static_cast(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(shift - num_tz); Int new_num = (num_tz > 0) ? (lhs.numerator_ >> static_cast(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(exponent)); result.denominator_ = sangi::pow(denominator_, static_cast(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(-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(num.bitLength()); int denBits = static_cast(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(shift), isNegative()); result.setPrecision(precision); return result; } //========================================================================== // Continued fractions //========================================================================== std::vector Rational::toContinuedFraction() const { std::vector 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& 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 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