// Copyright (C) 2026 Kiyotsugu Arai // SPDX-License-Identifier: LGPL-3.0-or-later // IntOps.cpp // Implementation of utilities for multi-precision integer operations #include #include #include #include #include #include #include #include // With C++20 , the __builtin_clzll compatibility macro is no longer needed // Use std::countl_zero / std::countr_zero namespace sangi { // Comparison operations bool IntOps::compareAbsLess(const Int& lhs, const Int& rhs) { // Handle special states if (lhs.isSpecialState() || rhs.isSpecialState()) { return IntSpecialStates::compareLessThan( lhs.getSign() < 0 ? -lhs : lhs, rhs.getSign() < 0 ? -rhs : rhs ); } // Compare by word count if (lhs.m_words.size() != rhs.m_words.size()) { return lhs.m_words.size() < rhs.m_words.size(); } // Same word count: compare from the most significant word for (int i = static_cast(lhs.m_words.size()) - 1; i >= 0; --i) { if (lhs.m_words[i] != rhs.m_words[i]) { return lhs.m_words[i] < rhs.m_words[i]; } } // Completely equal return false; } bool IntOps::compareAbsGreater(const Int& lhs, const Int& rhs) { // Handle special states if (lhs.isSpecialState() || rhs.isSpecialState()) { return IntSpecialStates::compareLessThan( rhs.getSign() < 0 ? -rhs : rhs, lhs.getSign() < 0 ? -lhs : lhs ); } // Compare by word count if (lhs.m_words.size() != rhs.m_words.size()) { return lhs.m_words.size() > rhs.m_words.size(); } // Same word count: compare from the most significant word for (int i = static_cast(lhs.m_words.size()) - 1; i >= 0; --i) { if (lhs.m_words[i] != rhs.m_words[i]) { return lhs.m_words[i] > rhs.m_words[i]; } } // Completely equal return false; } bool IntOps::compareAbsEqual(const Int& lhs, const Int& rhs) { // Handle special states if (lhs.isSpecialState() || rhs.isSpecialState()) { return lhs.getState() == rhs.getState(); } // Direct comparison of the word arrays return lhs.m_words == rhs.m_words; } // Addition operation void IntOps::addAbsolute(Int& result, const Int& other) { // Handle special states if (result.isSpecialState() || other.isSpecialState()) { result = IntSpecialStates::handleAddition(result, other); return; } // Special handling for addition with zero if (other.getSign() == 0) { return; } if (result.getSign() == 0) { result.m_words = other.m_words; result.setSign(other.getSign()); return; } // Addition processing size_t resultSize = result.m_words.size(); size_t otherSize = other.m_words.size(); size_t maxSize = std::max(resultSize, otherSize); // Expand the result vector size result.m_words.resize(maxSize, 0); uint64_t carry = mpn::add(result.m_words.data(), result.m_words.data(), maxSize, other.m_words.data(), otherSize); if (carry) { result.m_words.push_back(carry); } // normalize unnecessary: both inputs are already normalized → most significant word is non-zero } // The 3-argument addAbsolute/subtractAbsolute are defined inline in IntOps.hpp // Increment/decrement: use mpn::add_1/sub_1 to avoid constructing Int(1) void IntOps::addDelta(Int& value, int delta) { if (value.m_sign == 0) { value.m_words.resize_uninitialized(1); value.m_words[0] = 1; value.m_sign = delta; value.m_state = NumericState::Normal; return; } // Positive delta with positive value (or negative delta with negative value): add 1 to the magnitude if ((value.m_sign > 0) == (delta > 0)) { uint64_t carry = mpn::add_1(value.m_words.data(), value.m_words.size(), 1); if (carry) value.m_words.push_back(carry); } else { // Opposite signs: subtract 1 from the magnitude if (value.m_words.size() == 1 && value.m_words[0] == 1) { // |value| == 1 → result is 0 value.m_words.clear(); value.m_sign = 0; return; } mpn::sub_1(value.m_words.data(), value.m_words.size(), 1); if (value.m_words.back() == 0) value.m_words.pop_back(); } } // Single-word addition: result += word (signed) void IntOps::addWord(Int& result, uint64_t word) { if (word == 0) return; if (result.m_sign == 0) { result.m_words.resize_uninitialized(1); result.m_words[0] = word; result.m_sign = 1; return; } if (result.m_sign > 0) { uint64_t carry = mpn::add_1(result.m_words.data(), result.m_words.size(), word); if (carry) result.m_words.push_back(carry); } else { // result < 0: -(|result| - word) if (result.m_words.size() == 1) { uint64_t v = result.m_words[0]; if (v <= word) { uint64_t diff = word - v; if (diff == 0) { result.m_words.clear(); result.m_sign = 0; } else { result.m_words[0] = diff; result.m_sign = 1; } return; } } mpn::sub_1(result.m_words.data(), result.m_words.size(), word); if (result.m_words.back() == 0) result.m_words.pop_back(); if (result.m_words.empty()) result.m_sign = 0; } } // Single-word subtraction: result -= word (signed) void IntOps::subWord(Int& result, uint64_t word) { if (word == 0) return; if (result.m_sign == 0) { result.m_words.resize_uninitialized(1); result.m_words[0] = word; result.m_sign = -1; return; } if (result.m_sign < 0) { // result < 0: -(|result| + word) uint64_t carry = mpn::add_1(result.m_words.data(), result.m_words.size(), word); if (carry) result.m_words.push_back(carry); } else { // result > 0: |result| - word if (result.m_words.size() == 1) { uint64_t v = result.m_words[0]; if (v <= word) { uint64_t diff = word - v; if (diff == 0) { result.m_words.clear(); result.m_sign = 0; } else { result.m_words[0] = diff; result.m_sign = -1; } return; } } mpn::sub_1(result.m_words.data(), result.m_words.size(), word); if (result.m_words.back() == 0) result.m_words.pop_back(); if (result.m_words.empty()) result.m_sign = 0; } } // Single-word division: result /= word, returns the remainder uint64_t IntOps::divWord(Int& result, uint64_t word) { if (result.m_sign == 0) return 0; // Division by a power of two → right shift (faster than divmod_1) if (word != 0 && (word & (word - 1)) == 0) { unsigned long shift; _BitScanForward64(&shift, word); // The remainder is the low shift bits uint64_t rem = 0; if (shift > 0 && !result.m_words.empty()) { uint64_t mask = (uint64_t(1) << shift) - 1; rem = result.m_words[0] & mask; } rightShift(result, static_cast(shift)); return rem; } uint64_t rem = mpn::divmod_1(result.m_words.data(), result.m_words.data(), result.m_words.size(), word); while (!result.m_words.empty() && result.m_words.back() == 0) result.m_words.pop_back(); if (result.m_words.empty()) result.m_sign = 0; return rem; } // Single-word exact division: result /= word // Precondition: word > 0 and word evenly divides result (guaranteed by caller) // No quotient estimation needed (determined from the low end via Hensel inverse), so about 2x faster than divmod_1. void IntOps::divExactWord(Int& result, uint64_t word) { if (result.m_sign == 0) return; if (word == 1) return; // Separate the even factor: word = 2^tz × d_odd int tz = std::countr_zero(word); if (tz > 0) { rightShift(result, tz); if (result.m_sign == 0) return; // all bits dropped } uint64_t d_odd = word >> tz; if (d_odd == 1) return; // word was a power of two uint64_t inv = mpn::hensel_inverse_u64(d_odd); size_t new_size = mpn::divexact_by_odd( result.m_words.data(), result.m_words.data(), result.m_words.size(), d_odd, inv); result.m_words.resize(new_size); if (result.m_words.empty()) result.m_sign = 0; } // Subtraction operation (in-place: result -= |other|, precondition: |result| >= |other|) void IntOps::subtractAbsoluteInPlace(Int& result, const Int& other) { size_t resultSize = result.m_words.size(); size_t otherSize = other.m_words.size(); mpn::sub(result.m_words.data(), result.m_words.data(), resultSize, other.m_words.data(), otherSize); result.normalize(); } // Multiplication operation void IntOps::multiplyAbsolute(Int& result, const Int& other) { // Handle special states if (result.isSpecialState() || other.isSpecialState()) { result = IntSpecialStates::handleMultiplication(result, other); return; } // Special handling for multiplication by zero if (result.getSign() == 0 || other.getSign() == 0) { result.m_words.clear(); result.setSign(0); return; } // Special handling for multiplication by 1 if (other.m_words.size() == 1 && other.m_words[0] == 1) { // Sign is handled by the caller return; } if (result.m_words.size() == 1 && result.m_words[0] == 1) { result.m_words = other.m_words; // Sign is handled by the caller return; } size_t an = result.m_words.size(); size_t bn = other.m_words.size(); // Small-size fast path: use a stack buffer to avoid creating an intermediate Int // (eliminates the overhead of ScratchScope, arena alloc, fromRawWords) if (an < mpn::KARATSUBA_THRESHOLD && bn < mpn::KARATSUBA_THRESHOLD) { size_t rn = an + bn; uint64_t rbuf[mpn::KARATSUBA_THRESHOLD * 2]; mpn::multiply(rbuf, result.m_words.data(), an, other.m_words.data(), bn, nullptr); // The product has at most rn limbs; the leading 0 is at most the single most significant word if (rbuf[rn - 1] == 0) --rn; result.m_words.assign(rbuf, rbuf + rn); return; } // Large size: via IntMultiplication (arena + Karatsuba/Toom/NTT) generalMultiply(result, other, result); } // Multiplication operation (3-argument version: references lhs, rhs directly, no copy needed) // Precondition: zero and special states are already checked by the caller // Note: sign and state are set by the caller (reduces overhead) void IntOps::multiplyAbsolute(const Int& lhs, const Int& rhs, Int& result) { // Special handling for ×1 if (rhs.m_words.size() == 1 && rhs.m_words[0] == 1) { result.m_words = lhs.m_words; return; } if (lhs.m_words.size() == 1 && lhs.m_words[0] == 1) { result.m_words = rhs.m_words; return; } size_t an = lhs.m_words.size(); size_t bn = rhs.m_words.size(); // 1×1 specialization: perform 128-bit multiplication directly (skips the entire function call chain) if (an == 1 && bn == 1) { UInt128 prod = UInt128::multiply(lhs.m_words[0], rhs.m_words[0]); if (prod.high != 0) { result.m_words.resize_uninitialized(2); result.m_words[0] = prod.low; result.m_words[1] = prod.high; } else { result.m_words.resize_uninitialized(1); result.m_words[0] = prod.low; } return; } // N×1 specialization: call mul_1 directly (skips the multiply→mul_basecase dispatch) if (an == 1 || bn == 1) { const uint64_t* long_data; size_t long_n; uint64_t short_val; if (bn == 1) { long_data = lhs.m_words.data(); long_n = an; short_val = rhs.m_words[0]; } else { long_data = rhs.m_words.data(); long_n = bn; short_val = lhs.m_words[0]; } result.m_words.resize_uninitialized(long_n + 1); result.m_words[long_n] = mpn::mul_1(result.m_words.data(), long_data, long_n, short_val); if (result.m_words.back() == 0) result.m_words.pop_back(); return; } // Small-size: write directly into result.m_words (eliminates intermediate copy) if (an < mpn::KARATSUBA_THRESHOLD && bn < mpn::KARATSUBA_THRESHOLD) { size_t rn = an + bn; result.m_words.resize_uninitialized(rn); mpn::multiply(result.m_words.data(), lhs.m_words.data(), an, rhs.m_words.data(), bn, nullptr); if (result.m_words.back() == 0) result.m_words.pop_back(); return; } // Large-size: via IntMultiplication generalMultiply(lhs, rhs, result); } // Multiplication by a single word (checked) void IntOps::multiplyWord(Int& result, uint64_t word) { if (result.isSpecialState()) [[unlikely]] return; multiplyWordUnchecked(result, word); } // Multiplication by a single word (unchecked) void IntOps::multiplyWordUnchecked(Int& result, uint64_t word) { // Special handling for multiplication by zero if (result.getSign() == 0 || word == 0) { result.m_words.clear(); result.setSign(0); return; } // Special handling for multiplication by 1 if (word == 1) { return; } // Multiplication by a power of two → left shift (faster than mpn::mul_1) if ((word & (word - 1)) == 0) { unsigned long shift; _BitScanForward64(&shift, word); leftShift(result, static_cast(shift)); return; } // Single-word multiplication via mpn::mul_1 (supports MASM 4x unrolling) uint64_t carry = mpn::mul_1(result.m_words.data(), result.m_words.data(), result.m_words.size(), word); if (carry > 0) { result.m_words.push_back(carry); } result.normalize(); } // 3-argument mul: reuses result's buffer (equivalent to GMP mpz_mul) void IntOps::mulUnchecked(const Int& lhs, const Int& rhs, Int& result) { if (lhs.getSign() == 0 || rhs.getSign() == 0) [[unlikely]] { result.m_words.clear(); result.m_sign = 0; result.m_state = NumericState::Normal; return; } int resultSign = lhs.getSign() * rhs.getSign(); // BUGFIX (cas-c5cc deep #4): alias safety. multiplyAbsolute/squareUnchecked // resize `result.m_words` before/while reading the operand limbs, so if // `result` aliases `lhs` or `rhs` the operand is clobbered mid-computation, // producing a silently wrong product. Callers do alias deliberately to // save temporaries (e.g. Rational operator+/-: `mulUnchecked(num, t, t)`), // which corrupted Rational add/subtract whenever the denominators shared a // factor (general path). Route aliased calls through a fresh temporary. if (&result == &lhs || &result == &rhs) [[unlikely]] { Int tmp; if (&lhs == &rhs) squareUnchecked(lhs, tmp); else multiplyAbsolute(lhs, rhs, tmp); tmp.m_sign = resultSign; tmp.m_state = NumericState::Normal; result = std::move(tmp); return; } if (&lhs == &rhs) { squareUnchecked(lhs, result); } else { multiplyAbsolute(lhs, rhs, result); } result.m_sign = resultSign; result.m_state = NumericState::Normal; } void IntOps::mul(const Int& lhs, const Int& rhs, Int& result) { if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] { result = lhs * rhs; return; } mulUnchecked(lhs, rhs, result); } // Square computation (checked) void IntOps::square(const Int& value, Int& result) { if (value.isSpecialState()) [[unlikely]] { result = IntSpecialStates::handleMultiplication(value, value); return; } squareUnchecked(value, result); } // Square computation (unchecked) void IntOps::squareUnchecked(const Int& value, Int& result) { // Special handling for the square of zero if (value.getSign() == 0) { result.m_words.clear(); result.m_sign = 0; result.m_state = NumericState::Normal; return; } // BUGFIX (cas-c5cc deep #4): alias safety (see mulUnchecked). resize of // result.m_words would clobber `value` when result aliases value. if (&result == &value) [[unlikely]] { Int tmp; squareUnchecked(value, tmp); result = std::move(tmp); return; } size_t n = value.m_words.size(); // Small size: ASM mul_basecase is faster than ASM sqr_basecase constexpr size_t SQR_MUL_FALLBACK = 8; if (n <= SQR_MUL_FALLBACK) { multiplyAbsolute(value, value, result); result.m_sign = 1; result.m_state = NumericState::Normal; return; } // Large size: dedicated squaring via mpn::square (sped up by exploiting symmetry) result.m_words.resize_uninitialized(2 * n); size_t scratch_sz = mpn::square_scratch_size(n); if (scratch_sz > 0) { ScratchScope scope; uint64_t* scratch = getThreadArena().alloc_limbs(scratch_sz); mpn::square(result.m_words.data(), value.m_words.data(), n, scratch); } else { mpn::square(result.m_words.data(), value.m_words.data(), n, nullptr); } if (result.m_words.back() == 0) result.m_words.pop_back(); result.m_sign = 1; result.m_state = NumericState::Normal; } // 3-argument divUnchecked: version that omits the isSpecialState check void IntOps::divUnchecked(const Int& dividend, const Int& divisor, Int& result) { // Division by zero (not a special state, but kept for safety) if (divisor.m_sign == 0) [[unlikely]] { result = dividend / divisor; return; } // Dividend is zero if (dividend.m_sign == 0) { result.m_words.clear(); result.m_sign = 0; result.m_state = NumericState::Normal; return; } size_t an = dividend.m_words.size(); size_t bn = divisor.m_words.size(); // |dividend| < |divisor| → quotient = 0 int abscmp = mpn::cmp(dividend.m_words.data(), an, divisor.m_words.data(), bn); if (abscmp < 0) { result.m_words.clear(); result.m_sign = 0; result.m_state = NumericState::Normal; return; } if (abscmp == 0) { result.m_words.resize(1); result.m_words[0] = 1; result.m_sign = dividend.m_sign * divisor.m_sign; result.m_state = NumericState::Normal; return; } int resultSign = dividend.m_sign * divisor.m_sign; // 1-limb divisor: divmod_1 fast path if (bn == 1) { size_t qn = an; result.m_words.resize_uninitialized(qn); mpn::divmod_1(result.m_words.data(), dividend.m_words.data(), an, divisor.m_words[0]); while (qn > 0 && result.m_words[qn - 1] == 0) qn--; result.m_words.resize(qn); result.m_sign = (qn == 0) ? 0 : resultSign; result.m_state = NumericState::Normal; return; } // Multi-word: quotient-only division (remainder not needed) { size_t scratch_sz = mpn::divide_q_scratch_size(an, bn); size_t qn = an - bn + 1; // Write the quotient directly into the result buffer (eliminates copy) result.m_words.resize_uninitialized(qn + 1); uint64_t* qp = result.m_words.data(); ScratchScope scope; uint64_t* scratch = getThreadArena().alloc_limbs(scratch_sz); qn = mpn::divide_q(qp, dividend.m_words.data(), an, divisor.m_words.data(), bn, scratch); result.m_words.resize(qn); result.m_sign = (qn == 0) ? 0 : resultSign; result.m_state = NumericState::Normal; } } // 3-argument div: checked wrapper void IntOps::div(const Int& dividend, const Int& divisor, Int& result) { if (dividend.isSpecialState() || divisor.isSpecialState()) [[unlikely]] { result = dividend / divisor; return; } divUnchecked(dividend, divisor, result); } // ================================================================ // addmul / submul: rop ±= a * b // ================================================================ // rop += a * b void IntOps::addmul(Int& rop, const Int& a, const Int& b) { // Special states if (a.isSpecialState() || b.isSpecialState() || rop.isSpecialState()) [[unlikely]] { rop = rop + a * b; return; } // Do nothing if a or b is zero if (a.getSign() == 0 || b.getSign() == 0) return; // If rop is zero, simply rop = a * b if (rop.getSign() == 0) { mul(a, b, rop); return; } // Compute a * b into a temporary buffer Int product; mul(a, b, product); // rop += product (3-argument add reuses the buffer) add(rop, product, rop); } // rop -= a * b void IntOps::submul(Int& rop, const Int& a, const Int& b) { if (a.isSpecialState() || b.isSpecialState() || rop.isSpecialState()) [[unlikely]] { rop = rop - a * b; return; } if (a.getSign() == 0 || b.getSign() == 0) return; if (rop.getSign() == 0) { mul(a, b, rop); rop.m_sign = -rop.m_sign; return; } Int product; mul(a, b, product); sub(rop, product, rop); } // rop += a * word (single-word version) // ★ If signs are equal, accumulate directly via mpn::addmul_1 (one pass, no temporary buffer) void IntOps::addmul(Int& rop, const Int& a, uint64_t word) { if (a.isSpecialState() || rop.isSpecialState()) [[unlikely]] { rop = rop + a * Int(word); return; } if (a.getSign() == 0 || word == 0) return; if (rop.getSign() == 0) { rop = a; multiplyWord(rop, word); return; } size_t an = a.m_words.size(); int prodSign = a.getSign(); if (rop.getSign() == prodSign) { // ★ Same sign: accumulate directly into rop via mpn::addmul_1 size_t rn = rop.m_words.size(); size_t new_size = std::max(rn, an) + 1; rop.m_words.resize(new_size); // expand + zero-fill uint64_t* rp = rop.m_words.data(); // addmul_1: rp += a * word uint64_t carry = mpn::addmul_1(rp, a.m_words.data(), an, word); // Propagate the carry for (size_t i = an; i < new_size && carry; i++) { uint64_t sum = rp[i] + carry; carry = (sum < rp[i]) ? 1ULL : 0ULL; rp[i] = sum; } if (carry) { rop.m_words.resize(new_size + 1); rop.m_words.data()[new_size] = carry; new_size++; } // Normalize (remove leading zeros) while (new_size > 1 && rop.m_words.data()[new_size - 1] == 0) new_size--; rop.m_words.resize(new_size); return; } // Opposite signs: via a temporary buffer (in-place subtraction is complex) size_t pn = an + 1; constexpr size_t BUF_LIMIT = 128; uint64_t stack_buf[BUF_LIMIT]; uint64_t* pbuf = (pn <= BUF_LIMIT) ? stack_buf : new uint64_t[pn]; uint64_t carry = mpn::mul_1(pbuf, a.m_words.data(), an, word); if (carry) { pbuf[an] = carry; } else { pn = an; } Int product; product.m_words.resize_uninitialized(pn); std::copy_n(pbuf, pn, product.m_words.data()); product.m_sign = prodSign; product.m_state = NumericState::Normal; if (pbuf != stack_buf) delete[] pbuf; add(rop, product, rop); } // rop -= a * word (single-word version) // submul is the sign-flipped version of addmul: rop -= a*w == rop += (-a)*w void IntOps::submul(Int& rop, const Int& a, uint64_t word) { if (a.isSpecialState() || rop.isSpecialState()) [[unlikely]] { rop = rop - a * Int(word); return; } if (a.getSign() == 0 || word == 0) return; if (rop.getSign() == 0) { rop = a; multiplyWord(rop, word); rop.m_sign = -rop.m_sign; return; } size_t an = a.m_words.size(); int prodSign = -a.getSign(); // submul: sign flip if (rop.getSign() == prodSign) { // Same sign (rop and -a*w): accumulate directly via mpn::addmul_1 size_t rn = rop.m_words.size(); size_t new_size = std::max(rn, an) + 1; rop.m_words.resize(new_size); uint64_t* rp = rop.m_words.data(); uint64_t carry = mpn::addmul_1(rp, a.m_words.data(), an, word); for (size_t i = an; i < new_size && carry; i++) { uint64_t sum = rp[i] + carry; carry = (sum < rp[i]) ? 1ULL : 0ULL; rp[i] = sum; } if (carry) { rop.m_words.resize(new_size + 1); rop.m_words.data()[new_size] = carry; new_size++; } while (new_size > 1 && rop.m_words.data()[new_size - 1] == 0) new_size--; rop.m_words.resize(new_size); return; } // Opposite signs: via a temporary buffer size_t pn = an + 1; constexpr size_t BUF_LIMIT = 128; uint64_t stack_buf[BUF_LIMIT]; uint64_t* pbuf = (pn <= BUF_LIMIT) ? stack_buf : new uint64_t[pn]; uint64_t carry = mpn::mul_1(pbuf, a.m_words.data(), an, word); if (carry) { pbuf[an] = carry; } else { pn = an; } Int product; product.m_words.resize_uninitialized(pn); std::copy_n(pbuf, pn, product.m_words.data()); product.m_sign = prodSign; product.m_state = NumericState::Normal; if (pbuf != stack_buf) delete[] pbuf; add(rop, product, rop); } // Division operation (checked) void IntOps::divideAbsolute(Int& result, const Int& divisor) { if (result.isSpecialState() || divisor.isSpecialState()) [[unlikely]] { result = IntSpecialStates::handleDivision(result, divisor); return; } if (divisor.getSign() == 0) [[unlikely]] { result.setState(NumericState::NaN, NumericError::DivideByZero); return; } divideAbsoluteUnchecked(result, divisor); } // Division operation (unchecked): special states and divisor==0 are already verified by the caller void IntOps::divideAbsoluteUnchecked(Int& result, const Int& divisor) { // Division of zero if (result.getSign() == 0) { return; // 0 / x = 0 (x ≠ 0) } // Division by the same value if (compareAbsEqual(result, divisor)) { result.m_words.resize(1); result.m_words[0] = 1; // Sign is handled by the caller return; } // Case |result| < |divisor| if (compareAbsLess(result, divisor)) { result.m_words.clear(); // 0 result.setSign(0); return; } // Select the algorithm based on the global setting switch (g_division_algorithm) { case DivisionAlgorithm::Knuth: { // Use Knuth Algorithm D // divKnuth handles signed Int, so save the sign result.setSign(1); // treat as the magnitude Int abs_divisor = divisor; abs_divisor.setSign(1); Int remainder; Int quotient = IntDivision::divKnuth(result, abs_divisor, remainder); result = std::move(quotient); // Sign is handled by the caller (stays positive since this is absolute-value division) if (!result.isZero()) { result.setSign(1); } return; } case DivisionAlgorithm::BitByBit: // Jump directly to bitwise division goto bit_by_bit_division; case DivisionAlgorithm::SingleWordOnly: // Optimize only single-word divisors; otherwise bitwise if (divisor.m_words.size() == 1) { goto single_word_division; } else { goto bit_by_bit_division; } case DivisionAlgorithm::Auto: default: // Automatic selection (optimization based on benchmark results) // 1. Power-of-two detection → PowerOfTwo (28-168x speedup) { bool is_power_of_two = true; size_t power = 0; bool found_nonzero = false; // Check whether the divisor is a power of two for (size_t i = 0; i < divisor.m_words.size(); ++i) { if (divisor.m_words[i] != 0) { // More than one non-zero word means it is not a power of two if (found_nonzero) { is_power_of_two = false; break; } found_nonzero = true; // Check whether exactly one bit is set uint64_t w = divisor.m_words[i]; if ((w & (w - 1)) != 0) { is_power_of_two = false; break; } // Compute the bit position for (int b = 0; b < 64; ++b) { if (w & (1ULL << b)) { power = i * 64 + b; break; } } } } if (is_power_of_two && power > 0) { // Use PowerOfTwo Int remainder_tmp; Int quotient_tmp = IntDivision::divPowerOfTwo(result, power, remainder_tmp); result = std::move(quotient_tmp); return; } } // 2. Single-word divisor → SingleWordOnly (1.67x speedup) if (divisor.m_words.size() == 1) { goto single_word_division; } // 3. Multi-word divisor → mpn::divide (BZ / schoolbook) if (divisor.m_words.size() >= 2) { const uint64_t* a = result.m_words.data(); size_t an = result.m_words.size(); const uint64_t* b = divisor.m_words.data(); size_t bn = divisor.m_words.size(); // --- 2-word divisor × 2-word dividend specialization (128÷128 bit) --- // Avoid the overhead of lshift/div_basecase/assign and // perform Knuth 3÷2 division inline (quotient is at most 1 word) if (bn == 2 && an == 2) { unsigned shift = std::countl_zero(b[1]); uint64_t nb1, nb0, u2, u1, u0; if (shift > 0) { nb1 = (b[1] << shift) | (b[0] >> (64 - shift)); nb0 = b[0] << shift; u2 = a[1] >> (64 - shift); u1 = (a[1] << shift) | (a[0] >> (64 - shift)); u0 = a[0] << shift; } else { nb1 = b[1]; nb0 = b[0]; u2 = 0; u1 = a[1]; u0 = a[0]; } // u2 < nb1 always holds (u2 < 2^shift, nb1 >= 2^63) auto [q, re] = UInt128::divmod_fast(u2, u1, nb1); // Knuth adjustment (at most 2 times) auto prod = UInt128::multiply(q, nb0); if (prod.high > re || (prod.high == re && prod.low > u0)) { q--; re += nb1; if (re >= nb1) { // no overflow → re-check prod = UInt128::multiply(q, nb0); if (prod.high > re || (prod.high == re && prod.low > u0)) { q--; } } } if (q == 0) { result.m_words.clear(); result.setSign(0); } else { result.m_words.resize(1); result.m_words[0] = q; } return; } // --- Small-size fast-path: div_basecase directly with a stack buffer --- // Avoid the overhead of ScratchScope/Arena allocation constexpr size_t SMALL_DIV_THRESHOLD = 16; // 1024 bits if (an <= SMALL_DIV_THRESHOLD && bn < mpn::BZ_THRESHOLD) { uint64_t na[SMALL_DIV_THRESHOLD + 2]; uint64_t nb[SMALL_DIV_THRESHOLD + 1]; uint64_t qq[SMALL_DIV_THRESHOLD + 1]; unsigned shift = std::countl_zero(b[bn - 1]); if (shift > 0) { mpn::lshift(nb, b, bn, shift); na[an] = mpn::lshift(na, a, an, shift); } else { std::memcpy(nb, b, bn * sizeof(uint64_t)); std::memcpy(na, a, an * sizeof(uint64_t)); na[an] = 0; } size_t nan = an + (na[an] ? 1 : 0); na[nan] = 0; // sentinel mpn::div_basecase(qq, na, nan, nb, bn); size_t qn = mpn::normalized_size(qq, nan - bn + 1); if (qn == 0) { result.m_words.clear(); result.setSign(0); } else { result.m_words.assign(qq, qq + qn); result.setSign(1); } return; } // --- Normal path: Arena + mpn::divide (BZ / schoolbook) --- ScratchScope scope; size_t qn_max = an - bn + 1; uint64_t* q_buf = getThreadArena().alloc_limbs(qn_max + 1); size_t scratch_size = mpn::divide_scratch_size(an, bn); uint64_t* scratch = getThreadArena().alloc_limbs(scratch_size); uint64_t* r_buf = getThreadArena().alloc_limbs(bn); std::memset(q_buf, 0, (qn_max + 1) * sizeof(uint64_t)); size_t qn = mpn::divide(q_buf, r_buf, a, an, b, bn, scratch); if (qn == 0) { result = Int::Zero(); } else { result = Int::fromRawWords(std::span(q_buf, qn), 1); } return; } // 4. Default → BitByBit (single-word fallback) break; } // Optimization for a single-word divisor (in-place, no heap alloc) single_word_division: if (divisor.m_words.size() == 1) { uint64_t d = divisor.m_words[0]; uint64_t rem = 0; // Process from the most significant word — write the quotient directly at the read position for (int i = static_cast(result.m_words.size()) - 1; i >= 0; --i) { auto [q, r] = UInt128::divmod_fast(rem, result.m_words[i], d); result.m_words[i] = q; rem = r; } // normalize: the leading word becomes 0 in at most 1 word if (!result.m_words.empty() && result.m_words.back() == 0) { result.m_words.pop_back(); } if (result.m_words.empty()) { result.setSign(0); } return; } // Bitwise long division (fallback when DivisionAlgorithm::BitByBit is specified) // Note: Auto mode uses Knuth Algorithm D (bug fixed in Phase 0) bit_by_bit_division: Int remainder = Int::Zero(); Int quotient = Int::Zero(); // Difference in number of digits int bitDiff = static_cast((result.m_words.size() - divisor.m_words.size()) * 64); // Compute the most significant bit position of the most significant word uint64_t resultMsw = result.m_words.back(); uint64_t divisorMsw = divisor.m_words.back(); int resultMsBit = 63; if (resultMsw > 0) { // Find the most significant bit position of resultMsw for (int i = 63; i >= 0; --i) { if (resultMsw & (1ULL << i)) { resultMsBit = i; break; } } } int divisorMsBit = 63; if (divisorMsw > 0) { // Find the most significant bit position of divisorMsw for (int i = 63; i >= 0; --i) { if (divisorMsw & (1ULL << i)) { divisorMsBit = i; break; } } } bitDiff += resultMsBit - divisorMsBit; // Shift the divisor to align it Int shiftedDivisor = divisor; leftShift(shiftedDivisor, bitDiff); // Adjust if the aligned divisor became larger than the dividend if (compareAbsGreater(shiftedDivisor, result)) { rightShift(shiftedDivisor, 1); bitDiff--; } // Division processing remainder = result; quotient.m_words.resize((bitDiff / 64) + 1, 0); while (bitDiff >= 0) { if (compareAbsGreater(remainder, shiftedDivisor) || compareAbsEqual(remainder, shiftedDivisor)) { subtractAbsolute(remainder, shiftedDivisor, remainder); // Set the quotient bit quotient.m_words[bitDiff / 64] |= (1ULL << (bitDiff % 64)); } rightShift(shiftedDivisor, 1); bitDiff--; } result = std::move(quotient); result.normalize(); // quotient is built from Int::Zero(), so its sign stays 0. // If non-zero after normalization, set the sign (positive since this is absolute-value division). if (!result.m_words.empty()) { result.setSign(1); } } // Modulo operation (checked) void IntOps::moduloAbsolute(Int& result, const Int& divisor) { if (result.isSpecialState() || divisor.isSpecialState()) [[unlikely]] { result = IntSpecialStates::handleModulo(result, divisor); return; } if (divisor.getSign() == 0) [[unlikely]] { result.setState(NumericState::NaN, NumericError::DivideByZero); return; } moduloAbsoluteUnchecked(result, divisor); } // Modulo operation (unchecked): special states and divisor==0 are already verified by the caller void IntOps::moduloAbsoluteUnchecked(Int& result, const Int& divisor) { // Modulo of zero if (result.getSign() == 0) { return; // 0 % x = 0 (x ≠ 0) } // Modulo by the same value if (compareAbsEqual(result, divisor)) { result.m_words.clear(); result.setSign(0); return; } // Case |result| < |divisor| if (compareAbsLess(result, divisor)) { return; // magnitude is already smaller, so no change } // Optimization for a single-word divisor if (divisor.m_words.size() == 1) { uint64_t d = divisor.m_words[0]; uint64_t remainder = 0; // Process from the most significant word for (int i = static_cast(result.m_words.size()) - 1; i >= 0; --i) { std::pair qr = UInt128::divmod_fast(remainder, result.m_words[i], d); remainder = qr.second; } if (remainder == 0) { result.m_words.clear(); result.setSign(0); } else { result.m_words.resize(1); result.m_words[0] = remainder; // Sign is handled by the caller } return; } // Normal modulo processing divmod(result, divisor, result); } // Simultaneous computation of quotient and remainder Int IntOps::divmod(const Int& dividend, const Int& divisor, Int& remainder) { // Handle special states if (dividend.isSpecialState() || divisor.isSpecialState()) { remainder = IntSpecialStates::handleModulo(dividend, divisor); return IntSpecialStates::handleDivision(dividend, divisor); } // Division-by-zero check if (divisor.getSign() == 0) { remainder.setState(NumericState::NaN, NumericError::DivideByZero); Int quotient = Int::NaN(); quotient.setState(NumericState::NaN, NumericError::DivideByZero); return quotient; } // Division of zero if (dividend.getSign() == 0) { remainder = Int::Zero(); return Int::Zero(); } // Determine the sign int quotientSign = (dividend.getSign() == divisor.getSign()) ? 1 : -1; // Use the absolute values Int absDividend = dividend; if (absDividend.getSign() < 0) { absDividend.setSign(1); } Int absDivisor = divisor; if (absDivisor.getSign() < 0) { absDivisor.setSign(1); } // Case |dividend| < |divisor| if (compareAbsLess(absDividend, absDivisor)) { remainder = dividend; // magnitude is already smaller, so the remainder is the dividend itself return Int::Zero(); } // Optimization for a single-word divisor if (absDivisor.m_words.size() == 1) { uint64_t d = absDivisor.m_words[0]; uint64_t r = 0; // Process from the most significant word std::vector quotient(absDividend.m_words.size(), 0); for (int i = static_cast(absDividend.m_words.size()) - 1; i >= 0; --i) { UInt128 current(r, absDividend.m_words[i]); std::pair qr = UInt128::divmod_fast(r, absDividend.m_words[i], d); quotient[i] = qr.first; r = qr.second; } // Set the remainder if (r == 0) { remainder = Int::Zero(); } else { remainder.m_words.resize(1); remainder.m_words[0] = r; remainder.setSign(dividend.getSign()); // remainder's sign matches the dividend remainder.setState(NumericState::Normal); } // Set the quotient Int result; result.m_words.assign(quotient.data(), quotient.data() + quotient.size()); result.setSign(quotientSign); result.setState(NumericState::Normal); result.normalize(); return result; } // Normal division processing: mpn::divide (BZ / schoolbook) if (absDivisor.m_words.size() >= 2) { const uint64_t* a = absDividend.m_words.data(); size_t an = absDividend.m_words.size(); const uint64_t* b = absDivisor.m_words.data(); size_t bn = absDivisor.m_words.size(); ScratchScope scope; size_t qn_max = an - bn + 1; uint64_t* q_buf = getThreadArena().alloc_limbs(qn_max + 1); uint64_t* r_buf = getThreadArena().alloc_limbs(bn); size_t scratch_size = mpn::divide_scratch_size(an, bn); uint64_t* scratch = getThreadArena().alloc_limbs(scratch_size); std::memset(q_buf, 0, (qn_max + 1) * sizeof(uint64_t)); size_t qn = mpn::divide(q_buf, r_buf, a, an, b, bn, scratch); // Convert the quotient to Int Int quotient; if (qn == 0) { quotient = Int::Zero(); } else { quotient = Int::fromRawWords(std::span(q_buf, qn), quotientSign); } // Convert the remainder to Int size_t rn = mpn::normalized_size(r_buf, bn); if (rn == 0) { remainder = Int::Zero(); } else { remainder = Int::fromRawWords(std::span(r_buf, rn), dividend.getSign()); } return quotient; } // Fallback: BitByBit (should not be reached here since single-word divisors are handled above) Int quotient = Int::Zero(); remainder = absDividend; // Difference in number of digits int bitDiff = static_cast((absDividend.m_words.size() - absDivisor.m_words.size()) * 64); // Compute the most significant bit position of the most significant word uint64_t dividendMsw = absDividend.m_words.back(); uint64_t divisorMsw = absDivisor.m_words.back(); int dividendMsBit = 63; if (dividendMsw > 0) { // Find the most significant bit position of dividendMsw for (int i = 63; i >= 0; --i) { if (dividendMsw & (1ULL << i)) { dividendMsBit = i; break; } } } int divisorMsBit = 63; if (divisorMsw > 0) { // Find the most significant bit position of divisorMsw for (int i = 63; i >= 0; --i) { if (divisorMsw & (1ULL << i)) { divisorMsBit = i; break; } } } bitDiff += dividendMsBit - divisorMsBit; // Shift the divisor to align it Int shiftedDivisor = absDivisor; leftShift(shiftedDivisor, bitDiff); // Adjust if the aligned divisor became larger than the dividend if (compareAbsGreater(shiftedDivisor, absDividend)) { rightShift(shiftedDivisor, 1); bitDiff--; } // Division processing quotient.m_words.resize((bitDiff / 64) + 1, 0); while (bitDiff >= 0) { if (compareAbsGreater(remainder, shiftedDivisor) || compareAbsEqual(remainder, shiftedDivisor)) { subtractAbsolute(remainder, shiftedDivisor, remainder); // Set the quotient bit quotient.m_words[bitDiff / 64] |= (1ULL << (bitDiff % 64)); } rightShift(shiftedDivisor, 1); bitDiff--; } // Adjust the sign quotient.setSign(quotientSign); if (!remainder.isZero()) { remainder.setSign(dividend.getSign()); // remainder's sign matches the dividend } quotient.normalize(); remainder.normalize(); return quotient; } // ========================================================================= // Bit operations (two's-complement semantics) // // Bitwise operations on signed integers are interpreted as two's-complement representation. // The internal representation is sign-magnitude, so convert using the following identities: // // Positive number +m: bit pattern = m (extended infinitely with 0 in the high bits) // Negative number -m: bit pattern = ~(m-1) (extended infinitely with 1 in the high bits) // // Formula for each operation, per sign combination: // AND: (+a)&(+b) = +(a&b) // (+a)&(-b) = +(a & ~(|b|-1)) // (-a)&(+b) = +(~(|a|-1) & b) // (-a)&(-b) = -(((|a|-1)|(|b|-1)) + 1) // // OR: (+a)|(+b) = +(a|b) // (+a)|(-b) = -(((|b|-1) & ~a) + 1) // (-a)|(+b) = -(((|a|-1) & ~b) + 1) // (-a)|(-b) = -(((|a|-1)&(|b|-1)) + 1) // // XOR: (+a)^(+b) = +(a^b) // (+a)^(-b) = -((a^(|b|-1)) + 1) // (-a)^(+b) = -(((|a|-1)^b) + 1) // (-a)^(-b) = +((|a|-1)^(|b|-1)) // // NOT: ~n = -(n+1) // ========================================================================= // Helper: subtract 1 from a word array (for computing |x|-1) // Precondition: the array is non-zero template static void words_sub1(Container& w) { mpn::sub_1(w.data(), w.size(), 1ULL); } // Helper: add 1 to a word array (for computing result + 1) template static void words_add1(Container& w) { uint64_t carry = mpn::add_1(w.data(), w.size(), 1ULL); if (carry) w.push_back(1); } // Helper: normalize — remove leading zero words template static void words_normalize(Container& w) { while (!w.empty() && w.back() == 0) { w.pop_back(); } } void IntOps::bitwiseAnd(Int& result, const Int& other) { if (result.isSpecialState() || other.isSpecialState()) [[unlikely]] { result.setState(NumericState::NaN, NumericError::InvalidBitOperation); return; } bitwiseAndUnchecked(result, other); } void IntOps::bitwiseAndUnchecked(Int& result, const Int& other) { int sa = result.getSign(); int sb = other.getSign(); // AND with zero always absorbs toward zero if (sa == 0) { return; } if (sb == 0) { result = Int::Zero(); return; } if (sa > 0 && sb > 0) { // (+a) & (+b) = +(a & b) size_t minSize = std::min(result.m_words.size(), other.m_words.size()); result.m_words.resize(minSize); for (size_t i = 0; i < minSize; ++i) { result.m_words[i] &= other.m_words[i]; } result.normalize(); } else if (sa > 0 && sb < 0) { // (+a) & (-b) = +(a & ~(|b|-1)) // Two's-complement bits of -b: ~(|b|-1). Extended with 1 in the high bits. SboWords bm = other.m_words; words_sub1(bm); // |b| - 1 size_t aSize = result.m_words.size(); size_t bSize = bm.size(); // Result is at most the size of a (the high -b bits are 1, so a remains as is) for (size_t i = 0; i < aSize; ++i) { uint64_t b_bits = (i < bSize) ? ~bm[i] : UINT64_MAX; result.m_words[i] &= b_bits; } // Result is positive result.normalize(); } else if (sa < 0 && sb > 0) { // (-a) & (+b) = +(~(|a|-1) & b) SboWords am = result.m_words; words_sub1(am); // |a| - 1 size_t aSize = am.size(); size_t bSize = other.m_words.size(); // Result is at most the size of b result.m_words.resize(bSize); for (size_t i = 0; i < bSize; ++i) { uint64_t a_bits = (i < aSize) ? ~am[i] : UINT64_MAX; result.m_words[i] = a_bits & other.m_words[i]; } result.setSign(1); result.normalize(); } else { // (-a) & (-b) = -(((|a|-1) | (|b|-1)) + 1) SboWords am = result.m_words; SboWords bm = other.m_words; words_sub1(am); words_sub1(bm); size_t maxSize = std::max(am.size(), bm.size()); am.resize(maxSize, 0); bm.resize(maxSize, 0); result.m_words.resize(maxSize); for (size_t i = 0; i < maxSize; ++i) { result.m_words[i] = am[i] | bm[i]; } words_add1(result.m_words); result.setSign(-1); result.normalize(); } } void IntOps::bitwiseOr(Int& result, const Int& other) { if (result.isSpecialState() || other.isSpecialState()) [[unlikely]] { result.setState(NumericState::NaN, NumericError::InvalidBitOperation); return; } bitwiseOrUnchecked(result, other); } void IntOps::bitwiseOrUnchecked(Int& result, const Int& other) { int sa = result.getSign(); int sb = other.getSign(); // OR with zero if (sa == 0) { result = other; return; } if (sb == 0) { return; } if (sa > 0 && sb > 0) { // (+a) | (+b) = +(a | b) size_t maxSize = std::max(result.m_words.size(), other.m_words.size()); result.m_words.resize(maxSize, 0); for (size_t i = 0; i < other.m_words.size(); ++i) { result.m_words[i] |= other.m_words[i]; } result.normalize(); } else if (sa > 0 && sb < 0) { // (+a) | (-b) = -(((|b|-1) & ~a) + 1) SboWords bm = other.m_words; words_sub1(bm); size_t aSize = result.m_words.size(); size_t bSize = bm.size(); size_t maxSize = std::max(aSize, bSize); result.m_words.resize(maxSize, 0); bm.resize(maxSize, 0); for (size_t i = 0; i < maxSize; ++i) { uint64_t a_val = result.m_words[i]; result.m_words[i] = bm[i] & ~a_val; } words_normalize(result.m_words); words_add1(result.m_words); result.setSign(-1); result.normalize(); } else if (sa < 0 && sb > 0) { // (-a) | (+b) = -(((|a|-1) & ~b) + 1) SboWords am = result.m_words; words_sub1(am); size_t aSize = am.size(); size_t bSize = other.m_words.size(); size_t maxSize = std::max(aSize, bSize); am.resize(maxSize, 0); result.m_words.resize(maxSize); for (size_t i = 0; i < maxSize; ++i) { uint64_t b_val = (i < bSize) ? other.m_words[i] : 0; result.m_words[i] = am[i] & ~b_val; } words_normalize(result.m_words); words_add1(result.m_words); result.setSign(-1); result.normalize(); } else { // (-a) | (-b) = -(((|a|-1) & (|b|-1)) + 1) SboWords am = result.m_words; SboWords bm = other.m_words; words_sub1(am); words_sub1(bm); size_t minSize = std::min(am.size(), bm.size()); result.m_words.resize(minSize); for (size_t i = 0; i < minSize; ++i) { result.m_words[i] = am[i] & bm[i]; } words_normalize(result.m_words); words_add1(result.m_words); result.setSign(-1); result.normalize(); } } void IntOps::bitwiseXor(Int& result, const Int& other) { if (result.isSpecialState() || other.isSpecialState()) [[unlikely]] { result.setState(NumericState::NaN, NumericError::InvalidBitOperation); return; } bitwiseXorUnchecked(result, other); } void IntOps::bitwiseXorUnchecked(Int& result, const Int& other) { int sa = result.getSign(); int sb = other.getSign(); // XOR with zero if (sa == 0) { result = other; return; } if (sb == 0) { return; } if (sa > 0 && sb > 0) { // (+a) ^ (+b) = +(a ^ b) size_t maxSize = std::max(result.m_words.size(), other.m_words.size()); result.m_words.resize(maxSize, 0); for (size_t i = 0; i < other.m_words.size(); ++i) { result.m_words[i] ^= other.m_words[i]; } result.normalize(); } else if (sa > 0 && sb < 0) { // (+a) ^ (-b) = -((a ^ (|b|-1)) + 1) SboWords bm = other.m_words; words_sub1(bm); size_t aSize = result.m_words.size(); size_t bSize = bm.size(); size_t maxSize = std::max(aSize, bSize); result.m_words.resize(maxSize, 0); bm.resize(maxSize, 0); for (size_t i = 0; i < maxSize; ++i) { result.m_words[i] ^= bm[i]; } words_normalize(result.m_words); words_add1(result.m_words); result.setSign(-1); result.normalize(); } else if (sa < 0 && sb > 0) { // (-a) ^ (+b) = -(((|a|-1) ^ b) + 1) SboWords am = result.m_words; words_sub1(am); size_t aSize = am.size(); size_t bSize = other.m_words.size(); size_t maxSize = std::max(aSize, bSize); am.resize(maxSize, 0); result.m_words.resize(maxSize); for (size_t i = 0; i < maxSize; ++i) { uint64_t b_val = (i < bSize) ? other.m_words[i] : 0; result.m_words[i] = am[i] ^ b_val; } words_normalize(result.m_words); words_add1(result.m_words); result.setSign(-1); result.normalize(); } else { // (-a) ^ (-b) = +((|a|-1) ^ (|b|-1)) SboWords am = result.m_words; SboWords bm = other.m_words; words_sub1(am); words_sub1(bm); size_t maxSize = std::max(am.size(), bm.size()); am.resize(maxSize, 0); bm.resize(maxSize, 0); result.m_words.resize(maxSize); for (size_t i = 0; i < maxSize; ++i) { result.m_words[i] = am[i] ^ bm[i]; } result.setSign(1); result.normalize(); } } void IntOps::bitwiseNot(Int& result) { if (result.isSpecialState()) [[unlikely]] { result.setState(NumericState::NaN, NumericError::InvalidBitOperation); return; } bitwiseNotUnchecked(result); } void IntOps::bitwiseNotUnchecked(Int& result) { // ~n = -(n + 1) (two's-complement identity) if (result.getSign() == 0) { // ~0 = -1 result = Int(-1); } else if (result.getSign() > 0) { // ~(+m) = -(m + 1) words_add1(result.m_words); result.setSign(-1); } else { // ~(-m) = m - 1 words_sub1(result.m_words); result.setSign(1); result.normalize(); // m-1 may become 0 (when m=1) } } void IntOps::leftShift(Int& value, int shift) { // Handle special states if (value.isSpecialState()) { return; } // A shift of 0 or a zero value does nothing if (value.getSign() == 0 || shift == 0) { return; } // Word-granularity shift (in 64-bit units) int wordShift = shift / 64; int bitShift = shift % 64; // Compute the result size size_t newSize = value.m_words.size() + wordShift + (bitShift > 0 ? 1 : 0); std::vector result(newSize, 0); // Bit-granularity shift if (bitShift > 0) { uint64_t carry = 0; for (size_t i = 0; i < value.m_words.size(); ++i) { uint64_t word = value.m_words[i]; result[i + wordShift] = (word << bitShift) | carry; carry = word >> (64 - bitShift); } result[value.m_words.size() + wordShift] = carry; } else { // No bit shift (word-granularity shift only) for (size_t i = 0; i < value.m_words.size(); ++i) { result[i + wordShift] = value.m_words[i]; } } value.m_words.assign(result.data(), result.data() + result.size()); value.normalize(); } void IntOps::rightShift(Int& value, int shift) { // Handle special states if (value.isSpecialState()) { return; } // A shift of 0 or a zero value does nothing if (value.getSign() == 0 || shift == 0) { return; } // Two's-complement semantics: right shift of a negative value is floor division // -7 >> 1 = -4 (floor(-7/2) = -4), not -3 (truncation) // If any bit shifted out of a negative value is non-zero, add 1 to the magnitude before shifting bool needRoundUp = false; if (value.getSign() < 0) { // Check whether any bit being shifted out contains a 1 for (int i = 0; i < std::min(static_cast(value.m_words.size()), shift / 64); ++i) { if (value.m_words[i] != 0) { needRoundUp = true; break; } } if (!needRoundUp && (shift % 64) > 0) { int wordIdx = shift / 64; if (wordIdx < static_cast(value.m_words.size())) { uint64_t mask = (1ULL << (shift % 64)) - 1; if (value.m_words[wordIdx] & mask) { needRoundUp = true; } } } } // Word-granularity shift (in 64-bit units) int wordShift = shift / 64; int bitShift = shift % 64; // When everything is shifted out if (wordShift >= static_cast(value.m_words.size())) { if (needRoundUp) { // Negative value with all bits shifted out → floor = -1 value.m_words.resize(1); value.m_words[0] = 1; value.setSign(-1); } else { value.m_words.clear(); value.setSign(0); } return; } if (bitShift == 0) { // Word-granularity only: in-place removal via erase (no alloc) value.m_words.erase(value.m_words.begin(), value.m_words.begin() + wordShift); } else { // Bit + word shift: handled in two stages // Step 1: word-granularity erase (memmove) if (wordShift > 0) { value.m_words.erase(value.m_words.begin(), value.m_words.begin() + wordShift); } // Step 2: in-place bit shift (wordShift is done, so no overlap) size_t n = value.m_words.size(); uint64_t carry = 0; for (int i = static_cast(n) - 1; i >= 0; --i) { uint64_t w = value.m_words[i]; value.m_words[i] = (w >> bitShift) | carry; carry = w << (64 - bitShift); } // Remove the MSB zero word while (!value.m_words.empty() && value.m_words.back() == 0) { value.m_words.pop_back(); } } // Floor correction for negative values: add 1 to the magnitude if (needRoundUp) { uint64_t c = mpn::add_1(value.m_words.data(), value.m_words.size(), 1ULL); if (c) { value.m_words.push_back(1); } } value.normalize(); } // 3-argument leftShift (checked) void IntOps::leftShift(const Int& value, int shift, Int& result) { if (value.isSpecialState()) [[unlikely]] { if (&result != &value) result = value; return; } leftShiftUnchecked(value, shift, result); } // 3-argument leftShift (unchecked): result = value << shift (buffer reuse) // // ★ Precondition: `&result != &value` (aliasing forbidden). // For the in-place case use `IntOps::leftShift(Int&, int)` (the 1-argument version). // // Reason (latent bug, 2026-05-08 audit): // The memset below zeroes `rp[0..wordShift-1]`, but when `&result == &value` // this memset simultaneously destroys `ap[0..wordShift-1]` (rp == ap). The immediately following // `mpn::lshift(rp + wordShift, ap, n, bitShift)` reads ap[0..n-1], so // the destroyed low wordShift limbs are mixed into the output as zeros. // This function is currently not called in-place in production (bench / IntSequence's // fibonacci2_intops is the only user, and both use separate buffers), so it has not triggered. void IntOps::leftShiftUnchecked(const Int& value, int shift, Int& result) { if (value.getSign() == 0 || shift == 0) { if (&result != &value) result = value; return; } size_t n = value.m_words.size(); int wordShift = shift / 64; int bitShift = shift % 64; size_t newSize = n + wordShift + (bitShift > 0 ? 1 : 0); result.m_words.resize_uninitialized(newSize); uint64_t* rp = result.m_words.data(); const uint64_t* ap = value.m_words.data(); // Zero-fill the low words (when aliased, also destroys the same region of ap — see the precondition above) std::memset(rp, 0, wordShift * sizeof(uint64_t)); if (bitShift == 0) { std::memcpy(rp + wordShift, ap, n * sizeof(uint64_t)); } else { // Shift directly via mpn::lshift (ASM-capable) uint64_t overflow = mpn::lshift(rp + wordShift, ap, n, bitShift); rp[wordShift + n] = overflow; } // Normalize (remove leading zeros) while (newSize > 0 && rp[newSize - 1] == 0) newSize--; result.m_words.resize(newSize); result.m_sign = value.m_sign; result.m_state = NumericState::Normal; } // 3-argument rightShift (checked) void IntOps::rightShift(const Int& value, int shift, Int& result) { if (value.isSpecialState()) [[unlikely]] { if (&result != &value) result = value; return; } rightShiftUnchecked(value, shift, result); } // 3-argument rightShift (unchecked): result = value >> shift (buffer reuse) void IntOps::rightShiftUnchecked(const Int& value, int shift, Int& result) { if (value.getSign() == 0 || shift == 0) { if (&result != &value) result = value; return; } int wordShift = shift / 64; int bitShift = shift % 64; size_t n = value.m_words.size(); // Completely shifted out if (wordShift >= static_cast(n)) { if (value.getSign() < 0) { // Floor division of a negative value: -1 result.m_words.resize(1); result.m_words[0] = 1; result.m_sign = -1; } else { result.m_words.clear(); result.m_sign = 0; } result.m_state = NumericState::Normal; return; } // Floor correction check for negative values bool needRoundUp = false; if (value.getSign() < 0) { for (int i = 0; i < std::min(static_cast(n), wordShift); ++i) { if (value.m_words[i] != 0) { needRoundUp = true; break; } } if (!needRoundUp && bitShift > 0 && wordShift < static_cast(n)) { uint64_t mask = (1ULL << bitShift) - 1; if (value.m_words[wordShift] & mask) needRoundUp = true; } } size_t srcN = n - wordShift; const uint64_t* ap = value.m_words.data() + wordShift; result.m_words.resize_uninitialized(srcN); uint64_t* rp = result.m_words.data(); if (bitShift == 0) { std::memcpy(rp, ap, srcN * sizeof(uint64_t)); } else { mpn::rshift(rp, ap, srcN, bitShift); } // Normalize while (srcN > 0 && rp[srcN - 1] == 0) srcN--; result.m_words.resize(srcN); if (srcN == 0) { result.m_sign = needRoundUp ? -1 : 0; if (needRoundUp) { result.m_words.resize(1); result.m_words[0] = 1; } } else { result.m_sign = value.m_sign; if (needRoundUp) { uint64_t c = mpn::add_1(result.m_words.data(), result.m_words.size(), 1ULL); if (c) result.m_words.push_back(1); } } result.m_state = NumericState::Normal; } // Power-of-two coefficient Int IntOps::pow2(uint64_t exponent) { // Word-granularity shift (in 64-bit units) int wordShift = (int)(exponent / 64); int bitShift = exponent % 64; // Initialize the result Int result; result.m_words.resize(wordShift + 1, 0); result.m_words[wordShift] = 1ULL << bitShift; result.setSign(1); return result; } // GCD (greatest common divisor) Int IntOps::gcd(const Int& a, const Int& b) { // Delegate to Lehmer GCD (IntGCD) return IntGCD::gcd(a, b); } // Implementation of factor removal // Equivalent to GMP's mpz_remove: with the ladder algorithm (repeated factor^(2^k) doubling), // it compresses naive successive division (count times) into log(count) divmods + squarings. // For factor=2 (a power of two), there is a specialized path via trailing zero scan. uint64_t IntOps::removeFactor(const Int& value, const Int& factor, Int& result) { // Special-state checks if (value.isNaN() || factor.isNaN()) { result = Int::NaN(); return 0; } if (value.isInfinite() || factor.isInfinite()) { result = Int::NaN(); return 0; } // factor of 0 or ±1 is invalid if (factor.isZero() || abs(factor).isOne()) { result = Int::NaN(); return 0; } // Case where value is 0 if (value.isZero()) { result = Int::Zero(); return 0; } // GMP-compatible: the factor is always treated by its absolute value Int abs_factor = abs(factor); const int sign = value.getSign(); // Fast path: if the absolute value of factor is a power of two, decide immediately via trailing zero scan // (popcount == 1 ⇔ 2^k form) if (abs_factor.popcount() == 1) { const size_t k = abs_factor.countTrailingZeros(); // factor = 2^k, k >= 1 const size_t tz = (sign < 0 ? (-value).countTrailingZeros() : value.countTrailingZeros()); const uint64_t count = static_cast(tz / k); const size_t shift = count * k; if (shift == 0) { result = value; return 0; } IntOps::rightShift(value, static_cast(shift), result); return count; } // General ladder algorithm (same structure as GMP mpz_remove) // Ladder up: while successively computing fpow[p] = factor^(2^p), // keep going as long as dest is divisible by fpow[p]. // Each time dest /= fpow[p] succeeds, increment p and prepare fpow[p+1] = fpow[p]^2. // Ladder down: descend from the p at which division last failed, // handling the remainder with fpow[p-1], fpow[p-2], ..., fpow[0]. // LADDER_MAX = 63 covers counts up to 2^64 - 1. constexpr int LADDER_MAX = 63; Int fpow[LADDER_MAX + 1]; fpow[0] = std::move(abs_factor); Int dest = (sign < 0 ? -value : value); Int rem; int p = 0; for (;;) { Int q = IntOps::divmod(dest, fpow[p], rem); if (!rem.isZero()) { break; } if (p + 1 >= LADDER_MAX + 1) { // Overflow guard: normally not reached dest = std::move(q); p++; break; } fpow[p + 1] = fpow[p] * fpow[p]; dest = std::move(q); p++; } if (p == 0) { // Not divisible from the start result = value; return 0; } // By this point, removed once each with fpow[0..p-1] → cumulative exponent = 1+2+...+2^(p-1) = 2^p - 1 uint64_t count = (static_cast(1) << p) - 1; // Ladder down: try fpow[p-1], fpow[p-2], ..., fpow[0] while (p > 0) { p--; Int q = IntOps::divmod(dest, fpow[p], rem); if (rem.isZero()) { dest = std::move(q); count += static_cast(1) << p; } } // Restore the sign (GMP-compatible: the factor is treated by its absolute value, but the original sign is kept in the result) result = std::move(dest); if (sign < 0) result = -result; return count; } // Square root (delegated to IntSqrt) Int IntOps::sqrt(const Int& value) { return IntSqrt::sqrt(value); } // Multiplication via the Karatsuba method void IntOps::karatsubaMultiply(const Int& a, const Int& b, Int& result) { // Use the IntMultiplication class's Karatsuba implementation result = IntMultiplication::karatsubaMultiply(a, b); } // Multiplication via Toom-Cook-3 void IntOps::toomCookMultiply(const Int& a, const Int& b, Int& result) { // Use the IntMultiplication class's Toom-Cook-3 implementation result = IntMultiplication::toomCook3Multiply(a, b); } // Multiplication via Toom-Cook-4 void IntOps::toomCook4Multiply(const Int& a, const Int& b, Int& result) { // Use the IntMultiplication class's Toom-Cook-4 implementation result = IntMultiplication::toomCook4Multiply(a, b); } // General-purpose multiplication (handles unbalanced operands) // Writes directly into result.m_words, eliminating the intermediate copy from fromRawWords void IntOps::generalMultiply(const Int& a, const Int& b, Int& result) { size_t an = a.m_words.size(); size_t bn = b.m_words.size(); ScratchScope scope; size_t scratch_size = mpn::multiply_scratch_size(an, bn); uint64_t* scratch = (scratch_size > 0) ? getThreadArena().alloc_limbs(scratch_size) : nullptr; size_t rn = an + bn; result.m_words.resize_uninitialized(rn); mpn::multiply(result.m_words.data(), a.m_words.data(), an, b.m_words.data(), bn, scratch); // normalize: remove leading zero words while (!result.m_words.empty() && result.m_words.back() == 0) result.m_words.pop_back(); result.setSign(1); result.setState(NumericState::Normal); } // Absolute-value multiplication with NTT cache: caches and reuses the forward NTT of rhs void IntOps::mulAbsCached(const Int& lhs, const Int& rhs, Int& result, prime_ntt::NttCache& cache) { size_t an = lhs.m_words.size(); size_t bn = rhs.m_words.size(); // Below the NTT threshold no cache is needed → fall back to normal multiplication if (an == 0 || bn == 0) { result.m_words.clear(); result.setSign(0); result.setState(NumericState::Normal); return; } if (std::min(an, bn) < mpn::PRIME_NTT_THRESHOLD) { multiplyAbsolute(lhs, rhs, result); result.setSign(1); result.setState(NumericState::Normal); return; } size_t rn = an + bn; result.m_words.resize_uninitialized(rn); prime_ntt::mul_prime_ntt_cached(result.m_words.data(), lhs.m_words.data(), an, rhs.m_words.data(), bn, cache); while (!result.m_words.empty() && result.m_words.back() == 0) result.m_words.pop_back(); result.setSign(1); result.setState(NumericState::Normal); } // YC-2: Fused multiply-add: result = a*b + c*d (including sign handling) void IntOps::mulAdd(const Int& a, const Int& b, const Int& c, const Int& d, Int& result) { int sa = a.getSign(), sb = b.getSign(); int sc = c.getSign(), sd = d.getSign(); int sign_ab = sa * sb; // sign of a*b int sign_cd = sc * sd; // sign of c*d // If either product is zero → plain multiplication if (sign_ab == 0 && sign_cd == 0) { result = Int(0); return; } if (sign_ab == 0) { mul(c, d, result); return; } if (sign_cd == 0) { mul(a, b, result); return; } size_t an = a.m_words.size(); size_t bn = b.m_words.size(); size_t cn = c.m_words.size(); size_t dn = d.m_words.size(); // When the signs differ: NTT fusion is not possible → fall back to normal multiplication + addition // (subtraction mod p is possible inside the NTT, but signed large numbers cannot be handled during CRT reconstruction) if (sign_ab != sign_cd) { result = a * b; result += c * d; return; } // From here, sign_ab == sign_cd (both positive or both negative) // Compute |a*b| + |c*d| via NTT fusion and apply the common sign // Fuse only when the smaller side of both is at or above the NTT threshold size_t min_ab = std::min(an, bn); size_t min_cd = std::min(cn, dn); if (min_ab < mpn::PRIME_NTT_THRESHOLD || min_cd < mpn::PRIME_NTT_THRESHOLD) { // Fallback result = a * b; result += c * d; return; } // NTT fusion path: |a|*|b| + |c|*|d| size_t rn = std::max(an + bn, cn + dn) + 1; // +1 for carry from addition result.m_words.resize_uninitialized(rn); prime_ntt::mul_add_prime_ntt(result.m_words.data(), rn, a.m_words.data(), an, b.m_words.data(), bn, c.m_words.data(), cn, d.m_words.data(), dn); while (!result.m_words.empty() && result.m_words.back() == 0) result.m_words.pop_back(); result.setSign(result.m_words.empty() ? 0 : sign_ab); result.setState(NumericState::Normal); } // Multiplication using FFT void IntOps::fftMultiply(const Int& a, const Int& b, Int& result) { // Use the IntMultiplication class's FFT implementation result = IntMultiplication::fftMultiply(a, b); } // Division via Newton's method void IntOps::newtonDivision(const Int& dividend, const Int& divisor, Int& result) { // Use fast division via Newton's method (for very large numbers) Int remainder; result = IntDivision::divNewton(dividend, divisor, remainder); } // Implementation of floor division Int IntOps::floorDiv(const Int& dividend, const Int& divisor) { // Special-state checks if (dividend.isNaN() || divisor.isNaN()) { return Int::NaN(); } if (divisor.isZero()) { return Int::NaN(); // division by zero } if (dividend.isInfinite() || divisor.isInfinite()) { // ∞ / finite = ∞, finite / ∞ = 0, ∞ / ∞ = NaN if (dividend.isInfinite() && !divisor.isInfinite()) { return (dividend.getSign() * divisor.getSign() > 0) ? Int::PositiveInfinity() : Int::NegativeInfinity(); } else if (!dividend.isInfinite() && divisor.isInfinite()) { return Int::Zero(); } else { return Int::NaN(); } } // Normal division Int quotient = dividend / divisor; Int remainder = dividend % divisor; // If the remainder is 0, it is the same as normal division if (remainder.isZero()) { return quotient; } // If the signs differ, decrement the quotient by 1 (round toward negative infinity) if ((dividend.getSign() > 0 && divisor.getSign() < 0) || (dividend.getSign() < 0 && divisor.getSign() > 0)) { quotient = quotient - Int::One(); } return quotient; } // Return floor division and remainder simultaneously Int IntOps::floorDivMod(const Int& dividend, const Int& divisor, Int& remainder) { if (dividend.isNaN() || divisor.isNaN()) { remainder = Int::NaN(); return Int::NaN(); } if (divisor.isZero()) { remainder = Int::NaN(); return Int::NaN(); } Int quotient = divmod(dividend, divisor, remainder); // If the remainder is non-zero and the signs differ, adjust if (!remainder.isZero() && ((dividend.getSign() > 0 && divisor.getSign() < 0) || (dividend.getSign() < 0 && divisor.getSign() > 0))) { quotient = quotient - Int::One(); remainder = remainder + divisor; } return quotient; } // Implementation of floor modulo Int IntOps::floorMod(const Int& dividend, const Int& divisor) { Int remainder; floorDivMod(dividend, divisor, remainder); return remainder; } // Implementation of ceiling division Int IntOps::ceilDiv(const Int& dividend, const Int& divisor) { // Special-state checks if (dividend.isNaN() || divisor.isNaN()) { return Int::NaN(); } if (divisor.isZero()) { return Int::NaN(); // division by zero } if (dividend.isInfinite() || divisor.isInfinite()) { // ∞ / finite = ∞, finite / ∞ = 0, ∞ / ∞ = NaN if (dividend.isInfinite() && !divisor.isInfinite()) { return (dividend.getSign() * divisor.getSign() > 0) ? Int::PositiveInfinity() : Int::NegativeInfinity(); } else if (!dividend.isInfinite() && divisor.isInfinite()) { return Int::Zero(); } else { return Int::NaN(); } } // Normal division Int quotient = dividend / divisor; Int remainder = dividend % divisor; // If the remainder is 0, it is the same as normal division if (remainder.isZero()) { return quotient; } // If the signs are the same, increment the quotient by 1 (round toward positive infinity) if ((dividend.getSign() > 0 && divisor.getSign() > 0) || (dividend.getSign() < 0 && divisor.getSign() < 0)) { quotient = quotient + Int::One(); } return quotient; } // Implementation of exact division (multi-limb divisor: light via Hensel lifting) Int IntOps::divExact(const Int& dividend, const Int& divisor) { // Special-state checks if (dividend.isNaN() || divisor.isNaN()) { return Int::NaN(); } if (divisor.isZero()) { return Int::NaN(); // division by zero } if (dividend.isInfinite() || divisor.isInfinite()) { if (dividend.isInfinite() && !divisor.isInfinite()) { return (dividend.getSign() * divisor.getSign() > 0) ? Int::PositiveInfinity() : Int::NegativeInfinity(); } else if (!dividend.isInfinite() && divisor.isInfinite()) { return Int::Zero(); } else { return Int::NaN(); } } Int result = dividend; IntOps::divExactInPlace(result, divisor); return result; } // In-place exact division: result /= divisor (precondition: divisor evenly divides result) // 1-limb divisor: divExactWord (Hensel) / multi-limb: mpn::divexact void IntOps::divExactInPlace(Int& result, const Int& divisor) { if (result.isSpecialState() || divisor.isSpecialState()) { return; } if (divisor.isZero()) { result = Int::NaN(); return; } if (result.isZero()) { return; } // Result sign = result.sign × divisor.sign int new_sign = result.getSign() * divisor.getSign(); // 1-limb divisor: delegate to the existing divExactWord if (divisor.size() == 1) { divExactWord(result, divisor.word(0)); // divExactWord preserves the sign (word is treated as unsigned) // Flip the sign only if divisor is negative if (divisor.getSign() < 0 && result.m_sign != 0) { result.m_sign = -result.m_sign; } return; } // multi-limb divisor: mpn::divexact (Hensel lifting) size_t an = result.m_words.size(); size_t bn = divisor.m_words.size(); if (an < bn) { // Cannot happen for true exact division, but err on the safe side with 0 result.m_words.clear(); result.m_sign = 0; return; } ScratchScope scope; size_t qn_max = an - bn + 1; uint64_t* q_buf = getThreadArena().alloc_limbs(qn_max + 1); uint64_t* scratch = getThreadArena().alloc_limbs(mpn::divexact_scratch_size(an, bn)); size_t qn = mpn::divexact(q_buf, result.m_words.data(), an, divisor.m_words.data(), bn, scratch); if (qn == 0) { result.m_words.clear(); result.m_sign = 0; } else { result.m_words.assign(q_buf, q_buf + qn); result.m_sign = new_sign; } } } // namespace sangi