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

// IntOps.cpp
// Implementation of utilities for multi-precision integer operations

#include <math/core/mp/Int/IntOps.hpp>
#include <math/core/mp/Int/IntSpecialStates.hpp>
#include <math/core/mp/Int/IntDivision.hpp>
#include <math/core/mp/Int/IntMultiplication.hpp>
#include <math/core/mp/Int/IntSqrt.hpp>
#include <math/core/mp/Int/IntGCD.hpp>
#include <algorithm>
#include <stdexcept>

// With C++20 <bit>, 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<int>(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<int>(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<int>(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<int>(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<const uint64_t>(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<int>(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<int>((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<int>(result.m_words.size()) - 1; i >= 0; --i) {
                std::pair<uint64_t, uint64_t> 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<uint64_t> quotient(absDividend.m_words.size(), 0);
            for (int i = static_cast<int>(absDividend.m_words.size()) - 1; i >= 0; --i) {
                UInt128 current(r, absDividend.m_words[i]);
                std::pair<uint64_t, uint64_t> 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<const uint64_t>(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<const uint64_t>(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<int>((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<typename Container>
    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<typename Container>
    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<typename Container>
    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<uint64_t> 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<int>(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<int>(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<int>(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<int>(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<int>(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<int>(n), wordShift); ++i) {
                if (value.m_words[i] != 0) { needRoundUp = true; break; }
            }
            if (!needRoundUp && bitShift > 0 && wordShift < static_cast<int>(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<uint64_t>(tz / k);
            const size_t shift = count * k;
            if (shift == 0) {
                result = value;
                return 0;
            }
            IntOps::rightShift(value, static_cast<int>(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<uint64_t>(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<uint64_t>(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