// Copyright (C) 2026 Kiyotsugu Arai // SPDX-License-Identifier: LGPL-3.0-or-later // Rational.cpp // 多倍長有理数クラスの実装 // // 移植元: lib++20 有理数 (Allisone::有理数) // 改良点: // - std::bit_cast による IEEE 754 分解 (C++23) // - 非正規化数の指数計算バグ修正 // - 分母=0 チェック追加 // - move コンストラクタで std::move 使用 // - toDecimal の四捨五入バグ修正 #include #include #include #include #include #include #include // std::gcd #include // std::abs, std::max #include // _umul128 namespace { // Double-word 高速パス用ヘルパー // 分子・分母がともに 1 limb に収まるかチェック inline bool isSmallRational(const calx::Int& num, const calx::Int& den) { return num.size() <= 1 && den.size() <= 1 && num.isNormal() && den.isNormal(); } // 分子・分母がともに 2 limb 以下 (かつ 1 limb ではない) かチェック inline bool isMediumRationalPair(const calx::Rational&, const calx::Rational&) { return false; // 2-limb パスを無効化 (退行調査) } // Int から符号なし 64-bit 値を取得 (0 の場合は 0) inline uint64_t toU64(const calx::Int& x) { return x.isZero() ? 0 : x.word(0); } // ========================================================================= // 128-bit unsigned 演算ヘルパー (2-limb fast path 用) // ========================================================================= struct u128 { uint64_t lo, hi; bool isZero() const { return lo == 0 && hi == 0; } }; inline u128 u128_from_int(const calx::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 → 下位 256 bit のうち下位 128 bit のみ (gcd 用には十分) // ただし Rational の乗算では 256 bit が必要なので別関数 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) { // (a.hi:a.lo) × (b.hi:b.lo) の完全 256-bit 積 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 → 加算先 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 → 加算先 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 → 加算先 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 余り (安全版) inline uint64_t u128_mod64(u128 a, uint64_t d) { // (a.hi:a.lo) mod d // _udiv128 の前提: hi < d uint64_t r = a.hi % d; _udiv128(r, a.lo, d, &r); return r; } // 128-bit GCD (ユークリッド互除法) inline u128 u128_gcd(u128 a, u128 b) { // 高速パス: 両方が 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 なので a mod b = a → swap r = a; } else { // 両方 128-bit: a mod b を反復減算で計算 // (GCD の再帰深度は O(log(min(a,b))) なので両方 128-bit は稀、 // かつ最初の 1-2 反復で片方が 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 除算 → 128-bit 商 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 → 概算不要、GCD 用なので余りだけ必要な場合は別途 // 256-bit から u128 への mod: 約分用 inline u128 u256_mod128(u256 a, u128 m) { // 上位から順に mod を取る // a = (w3:w2:w1:w0), m = 128-bit // (w3:w2) mod m → r2, (r2:w1) mod m → r1, (r1:w0) mod m → result // 各ステップは 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); // ← これは gcd で mod ではない。正しく実装 // 実際は mod が必要だが、Rational の約分は gcd(num, g) で g は元の gcd なので // g は常に 64-bit に収まる。128-bit mod 128-bit は不要。 } // 実用上、Rational の約分で必要な gcd は gcd(numerator, g) で g ≤ 64-bit // なので u256_mod128 は使わない。代わに u256_mod64 を実装。 return {0, 0}; // placeholder } // 256-bit mod 64-bit → 64-bit 余り inline uint64_t u256_mod64(u256 a, uint64_t d) { // (w3:w2:w1:w0) mod d を上位から計算 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 商 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 加算 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 減算 (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 比較 (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]; } // u256 → Rational 構築 (最大 4 limb) // public コンストラクタ Rational(Int, Int, bool reduce) を使用 inline calx::Rational makeRationalFromU256(u256 num, int num_sign, u256 den) { if (num.isZero()) return calx::Rational(); size_t nn = num.wordCount(), dn = den.wordCount(); calx::Int n_i = calx::Int::fromRawWordsPreNormalized( std::span(num.w, nn), num_sign); calx::Int d_i = calx::Int::fromRawWordsPreNormalized( std::span(den.w, dn), 1); return calx::Rational(std::move(n_i), std::move(d_i), false); } // 符号付きで Rational を構築 (1 limb の分子・分母から) inline calx::Rational makeSmallRational(uint64_t num, int num_sign, uint64_t den) { if (num == 0) return calx::Rational(); calx::Int n(num); if (num_sign < 0) n = -n; return calx::Rational(std::move(n), calx::Int(den), false); } } // anonymous namespace namespace calx { //========================================================================== // コンストラクタ //========================================================================== 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_)) { } // 文字列コンストラクタ // 対応形式: "3/7", "42", "-1/3", "1.5" (10進のみ小数対応) Rational::Rational(std::string_view str, int base) : numerator_(Int(0)) , denominator_(Int(1)) { if (str.empty()) { numerator_ = Int::NaN(); return; } // "/" を探す → 分数形式 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; } // "." を探す → 小数形式 (10進のみ) 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 = calx::pow(Int(10), static_cast(n)); numerator_ = intVal * den + decVal; denominator_ = den; if (negative) numerator_ = -numerator_; reduce(); return; } // 整数形式 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; } // 多倍長浮動小数点数 → 有理数（正確変換） // Float = mantissa * 2^exponent (mantissa は正、符号は is_negative_ で管理) Rational::Rational(const Float& value) : numerator_(Int(0)) , denominator_(Int(1)) { // 特殊状態 if (value.isNaN() || value.isInfinity()) { numerator_ = Int::NaN(); return; } if (value.mantissa().isZero()) { return; // 0/1 } Int m = value.mantissa(); // 常に正 int64_t e = value.exponent(); if (e >= 0) { // mantissa * 2^e → 整数 numerator_ = m << static_cast(e); denominator_ = Int(1); } else { // mantissa / 2^(-e) → 有理数 numerator_ = m; denominator_ = Int(1) << static_cast(-e); reduce(); } if (value.isNegative()) { numerator_ = -numerator_; } } //========================================================================== // 代入演算子 //========================================================================== 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 倍精度浮動小数点数 → 有理数 // 【出典】http://ja.wikipedia.org/wiki/倍精度 // 【改良】std::bit_cast 使用 (C++23), 非正規化数のバグ修正 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 または NaN if (exponent == 2047) { numerator_ = Int::NaN(); denominator_ = Int(1); return *this; } int e; if (exponent != 0) { // 正規化数: 隠れビットを付加 mantissa |= 1ULL << 52; e = exponent - 1023 - 52; } else { // 非正規化数: 実効指数は 1 - 1023 (lib++20 のバグを修正) 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) { // mantissa と 2^(-e) の GCD は 2 の冪のみなので直接計算 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 単精度浮動小数点数 → 有理数 // 【出典】http://ja.wikipedia.org/wiki/単精度 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 または 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) { // mantissa と 2^(-e) の GCD は 2 の冪のみなので直接計算 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; } //========================================================================== // 単項演算子 //========================================================================== Rational Rational::operator-() const& { Rational result; result.numerator_ = -numerator_; result.denominator_ = denominator_; return result; } Rational Rational::operator-() && { numerator_ = -numerator_; return std::move(*this); } //========================================================================== // 約分 //========================================================================== bool Rational::reduce() { // 特殊状態チェック if (numerator_.isSpecialState() || denominator_.isSpecialState()) { return false; } // 分母=0 → NaN if (denominator_.isZero()) { numerator_ = Int::NaN(); denominator_ = Int(1); return true; } // 分母が正になるように符号を調整 if (denominator_.isNegative()) { numerator_ = -numerator_; denominator_ = -denominator_; } // 分子=0 → 0/1 に正規化 if (numerator_.isZero()) { numerator_.setSign(0); // 負のゼロを正規化 denominator_ = Int(1); return false; } // GCD で約分 (gcd() は内部で絶対値を扱うので abs() は不要) Int g = gcd(numerator_, denominator_); if (!g.isOne()) { numerator_ /= g; denominator_ /= g; return true; } return false; } //========================================================================== // 変換 //========================================================================== 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(); } //========================================================================== // 文字列変換 //========================================================================== 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); } // 小数点を含む文字列に変換する // 【移植元】lib++20 有理数::ToNumeric // 【修正】四捨五入の繰り上がりバグ (t[nn] → s[nn]) std::string Rational::toDecimal(int digits) const { if (isNaN()) return "NaN"; if (isZero()) return "0"; if (digits <= 0) digits = 0; // 整数部 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; } // 小数部を1桁ずつ計算 std::string decPart; Int ten(10); for (int i = 0; i <= digits; i++) { // 1桁余分に計算 (四捨五入用) r = r * ten; Int digit = r / denominator_; r = r % denominator_; decPart += digit.toString(); } // 四捨五入 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) { // 整数部に繰り上がり q += 1; intPart = q.toString(); } } return sign + intPart + "." + decPart; } //========================================================================== // ストリーム出力 //========================================================================== 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; } //========================================================================== // 比較演算子 //========================================================================== std::partial_ordering operator<=>(const Rational& lhs, const Rational& rhs) { if (lhs.isNaN() || rhs.isNaN()) { return std::partial_ordering::unordered; } // 符号による早期判定 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; } // 両方ゼロ if (lsign == 0) return std::partial_ordering::equivalent; // 同符号: サイズによる早期判定 (GMP mpq_cmp と同じ手法) // a/b vs c/d → a*d vs c*b // サイズ (limb 数) で乗算結果の大きさを近似 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) { // 正の場合: ad_size > cb_size → a*d > c*b → lhs > rhs // 負の場合: 逆 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 ファストパス: Int 構築を回避して raw 128-bit 比較 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); // a*d vs c*b を 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 ファストパス: 128x128 → 256-bit 比較 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 } // 浮動小数点近似による早期判定 // 上位数 limb から double 近似を計算し、十分離れていれば乗算を回避 { 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; // 相対差が十分大きければ確定 (double の精度: ~15桁) 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; } } } // 一般ケース: クロス乗算で正確に比較 (3引数版でバッファ再利用) 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; // 既約分数同士の比較: 分母子がそれぞれ等しいかチェック 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 // 符号による早期判定 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; // サイズによる早期判定 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 ファストパス 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 は 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; } // 一般ケース 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; // 既約分数: 分母が 1 かつ分子が等しい 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); } //========================================================================== // 四則演算 (Rational <-> Rational) // 加算: 大澤川の方法 (GCDで中間桁数削減) // 【出典】http://www.oishi.info.waseda.ac.jp/~samukawa/LongintRational.pdf p47-48 //========================================================================== Rational operator+(const Rational& lhs, const Rational& rhs) { // NaN 伝播 if (lhs.isNaN() || rhs.isNaN()) [[unlikely]] { Rational r; r.numerator_ = Int::NaN(); r.denominator_ = Int(1); return r; } // Double-word 高速パス: 分子・分母がすべて 1 limb かつ // 中間値 (a*d_g, c*b_g) も 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 // a * (d/g) と c * (b/g) を 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); // 128-bit 加減算で分子を計算 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 // (1 limb 入力なので num ≤ 2 * (2^64-1)^2 / g < 2^129, 実際はほぼ起きない) // num_hi の最上位ビットが立つ可能性はあるが 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(); // 分母 = (b/g) * d (128-bit) uint64_t den_hi; uint64_t den_lo = _umul128(b_g, d, &den_hi); // 約分: f = gcd(numerator, g) // g は 64-bit なので、numerator mod g を 128-bit で計算してから 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 なので _udiv128 の前提条件を満たす 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 除算 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) → 商を返す (hi < div が前提) num_lo = _udiv128(r_hi, num_lo, f, &rem); num_hi = q_hi; } // den = b_g * (d / f) — d/f は 64-bit, b_g は 64-bit den_lo = _umul128(b_g, d / f, &den_hi); } } // 結果を Rational として構築 (2 limb まで対応) if (num_hi == 0 && den_hi == 0) { return makeSmallRational(num_lo, num_sign, den_lo); } // 2 limb: Int を直接構築 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 高速パス: 分子・分母が 2 limb 以下 (20-38 桁) // 128-bit 演算で GCD/乗算を直接計算し、Int 構築オーバーヘッドを回避 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 → 簡略化: g は通常小さい) 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 { // g が 128-bit の場合、d/g と b/g は小さい → 一般パスにフォールバック goto general_add; } // a * d_g と c * b_g を 256-bit で計算 u256 ad = u256_mul128(abs_a, d_g); u256 cb = u256_mul128(abs_c, b_g); // 256-bit 加減算で分子を計算 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(); // 分母 = b_g * d (256-bit) u256 den256 = u256_mul128(b_g, d); // 約分: f = gcd(num256 mod g, g) — g は 64-bit なので 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_) { // 分母が同じ: (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()) { // 大澤川の方法 (3引数版で一時 Int 削減) 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: 結果は既約 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 伝播 if (lhs.isNaN() || rhs.isNaN()) [[unlikely]] { Rational r; r.numerator_ = Int::NaN(); r.denominator_ = Int(1); return r; } // Double-word 高速パス (加算と同構造、rhs の符号を反転) 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()); // 符号反転 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); // 128-bit 加減算で分子を計算 (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 高速パス (operator+ と同構造、rhs 符号反転) 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()); // 符号反転 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; } 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()) { // 大澤川の方法 (3引数版で一時 Int 削減) — 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: 結果は既約 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 伝播 if (lhs.isNaN() || rhs.isNaN()) [[unlikely]] { Rational r; r.numerator_ = Int::NaN(); r.denominator_ = Int(1); return r; } // Double-word 高速パス: 分子・分母がすべて 1 limb なら 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(); // 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 乗算 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); } // 2 limb に収まる場合: Int を直接構築 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 高速パス: 2 limb 以下の cross-GCD + 乗算 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); // g1, g2 が 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); } // 結果: num = a*c, den = b*d (最大 256-bit) u256 num256 = u256_mul128(a, c); u256 den256 = u256_mul128(b, d); return makeRationalFromU256(num256, sign, den256); } } } // GCD最適化: 交差約分してから乗算 // gcd() は内部で絶対値を扱うので abs() は不要 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); } // 符号を調整 (分母が負になる場合) 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 伝播 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 高速パス 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) と 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 高速パス: 2 limb 以下 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); } } } // 逆数を掛ける (交差約分最適化) // gcd() は内部で絶対値を扱うので abs() は不要 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 = -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) { // a % b = a - floor(a/b) * b Int m = floor(lhs / rhs); return lhs - Rational(m) * rhs; } //========================================================================== // 四則演算 (Rational <-> Int) //========================================================================== Rational operator+(const Rational& lhs, const Int& rhs) { // a/b + n = (a + n*b) / b (約分不要: gcd(a,b)==1 かつ gcd(n*b,b)==b) return Rational(lhs.numerator_ + rhs * lhs.denominator_, lhs.denominator_, false); } Rational operator+(const Int& lhs, const Rational& rhs) { return rhs + lhs; } Rational operator-(const Rational& lhs, const Int& rhs) { return Rational(lhs.numerator_ - rhs * lhs.denominator_, lhs.denominator_, false); } Rational operator-(const Int& lhs, const Rational& rhs) { return Rational(lhs * rhs.denominator_ - rhs.numerator_, rhs.denominator_, false); } Rational operator*(const Rational& lhs, const Int& rhs) { // 2のべき乗 → 分子を左シフト or 分母を右シフト (GCD 計算を回避) if (!rhs.isZero() && !rhs.isNegative() && rhs.bitLength() <= 64) { uint64_t w = rhs.data()[0]; if ((w & (w - 1)) == 0) { unsigned long shift; _BitScanForward64(&shift, w); // doubled() と同じロジックを一般化: 分母の下位ビットが 0 なら右シフト size_t den_tz = lhs.denominator_.countTrailingZeros(); if (den_tz >= shift) { return Rational(lhs.numerator_, lhs.denominator_ >> shift, false); } else { // 分母を den_tz だけ右シフト、残りを分子に左シフト 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: gcd(n, b) で事前約分 Int g = gcd(rhs, lhs.denominator_); if (g.isOne()) { return Rational(lhs.numerator_ * rhs, lhs.denominator_, false); } else { return Rational(lhs.numerator_ * (rhs / g), lhs.denominator_ / g, 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; } // 2のべき乗 → 分子を右シフト or 分母を左シフト (GCD 計算を回避) if (!rhs.isNegative() && rhs.bitLength() <= 64) { uint64_t w = rhs.data()[0]; if ((w & (w - 1)) == 0) { unsigned long shift; _BitScanForward64(&shift, w); // 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: gcd(a, n) で事前約分 Int g = gcd(lhs.numerator_, rhs); if (g.isOne()) { Rational result; result.numerator_ = lhs.numerator_; result.denominator_ = lhs.denominator_ * rhs; if (result.denominator_.isNegative()) { result.numerator_ = -result.numerator_; result.denominator_ = -result.denominator_; } return result; } else { Rational result; result.numerator_ = lhs.numerator_ / g; result.denominator_ = lhs.denominator_ * (rhs / g); if (result.denominator_.isNegative()) { result.numerator_ = -result.numerator_; result.denominator_ = -result.denominator_; } 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: gcd(n, a) で事前約分 Int g = gcd(lhs, rhs.numerator_); if (g.isOne()) { Rational result; result.numerator_ = lhs * rhs.denominator_; result.denominator_ = rhs.numerator_; if (result.denominator_.isNegative()) { result.numerator_ = -result.numerator_; result.denominator_ = -result.denominator_; } return result; } else { Rational result; result.numerator_ = (lhs / g) * rhs.denominator_; result.denominator_ = rhs.numerator_ / g; if (result.denominator_.isNegative()) { result.numerator_ = -result.numerator_; result.denominator_ = -result.denominator_; } return result; } } //========================================================================== // 四則演算 (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; } //========================================================================== // 複合代入 //========================================================================== 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: gcd(n, b) で事前約分 Int g = gcd(rhs, denominator_); if (g.isOne()) { numerator_ *= rhs; } else { numerator_ *= (rhs / g); denominator_ /= g; } return *this; } Rational& Rational::operator/=(const Int& rhs) { if (rhs.isZero()) { numerator_ = Int::NaN(); denominator_ = Int(1); return *this; } // (a/b) / n: gcd(a, n) で事前約分 Int g = gcd(numerator_, rhs); if (g.isOne()) { denominator_ *= rhs; } else { numerator_ /= g; denominator_ *= (rhs / g); } if (denominator_.isNegative()) { numerator_ = -numerator_; denominator_ = -denominator_; } 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); } //========================================================================== // インクリメント / デクリメント //========================================================================== Rational& Rational::operator++() { // a/b + 1 = (a+b)/b (約分不要: gcd(a,b)==1 なら 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; } //========================================================================== // 冪乗 //========================================================================== Rational Rational::pow(int exponent) const { Rational result; if (exponent > 0) { result.numerator_ = calx::pow(numerator_, static_cast(exponent)); result.denominator_ = calx::pow(denominator_, static_cast(exponent)); } else if (exponent < 0) { // 0 の負の冪は未定義 → NaN if (isZero()) { result.numerator_ = Int::NaN(); result.denominator_ = Int(1); return result; } unsigned int absExp = static_cast(-exponent); result.numerator_ = calx::pow(denominator_, absExp); result.denominator_ = calx::pow(numerator_, absExp); // 分母が負になった場合の符号調整 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; } //========================================================================== // シフト (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: フリー関数 //========================================================================== 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 (Wolfram 定義) // gcd(r1, r2, ...) = すべての ri/r が整数となる最大の有理数 r 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_); } //========================================================================== // ユーティリティ関数: inv, trunc, sgn, frac, mediant //========================================================================== Rational inv(const Rational& x) { if (x.isNaN()) return x; // Rational(den, num) のコンストラクタ内で: // num==0 → NaN (ゼロの逆数) // 符号調整 + 約分は reduce() が処理 return Rational(x.denominator(), x.numerator()); } Int trunc(const Rational& x) { // ゼロ方向への切り捨て = 整数除算 (C++ の除算はゼロ方向) return x.numerator() / x.denominator(); } Int round(const Rational& x) { // 最近接整数への丸め (half-away-from-zero: 0.5 → 1, -0.5 → -1) // |x| + 1/2 の floor を取り、元の符号を付与 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) { // 小数部分: x - floor(x) (常に 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()); } //========================================================================== // 便利関数 //========================================================================== Rational Rational::doubled() const { // x*2: 分母が偶数なら分母を半分、そうでなければ分子を2倍 if (!denominator_.getBit(0)) { // 偶数 (最下位ビット==0) return Rational(numerator_, denominator_ >> 1, false); } else { return Rational(numerator_ << 1, denominator_, false); } } Rational Rational::halved() const { // x/2: 分子が偶数なら分子を半分、そうでなければ分母を2倍 // getBit(0) は符号に依存しないので numerator_ に直接使える if (!numerator_.getBit(0)) { // 偶数 (最下位ビット==0) return Rational(numerator_ >> 1, denominator_, false); } else { return Rational(numerator_, denominator_ << 1, false); } } //========================================================================== // Rational → Float 変換（正確変換） //========================================================================== Float Rational::toMpFloat(int precision) const { if (isNaN()) return Float::nan(); // precision == 0 はデフォルト精度を使用 if (precision <= 0) { precision = Float::defaultPrecision(); } if (isZero()) return Float::zero(precision); // 分子を十分にシフトしてから除算し、precision ビット精度の仮数を得る Int num = abs(numerator_); Int den = denominator_; int numBits = static_cast(num.bitLength()); int denBits = static_cast(den.bitLength()); // precision + ガードビット分のビットを確保 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; } //========================================================================== // 連分数 //========================================================================== std::vector Rational::toContinuedFraction() const { std::vector result; if (isNaN()) return result; Int num = numerator_; // 負の場合あり Int den = denominator_; // 常に正 while (!den.isZero()) { // Floor division: q = floor(num / den) Int q = num / den; Int r = num - q * den; // C++ の truncation 剰余 // Floor に補正: r < 0 なら 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); // 後ろから計算: [a0; a1, a2, ..., an] // = a0 + 1/(a1 + 1/(a2 + ... + 1/an)) // 効率的な方法: 収束分数を前方反復で計算 // 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) //========================================================================== // 拡張ユークリッド互除法で a mod m から p/q を復元 // |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 的アプローチ: // r_0 = m, r_1 = a mod m // 互除法を |r_i| ≤ N になるまで実行 // その時点の (t_i が q, r_i が 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; } // q ≤ D のチェック if (q > D) return Rational(); // 復元失敗 return Rational(p, q); } //========================================================================== // FLINT-24: 調和数 (Harmonic Numbers) //========================================================================== // H_n = 1 + 1/2 + ... + 1/n // 分割統治法: H(lo, hi) で log(n) 深さの再帰、乗算回数を削減 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) // // 互恵律 (Reciprocity law) による高速計算: // s(h,k) + s(k,h) = (h/k + k/h + 1/(hk)) / 12 - 1/4 // 連分数展開で O(log(max(h,k))) に削減 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); // 互恵律による再帰: // 12·k·s(h,k) = h(k-1)(2k-1) - k·Σ floor(ih/k) ... 直接計算は O(k) // 連分数展開で O(log k) に削減: // s(h,k) は h/k の連分数 [a_0; a_1, ..., a_n] から計算可能 // // 直接定義で計算 (k が小さい場合に使用): 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 は整数でない) // ((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; } // 大きい k: 互恵律による再帰 // s(h,k) = (h²+k²+1)/(12hk) - 1/4 - s(k, h mod k) ... (h,k互いに素の場合) // まず gcd を取る Int g = IntGCD::gcd(hh, k); if (!g.isOne()) { // s(h,k) = s(h/g, k/g) (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) //========================================================================== // Calkin-Wilf 列で次の正有理数 // x = p/q → 次は 1/(2·floor(p/q) + 1 - p/q) = q/(2·floor(p/q)·q - p + q) // 簡潔に: 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); } // Stern-Brocot 順 (高さ昇順) で次の正有理数 // 高さ = max(|p|, q)。同じ高さ内では昇順。 Rational nextMinimal(const Rational& x) { Int p = x.numerator(), q = x.denominator(); if (p <= 0 || q <= 0) return Rational(1); Int h = height(x); // 同じ高さの次: p/q の次を探す // 高さ h の有理数: max(p, q) = h // q = h の場合: p を 1 増やす (p < h, gcd(p+1, h) == 1) // p = h の場合: q を 1 増やす (q < h) if (q == h) { // 分母 = h: 分子を増やす for (Int np = p + 1; np < h; np += 1) { if (IntGCD::gcd(np, h).isOne()) return Rational(np, h); } // q = h の最後 → p = h, q = 1 (つまり h/1 = h) // → p/q = (h-1)/h の次: 分子=h の系列へ 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) { // 分子 = h: 分母を増やす → 値は小さくなる for (Int nq = q + 1; nq < h; nq += 1) { if (IntGCD::gcd(h, nq).isOne()) return Rational(h, nq); } } // 次の高さへ Int nh = h + 1; // 高さ nh の最小有理数: 1/nh return Rational(1, nh); } //========================================================================== // FLINT-27: ファレイ近傍・最小高さ有理数 //========================================================================== // Farey 数列 F_Q における x = p/q の左右近傍 // 左近傍 a/b: bp - aq = 1, b ≤ Q (最大の b を選択) // 右近傍 c/d: cp - dq = 1 ... ではなく 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)}; // 左近傍: a/b < p/q で bp - aq = 1, b ≤ Q // 拡張GCDで p·s + q·t = 1 を求める → a = -t, b = s (調整が必要) Int s, t; IntGCD::extendedGcd(p, q, s, t); // s·p + t·q = 1 → (-t)·p - (-s)·q = ... ではなく // a/b = (s·p - 1) / (s·q)... いや、直接: // 左近傍 a/b: b·p - a·q = 1 // b = s mod q (正に調整), a = (b·p - 1) / q Int b_left = s; if (b_left.isNegative()) b_left += q; // b ≤ Q になるよう調整: 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 // 最大の b ≤ Q: b_left の q ステップで最大値 while (b_left + q <= Q) b_left += q; Int a_left = (b_left * p - 1) / q; Rational left(a_left, b_left); // 右近傍: c/d > p/q で 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}; } // 開区間 (a, b) 内の最小高さ有理数 (Stern-Brocot 木探索) Rational simplestBetween(const Rational& lo, const Rational& hi) { if (lo >= hi) return Rational(); // Stern-Brocot 木を使う: 0/1 と 1/0 から開始 // 中点 = mediant。lo, hi との比較で左右に進む。 Int lp(0), lq(1); // 左境界 0/1 Int rp(1), rq(0); // 右境界 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 //========================================================================== // (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 (q と m が互いに素でない場合は 0 を返す) Int g = IntGCD::gcd(q, m); if (!g.isOne()) return Int(0); // 逆元なし 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: 大指数冪乗 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); } // 二分冪乗 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 calx