#include #include #include #include // std::popcount
#include // std::unique_ptr

namespace sangi {

    //------------------------------------------------------------------------------
    // Arithmetic operators
    //------------------------------------------------------------------------------

    Int operator+(const Int& lhs, const Int& rhs) {
        // Handle special states
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] {
            return IntSpecialStates::handleAddition(lhs, rhs);
        }

        // Optimization for addition with zero
        if (lhs.getSign() == 0) return rhs;
        if (rhs.getSign() == 0) return lhs;

        Int result;
        if (lhs.getSign() == rhs.getSign()) {
            // Same sign: add absolute values (3-argument version avoids a copy)
            IntOps::addAbsolute(lhs, rhs, result);
            result.setSign(lhs.getSign());
        }
        else {
            // Different sign: subtract absolute values (pass lhs, rhs directly — m_words is the absolute value)
            if (IntOps::compareAbsLess(lhs, rhs)) {
                IntOps::subtractAbsolute(rhs, lhs, result);
                result.setSign(rhs.getSign());
            }
            else {
                IntOps::subtractAbsolute(lhs, rhs, result);
                result.setSign(lhs.getSign());
            }
        }

        return result;
    }

    Int operator+(Int&& lhs, const Int& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]]
            return IntSpecialStates::handleAddition(lhs, rhs);
        if (lhs.getSign() == 0) return rhs;
        if (rhs.getSign() == 0) return std::move(lhs);

        if (lhs.getSign() == rhs.getSign()) {
            IntOps::addAbsolute(lhs, rhs);
            return std::move(lhs);
        } else {
            if (IntOps::compareAbsLess(lhs, rhs)) {
                Int result;
                IntOps::subtractAbsolute(rhs, lhs, result);
                result.setSign(rhs.getSign());
                return result;
            } else {
                IntOps::subtractAbsoluteInPlace(lhs, rhs);
                return std::move(lhs);
            }
        }
    }

    Int operator+(const Int& lhs, Int&& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]]
            return IntSpecialStates::handleAddition(lhs, rhs);
        if (rhs.getSign() == 0) return lhs;
        if (lhs.getSign() == 0) return std::move(rhs);

        if (lhs.getSign() == rhs.getSign()) {
            IntOps::addAbsolute(rhs, lhs);
            return std::move(rhs);
        } else {
            if (IntOps::compareAbsLess(rhs, lhs)) {
                Int result;
                IntOps::subtractAbsolute(lhs, rhs, result);
                result.setSign(lhs.getSign());
                return result;
            } else {
                IntOps::subtractAbsoluteInPlace(rhs, lhs);
                return std::move(rhs);
            }
        }
    }

    Int operator+(Int&& lhs, Int&& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]]
            return IntSpecialStates::handleAddition(lhs, rhs);
        if (lhs.getSign() == 0) return std::move(rhs);
        if (rhs.getSign() == 0) return std::move(lhs);

        if (lhs.getSign() == rhs.getSign()) {
            // Reuse the larger buffer (minimize resizes)
            if (lhs.size() >= rhs.size()) {
                IntOps::addAbsolute(lhs, rhs);
                return std::move(lhs);
            } else {
                IntOps::addAbsolute(rhs, lhs);
                return std::move(rhs);
            }
        } else {
            if (IntOps::compareAbsLess(lhs, rhs)) {
                IntOps::subtractAbsoluteInPlace(rhs, lhs);
                return std::move(rhs);
            } else {
                IntOps::subtractAbsoluteInPlace(lhs, rhs);
                return std::move(lhs);
            }
        }
    }

    Int operator-(const Int& lhs, const Int& rhs) {
        // Handle special states
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] {
            return IntSpecialStates::handleSubtraction(lhs, rhs);
        }

        // Optimization for subtraction with zero
        if (lhs.getSign() == 0) {
            Int result = rhs;
            if (result.getSign() != 0) {
                result.setSign(-result.getSign());
            }
            return result;
        }
        if (rhs.getSign() == 0) {
            return lhs;
        }

        // Compute lhs - rhs directly (avoid copying rhs + going through operator+)
        // Effective sign when rhs's sign is negated
        int rhsEffSign = -rhs.getSign();

        Int result;
        if (lhs.getSign() == rhsEffSign) {
            // Same sign: add absolute values
            IntOps::addAbsolute(lhs, rhs, result);
            result.setSign(lhs.getSign());
        }
        else {
            // Different sign: subtract absolute values
            if (IntOps::compareAbsLess(lhs, rhs)) {
                IntOps::subtractAbsolute(rhs, lhs, result);
                result.setSign(rhsEffSign);
            }
            else {
                IntOps::subtractAbsolute(lhs, rhs, result);
                result.setSign(lhs.getSign());
            }
        }

        return result;
    }

    Int operator-(Int&& lhs, const Int& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]]
            return IntSpecialStates::handleSubtraction(lhs, rhs);
        if (lhs.getSign() == 0) {
            Int result = rhs;
            if (result.getSign() != 0) result.setSign(-result.getSign());
            return result;
        }
        if (rhs.getSign() == 0) return std::move(lhs);

        int rhsEffSign = -rhs.getSign();
        if (lhs.getSign() == rhsEffSign) {
            // Same sign (e.g. 5-(-3)=8): add absolute values
            IntOps::addAbsolute(lhs, rhs);
            return std::move(lhs);
        } else {
            // Different sign (e.g. 5-3=2): subtract absolute values
            if (IntOps::compareAbsLess(lhs, rhs)) {
                Int result;
                IntOps::subtractAbsolute(rhs, lhs, result);
                result.setSign(rhsEffSign);
                return result;
            } else {
                IntOps::subtractAbsoluteInPlace(lhs, rhs);
                return std::move(lhs);
            }
        }
    }

    Int operator-(const Int& lhs, Int&& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]]
            return IntSpecialStates::handleSubtraction(lhs, rhs);
        if (lhs.getSign() == 0) {
            // -rhs: negate the sign and reuse the buffer
            if (rhs.getSign() != 0) rhs.setSign(-rhs.getSign());
            return std::move(rhs);
        }
        if (rhs.getSign() == 0) return lhs;

        int rhsEffSign = -rhs.getSign();
        if (lhs.getSign() == rhsEffSign) {
            // Same sign: add absolute values, sign comes from lhs
            IntOps::addAbsolute(rhs, lhs);
            rhs.setSign(lhs.getSign());
            return std::move(rhs);
        } else {
            // Different sign: subtract absolute values
            if (IntOps::compareAbsLess(lhs, rhs)) {
                // |lhs| < |rhs|: reuse the rhs buffer
                IntOps::subtractAbsoluteInPlace(rhs, lhs);
                rhs.setSign(rhsEffSign);
                return std::move(rhs);
            } else {
                Int result;
                IntOps::subtractAbsolute(lhs, rhs, result);
                result.setSign(lhs.getSign());
                return result;
            }
        }
    }

    Int operator-(Int&& lhs, Int&& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]]
            return IntSpecialStates::handleSubtraction(lhs, rhs);
        if (lhs.getSign() == 0) {
            if (rhs.getSign() != 0) rhs.setSign(-rhs.getSign());
            return std::move(rhs);
        }
        if (rhs.getSign() == 0) return std::move(lhs);

        int rhsEffSign = -rhs.getSign();
        if (lhs.getSign() == rhsEffSign) {
            // Same sign: add absolute values
            if (lhs.size() >= rhs.size()) {
                IntOps::addAbsolute(lhs, rhs);
                return std::move(lhs);
            } else {
                IntOps::addAbsolute(rhs, lhs);
                rhs.setSign(lhs.getSign());
                return std::move(rhs);
            }
        } else {
            // Different sign: subtract absolute values
            if (IntOps::compareAbsLess(lhs, rhs)) {
                IntOps::subtractAbsoluteInPlace(rhs, lhs);
                rhs.setSign(rhsEffSign);
                return std::move(rhs);
            } else {
                IntOps::subtractAbsoluteInPlace(lhs, rhs);
                return std::move(lhs);
            }
        }
    }

    Int operator*(const Int& lhs, const Int& rhs) {
        // Handle special states
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] {
            return IntSpecialStates::handleMultiplication(lhs, rhs);
        }

        // Multiplication with zero
        if (lhs.getSign() == 0 || rhs.getSign() == 0) [[unlikely]] {
            return Int::Zero();
        }

        // Determine the sign
        int resultSign = lhs.getSign() * rhs.getSign();

        // Squaring detection: a * a → dedicated squaring path (faster by exploiting symmetry)
        // Special states already handled at entry → call Unchecked directly
        Int result;
        if (&lhs == &rhs) {
            IntOps::squareUnchecked(lhs, result);
        } else {
            IntOps::multiplyAbsolute(lhs, rhs, result);
        }
        result.setSign(resultSign);

        return result;
    }

    // Common edge-case handling for division (special states and division by zero)
    // Return value: true = result is finalized (written to out), false = proceed to the normal path
    static bool handleDivisionEdgeCases(const Int& lhs, const Int& rhs, Int& out) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] {
            out = IntSpecialStates::handleDivision(lhs, rhs);
            return true;
        }
        if (rhs.getSign() == 0) [[unlikely]] {
            if (lhs.getSign() == 0) {
                // 0/0 = NaN
                out = Int::NaN();
                out.setState(NumericState::NaN, NumericError::DivideByZero);
            } else {
                // x/0 (x≠0) = ±∞ (matches IntSpecialStates::handleDivision)
                out = (lhs.getSign() > 0) ? Int::PositiveInfinity() : Int::NegativeInfinity();
            }
            return true;
        }
        if (lhs.getSign() == 0) {
            out = Int::Zero();
            return true;
        }
        return false;
    }

    Int operator/(const Int& lhs, const Int& rhs) {
        Int early;
        if (handleDivisionEdgeCases(lhs, rhs, early)) return early;

        int resultSign = lhs.getSign() * rhs.getSign();
        Int result = lhs;
        result.setSign(1);
        // Edge cases already handled at entry → call Unchecked directly
        IntOps::divideAbsoluteUnchecked(result, rhs);
        if (result.getSign() != 0) result.setSign(resultSign);
        return result;
    }

    Int operator/(Int&& lhs, const Int& rhs) {
        Int early;
        if (handleDivisionEdgeCases(lhs, rhs, early)) return early;

        int resultSign = lhs.getSign() * rhs.getSign();
        lhs.setSign(1);
        IntOps::divideAbsoluteUnchecked(lhs, rhs);
        if (lhs.getSign() != 0) lhs.setSign(resultSign);
        return std::move(lhs);
    }

    Int operator%(const Int& lhs, const Int& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] {
            return IntSpecialStates::handleModulo(lhs, rhs);
        }
        if (rhs.getSign() == 0) [[unlikely]] {
            Int result = Int::NaN();
            result.setState(NumericState::NaN, NumericError::DivideByZero);
            return result;
        }
        if (lhs.getSign() == 0) return Int::Zero();

        Int result = lhs;
        result.setSign(1);
        IntOps::moduloAbsoluteUnchecked(result, rhs);
        if (lhs.getSign() < 0 && result.getSign() != 0) result.setSign(-1);
        return result;
    }

    Int operator%(Int&& lhs, const Int& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] {
            return IntSpecialStates::handleModulo(lhs, rhs);
        }
        if (rhs.getSign() == 0) [[unlikely]] {
            Int result = Int::NaN();
            result.setState(NumericState::NaN, NumericError::DivideByZero);
            return result;
        }
        if (lhs.getSign() == 0) return Int::Zero();

        int origSign = lhs.getSign();
        lhs.setSign(1);
        IntOps::moduloAbsoluteUnchecked(lhs, rhs);
        if (origSign < 0 && lhs.getSign() != 0) lhs.setSign(-1);
        return std::move(lhs);
    }

    //------------------------------------------------------------------------------
    // Increment/decrement operators
    //------------------------------------------------------------------------------

    Int& operator++(Int& value) {
        if (value.isSpecialState()) return value;
        IntOps::addDelta(value, 1);
        return value;
    }

    Int operator++(Int& value, int) {
        if (value.isSpecialState()) return value;
        Int old = value;
        IntOps::addDelta(value, 1);
        return old;
    }

    Int& operator--(Int& value) {
        if (value.isSpecialState()) return value;
        IntOps::addDelta(value, -1);
        return value;
    }

    Int operator--(Int& value, int) {
        if (value.isSpecialState()) return value;
        Int old = value;
        IntOps::addDelta(value, -1);
        return old;
    }

    //------------------------------------------------------------------------------
    // Bitwise operators
    //------------------------------------------------------------------------------

    // Helper that produces a NaN (InvalidBitOperation)
    static Int makeBitOpNaN() {
        Int result = Int::NaN();
        result.setState(NumericState::NaN, NumericError::InvalidBitOperation);
        return result;
    }

    Int operator&(const Int& lhs, const Int& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] return makeBitOpNaN();
        Int result = lhs;
        IntOps::bitwiseAndUnchecked(result, rhs);
        return result;
    }

    Int operator&(Int&& lhs, const Int& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] return makeBitOpNaN();
        IntOps::bitwiseAndUnchecked(lhs, rhs);
        return std::move(lhs);
    }

    Int operator&(const Int& lhs, Int&& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] return makeBitOpNaN();
        IntOps::bitwiseAndUnchecked(rhs, lhs);  // AND is commutative
        return std::move(rhs);
    }

    Int operator&(Int&& lhs, Int&& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] return makeBitOpNaN();
        IntOps::bitwiseAndUnchecked(lhs, rhs);
        return std::move(lhs);
    }

    Int operator|(const Int& lhs, const Int& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] return makeBitOpNaN();
        Int result = lhs;
        IntOps::bitwiseOrUnchecked(result, rhs);
        return result;
    }

    Int operator|(Int&& lhs, const Int& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] return makeBitOpNaN();
        IntOps::bitwiseOrUnchecked(lhs, rhs);
        return std::move(lhs);
    }

    Int operator|(const Int& lhs, Int&& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] return makeBitOpNaN();
        IntOps::bitwiseOrUnchecked(rhs, lhs);  // OR is commutative
        return std::move(rhs);
    }

    Int operator|(Int&& lhs, Int&& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] return makeBitOpNaN();
        IntOps::bitwiseOrUnchecked(lhs, rhs);
        return std::move(lhs);
    }

    // Exponentiation operator: a ^ b = pow(a, b)
    Int operator^(const Int& lhs, const Int& rhs) {
        return pow(lhs, rhs);
    }

    // Bitwise XOR (function form)
    Int bitwiseXor(const Int& lhs, const Int& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] return makeBitOpNaN();
        Int result = lhs;
        IntOps::bitwiseXorUnchecked(result, rhs);
        return result;
    }

    Int bitwiseXor(Int&& lhs, const Int& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] return makeBitOpNaN();
        IntOps::bitwiseXorUnchecked(lhs, rhs);
        return std::move(lhs);
    }

    Int bitwiseXor(const Int& lhs, Int&& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] return makeBitOpNaN();
        IntOps::bitwiseXorUnchecked(rhs, lhs);  // XOR is commutative
        return std::move(rhs);
    }

    Int bitwiseXor(Int&& lhs, Int&& rhs) {
        if (lhs.isSpecialState() || rhs.isSpecialState()) [[unlikely]] return makeBitOpNaN();
        IntOps::bitwiseXorUnchecked(lhs, rhs);
        return std::move(lhs);
    }

    Int operator~(const Int& value) {
        if (value.isSpecialState()) return makeBitOpNaN();
        // Mathematical definition of bit inversion: ~x = -(x+1)
        Int result = value + 1;
        result = -result;
        return result;
    }

    Int operator~(Int&& value) {
        if (value.isSpecialState()) return makeBitOpNaN();
        // Mathematical definition of bit inversion: ~x = -(x+1)
        value += 1;
        value = -value;
        return std::move(value);
    }

    Int operator<<(const Int& value, int shift) {
        if (value.isSpecialState()) return value;
        if (shift < 0) return value >> (-shift);
        Int result = value;
        IntOps::leftShift(result, shift);
        return result;
    }

    Int operator<<(Int&& value, int shift) {
        if (value.isSpecialState()) return std::move(value);
        if (shift < 0) return std::move(value) >> (-shift);
        IntOps::leftShift(value, shift);
        return std::move(value);
    }

    Int operator>>(const Int& value, int shift) {
        if (value.isSpecialState()) return value;
        if (shift < 0) return value << (-shift);
        Int result = value;
        IntOps::rightShift(result, shift);
        return result;
    }

    Int operator>>(Int&& value, int shift) {
        if (value.isSpecialState()) return std::move(value);
        if (shift < 0) return std::move(value) << (-shift);
        IntOps::rightShift(value, shift);
        return std::move(value);
    }

    //------------------------------------------------------------------------------
    // Comparison operators
    //------------------------------------------------------------------------------

    std::partial_ordering operator<=>(const Int& lhs, const Int& rhs) {
        // Fast path: both Normal (no special-state check needed)
        // Equivalent to GMP's mpz_cmp: signed size comparison → word comparison
        if (lhs.m_state == NumericState::Normal &&
            rhs.m_state == NumericState::Normal) [[likely]] {
            // Sign comparison
            if (lhs.m_sign != rhs.m_sign)
                return lhs.m_sign <=> rhs.m_sign;
            if (lhs.m_sign == 0)
                return std::partial_ordering::equivalent;

            // Size comparison (absolute value)
            size_t ln = lhs.m_words.size();
            size_t rn = rhs.m_words.size();
            if (ln != rn) {
                bool less = (ln < rn);
                if (lhs.m_sign < 0) less = !less;
                return less ? std::partial_ordering::less : std::partial_ordering::greater;
            }

            // Word comparison (from the MSB)
            const uint64_t* lp = lhs.m_words.data();
            const uint64_t* rp = rhs.m_words.data();
            for (int i = static_cast (ln) - 1; i >= 0; --i) {
                if (lp[i] != rp[i]) {
                    bool less = (lp[i] < rp[i]);
                    if (lhs.m_sign < 0) less = !less;
                    return less ? std::partial_ordering::less : std::partial_ordering::greater;
                }
            }
            return std::partial_ordering::equivalent;
        }

        // Slow path: special states (NaN, Infinity)
        if (lhs.isNaN() || rhs.isNaN())
            return std::partial_ordering::unordered;

        if (lhs.getState() == NumericState::NegativeInfinity) {
            return (rhs.getState() == NumericState::NegativeInfinity)
                ? std::partial_ordering::equivalent : std::partial_ordering::less;
        }
        if (lhs.getState() == NumericState::PositiveInfinity) {
            return (rhs.getState() == NumericState::PositiveInfinity)
                ? std::partial_ordering::equivalent : std::partial_ordering::greater;
        }
        if (rhs.getState() == NumericState::NegativeInfinity)
            return std::partial_ordering::greater;
        if (rhs.getState() == NumericState::PositiveInfinity)
            return std::partial_ordering::less;

        // We reach here when one side is Normal and the other is a special state
        // (NaN/Inf are already handled by the conditions above)
        return std::partial_ordering::equivalent;
    }

    bool operator==(const Int& lhs, const Int& rhs) {
        if (lhs.isNaN() || rhs.isNaN()) {
            return false;
        }
        return (lhs <=> rhs) == std::partial_ordering::equivalent;
    }

    //------------------------------------------------------------------------------
    // Utility functions
    //------------------------------------------------------------------------------

    Int abs(const Int& value) {
        // Fast path: Normal state (exclude special states in a single branch).
        // The old code chained 4 branches isNaN / isInfinite / isZero / getSign() < 0
        // + an Int copy, taking 18ns for a 1-limb value, 4.4x slower than GMP mpz_abs (4ns).
        if (value.getState() == NumericState::Normal) [[likely]] {
            Int result = value;
            // Write the sign directly (preserve sign=0 for empty mantissa: if the original
            // sign is 0 it stays 0 after the copy, so setSign's zero guard is unnecessary)
            if (result.getSign() < 0) {
                // Writing m_sign directly is private, so we go through setSign via friend, but
                // the combination of sign != 0 and an empty mantissa does not depend on the source state:
                // original sign < 0 → original mantissa non-empty → still non-empty after the copy → safe
                result.setSign(1);
            }
            return result;
        }
        // Slow path: special states
        if (value.isNaN()) return value;
        if (value.isInfinite()) {
            return (value.getState() == NumericState::NegativeInfinity)
                ? Int::PositiveInfinity() : value;
        }
        // What remains is PositiveZero / NegativeZero / Overflow / Underflow / ComplexInfinity
        // |0| is +0; abs(±Overflow/Underflow/ComplexInfinity) is the value itself
        if (value.getState() == NumericState::NegativeZero) return Int::Zero();
        return value;
    }

    Int abs(Int&& value) {
        // Fast path: Normal state
        if (value.getState() == NumericState::Normal) [[likely]] {
            if (value.getSign() < 0) value.setSign(1);
            return std::move(value);
        }
        // Slow path
        if (value.isNaN()) return std::move(value);
        if (value.isInfinite()) {
            if (value.getState() == NumericState::NegativeInfinity)
                return Int::PositiveInfinity();
            return std::move(value);
        }
        if (value.getState() == NumericState::NegativeZero) return Int::Zero();
        return std::move(value);
    }

    Int gcd(const Int& a, const Int& b) {
        // Handle special states
        if (a.isSpecialState() || b.isSpecialState()) {
            return Int::NaN();
        }

        // Handle zero
        if (a.getSign() == 0) {
            return abs(b);
        }
        if (b.getSign() == 0) {
            return abs(a);
        }

        // Use the optimized internal implementation
        return IntOps::gcd(a, b);
    }

    Int lcm(const Int& a, const Int& b) {
        // Handle special states
        if (a.isSpecialState() || b.isSpecialState()) {
            return Int::NaN();
        }

        // The least common multiple with zero is 0
        if (a.getSign() == 0 || b.getSign() == 0) {
            return Int(0);
        }

        Int gcd_value = gcd(a, b);
        // a and b are already special-state checked → use the Unchecked version to avoid a double check
        Int quot;
        IntOps::divUnchecked(a, gcd_value, quot);
        Int product;
        IntOps::mulUnchecked(quot, b, product);
        if (product.getSign() < 0) product.setSign(1);
        return product;
    }

    Int pow(const Int& base, unsigned int exponent) {
        // Handle special cases
        if (base.isSpecialState()) {
            if (exponent == 0) {
                return Int(1); // Normally special-value^0 = 1
            }
            if (base.isNaN()) {
                return Int::NaN();
            }
            if (base.isInfinite()) {
                if (exponent % 2 == 0) {
                    return Int::PositiveInfinity(); // ±∞^(2n) = +∞
                }
                else {
                    return base.getSign() > 0 ?
                        Int::PositiveInfinity() : Int::NegativeInfinity(); // ±∞^(2n+1) = ±∞
                }
            }
        }

        // Basic cases
        //[Note]0^0=1 holds
        //[Reference]In C++ pow(0, 0) is implementation-defined, but the standard library often returns 1

        if (exponent == 0) {
            return Int(1); // b^0 = 1 (any value to the 0th power is 1)
        }
        if (base.isZero()) {
            return Int::Zero(); // 0^n = 0 (when n is a positive integer)
        }
        if (exponent == 1) {
            return base; // b^1 = b
        }

        // Special handling for -1
        if (base == -1) {
            return (exponent % 2 == 0) ? Int(1) : Int(-1); // (-1)^even = 1, (-1)^odd = -1
        }

        // Special handling for 1
        if (base.isOne()) {
            return Int(1); // 1^n = 1
        }

        // Overflow protection for a huge exponent.
        // This is actually within the limits of the unsigned int type, so there is no case
        // extreme enough to warrant consideration in Int pow(const Int& base, const Int& exponent),
        // but it is included for optimization, taking computational cost into account.
        if (exponent > 1000000000) { // Overflow-prevention threshold
            Int absBase = abs(base);
            if (absBase > 1) {
                // When |base| > 1, the result becomes very large
                if (base.isNegative() && (exponent % 2 != 0)) {
                    return Int::NegativeInfinity(); // When a negative base is raised to an odd power
                }
                return Int::PositiveInfinity();
            }
            // When |base| == 1, normal processing (1 or -1)
            // When absBase == 1, normal processing (1 or -1)
        }

        // Left-to-right binary exponentiation
        // Same number of multiplications as the right-to-left method, but avoids an unnecessary trailing squaring.
        // Example: exp=6 (110₂) → x², x³=x²·x, x⁶=(x³)² with 3 multiplications
        int bits = 0;
        {
            unsigned int tmp = exponent;
            while (tmp > 1) { bits++; tmp >>= 1; }
        }

        Int result = base;
        Int temp;  // Reuse the temporary buffer
        for (int i = bits - 1; i >= 0; i--) {
            IntOps::mulUnchecked(result, result, temp);
            std::swap(result, temp);
            if ((exponent >> i) & 1) {
                IntOps::mulUnchecked(result, base, temp);
                std::swap(result, temp);
            }
        }

        return result;
    }

    Int pow(const Int& base, const Int& exponent) {
        // Handle special cases
        if (base.isSpecialState()) {
            if (exponent.isZero()) {
                return Int(1); // Normally special-value^0 = 1
            }
            if (base.isNaN() || exponent.isNaN()) {
                return Int::NaN();
            }
            if (base.isInfinite()) {
                if (exponent.isNegative()) {
                    return Int::Zero(); // ∞^(-n) = 0
                }
                // Check whether the exponent is even
                Int two(2);
                Int modResult = exponent % two;
                if (modResult.isZero()) {
                    return Int::PositiveInfinity(); // ±∞^(2n) = +∞
                }
                else {
                    return base.getSign() > 0 ?
                        Int::PositiveInfinity() : Int::NegativeInfinity(); // ±∞^(2n+1) = ±∞
                }
            }
        }

        // Handle the exponent being a special state
        if (exponent.isSpecialState()) {
            if (exponent.isNaN()) {
                return Int::NaN();
            }
            if (exponent.isInfinite()) {
                if (exponent.isNegative()) {
                    // Special handling when base is 1 or -1
                    if (base.isOne()) {
                        return Int(1); // 1^(-∞) = 1
                    }
                    else if (base == -1) {
                        // -1^(-∞) is regarded as undefined, so return NaN
                        return Int::NaN();
                    }
                    return Int::Zero(); // 0 if |base| > 1 or |base| < 1
                }
                else {
                    // When base is 1
                    if (base.isOne()) {
                        return Int(1); // 1^∞ = 1
                    }
                    // When base is -1
                    if (base == -1) {
                        // -1^∞ is regarded as undefined, so return NaN
                        return Int::NaN();
                    }
                    // ∞ if |base| > 1, 0 if |base| < 1
                    Int absBase = abs(base);
                    Int one(1);
                    if (absBase > one) {
                        return Int::PositiveInfinity();
                    }
                    else if (absBase < one) {
                        return Int::Zero();
                    }
                }
            }
        }

        // Basic cases
        //[Note]0^0=1 holds
        //[Reference]In C++ pow(0, 0) is implementation-defined, but the standard library often returns 1

        if (exponent.isZero()) {
            return Int(1);
        }
        if (base.isZero()) {
            if (exponent.isNegative()) {
                // 0^(-n) is undefined (division by zero)
                return Int::NaN();
            }
            return Int::Zero();
        }
        if (exponent.isOne()) {
            return base;
        }
        if (exponent == -1) {
            // b^(-1) is 1/b, but since Int is an integer it gets special handling
            if (base.isOne() || base == -1) {
                return base; // 1^(-1) = 1, (-1)^(-1) = -1
            }
            return Int::Zero(); // The reciprocal of an integer is fractional, so as an Int it is 0
        }

        // Handle a negative exponent
        if (exponent.isNegative()) {
            // Within the integer range, the result is 0 if |base| > 1
            // and ±1 if |base| = 1
            Int absBase = abs(base);
            if (absBase > 1) {
                return Int::Zero();
            }
            else if (absBase.isOne()) {
                // -1 if base is -1 and exponent is odd
                if (base.isNegative()) {
                    Int absExp = abs(exponent);
                    if (absExp.isOdd()) {
                        return Int(-1);
                    }
                }
                return Int(1);
            }
            return Int::Zero();
        }

        // When exponent is very large: treat as overflow
        const size_t MAX_UINT = std::numeric_limits ::max();
        if (exponent.bitLength() > 32 || exponent.toUInt64() > MAX_UINT) {
            // For a very large exponent, treat it as overflow
            if (abs(base) > 1) {
                return base.isNegative() && exponent.getBit(0) ?
                    Int::NegativeInfinity() : Int::PositiveInfinity();
            }
            else if (base.isOne()) {
                return Int(1);
            }
            else if (base == -1) {
                // Determine whether -1 is raised to an even or odd power
                return exponent.getBit(0) ? Int(-1) : Int(1);
            }
            else { // |base| < 1
                return Int::Zero();
            }
        }

        // Normal exponent computation (exponent is within the range of unsigned int)
        unsigned int expValue = static_cast (exponent.toUInt64());

        // Computation by repeated squaring (3-argument version reuses the buffer)
        Int result(1);
        Int power = base;
        Int temp;
        // Special states already handled at entry → use Unchecked to reduce redundant checks
        while (expValue > 0) {
            if (expValue & 1) {
                IntOps::mulUnchecked(result, power, temp);
                std::swap(result, temp);
            }
            IntOps::squareUnchecked(power, temp);
            std::swap(power, temp);
            expValue >>= 1;
        }
        return result;
    }

    Int powMod(const Int& base, const Int& exponent, const Int& modulus) {
        // Delegate to Montgomery modular exponentiation (handles both odd and even m)
        return IntModular::powerMod(base, exponent, modulus);
    }

    std::size_t bitCount(const Int& value) {
        // Count the number of set bits
        if (value.isNaN() || value.isInfinite() || value.getSign() < 0) {
            throw std::invalid_argument("bitCount requires non-negative normal value");
        }

        // When the value is 0
        if (value.getSign() == 0) {
            return 0;
        }

        std::size_t count = 0;
        for (std::size_t i = 0; i < value.size(); ++i) {
            count += static_cast (std::popcount(value.word(i)));
        }

        return count;
    }

    std::pair removeFactor(const Int& value, const Int& factor) {
        Int result;
        uint64_t count = IntOps::removeFactor(value, factor, result);
        return { std::move(result), count };
    }

    bool isProbablePrime(const Int& value, int iterations) {
        // Miller-Rabin primality test

        // Handle special cases
        if (value.isNaN() || value.isInfinite() || value.getSign() < 0) {
            return false;
        }

        // Handle small values
        if (value.isZero() || value.isOne()) {
            return false;
        }

        if (value == 2 || value == 3) {
            return true;
        }

        // Even check
        if (value.isEven()) {
            return false; // Even numbers other than 2 are not prime
        }

        // Basic implementation of the Miller-Rabin method
        // Decompose into the form value-1 = 2^s * d
        Int d = value - 1;
        int s = 0;
        while (d.isEven()) {
            d >>= 1;
            s++;
        }

        // It should have randomness, but this simple implementation uses fixed values
        Int a_values[] = { Int(2), Int(3), Int(5), Int(7), Int(11), Int(13), Int(17), Int(19), Int(23) };
        int num_tests = (iterations < 9) ? iterations : 9;

        for (int i = 0; i < num_tests; ++i) {
            Int a = a_values[i];

            // Skip if it is greater than or equal to value
            if (a >= value) {
                continue;
            }

            // Compute x = a^d mod value
            Int x = powMod(a, d, value);

            if (x.isOne() || x == value - 1) {
                continue;
            }

            bool is_witness = true;
            for (int j = 1; j < s; ++j) {
                // x and value are normal at entry → squareUnchecked is safe
                { Int sq; IntOps::squareUnchecked(x, sq); x = sq % value; }
                if (x == value - 1) {
                    is_witness = false;
                    break;
                }
            }

            if (is_witness) {
                return false; // Composite
            }
        }

        return true; // Probably prime
    }

    //--------------------------------------------------------------------------
    // Single-word arithmetic operators (equivalent to GMP _ui/_si)
    //--------------------------------------------------------------------------

    // === addScalar ===

    Int addScalar(const Int& lhs, uint64_t rhs) {
        if (lhs.isSpecialState()) [[unlikely]] return lhs;
        if (rhs == 0) return lhs;
        if (lhs.getSign() == 0) return Int(rhs);
        Int result(lhs);
        IntOps::addWord(result, rhs);
        return result;
    }

    Int addScalar(Int&& lhs, uint64_t rhs) {
        if (lhs.isSpecialState()) [[unlikely]] return std::move(lhs);
        if (rhs == 0) return std::move(lhs);
        if (lhs.getSign() == 0) return Int(rhs);
        IntOps::addWord(lhs, rhs);
        return std::move(lhs);
    }

    // === subScalar ===

    Int subScalar(const Int& lhs, uint64_t rhs) {
        if (lhs.isSpecialState()) [[unlikely]] return lhs;
        if (rhs == 0) return lhs;
        if (lhs.getSign() == 0) { Int r(rhs); r.setSign(-1); return r; }
        Int result(lhs);
        IntOps::subWord(result, rhs);
        return result;
    }

    Int subScalar(Int&& lhs, uint64_t rhs) {
        if (lhs.isSpecialState()) [[unlikely]] return std::move(lhs);
        if (rhs == 0) return std::move(lhs);
        if (lhs.getSign() == 0) {
            lhs = Int(rhs);
            lhs.setSign(-1);
            return std::move(lhs);
        }
        IntOps::subWord(lhs, rhs);
        return std::move(lhs);
    }

    Int rsubScalar(uint64_t lhs, const Int& rhs) {
        if (rhs.isSpecialState()) [[unlikely]] return -rhs;
        if (rhs.getSign() == 0) return Int(lhs);
        if (lhs == 0) return -rhs;
        Int result(rhs);
        result.setSign(-result.getSign());
        IntOps::addWord(result, lhs);
        return result;
    }

    // === mulScalar ===

    Int mulScalar(const Int& lhs, uint64_t rhs) {
        if (lhs.isSpecialState()) [[unlikely]] return lhs;
        if (lhs.getSign() == 0 || rhs == 0) return Int(0);
        if (rhs == 1) return lhs;
        Int result(lhs);
        IntOps::multiplyWordUnchecked(result, rhs);
        return result;
    }

    Int mulScalar(Int&& lhs, uint64_t rhs) {
        if (lhs.isSpecialState()) [[unlikely]] return std::move(lhs);
        if (lhs.getSign() == 0 || rhs == 0) return Int(0);
        if (rhs == 1) return std::move(lhs);
        IntOps::multiplyWordUnchecked(lhs, rhs);
        return std::move(lhs);
    }

    // === divScalar ===

    Int divScalar(const Int& lhs, uint64_t rhs) {
        if (rhs == 0) { Int r; r.setState(NumericState::NaN); return r; }
        if (lhs.isSpecialState()) [[unlikely]] return lhs;
        if (lhs.getSign() == 0) return Int(0);
        if (rhs == 1) return lhs;
        Int result(lhs);
        IntOps::divWord(result, rhs);
        return result;
    }

    Int divScalar(Int&& lhs, uint64_t rhs) {
        if (rhs == 0) { Int r; r.setState(NumericState::NaN); return r; }
        if (lhs.isSpecialState()) [[unlikely]] return std::move(lhs);
        if (lhs.getSign() == 0) return std::move(lhs);
        if (rhs == 1) return std::move(lhs);
        IntOps::divWord(lhs, rhs);
        return std::move(lhs);
    }

    // === modScalarU (remainder uint64_t) ===

    uint64_t modScalarU(const Int& lhs, uint64_t rhs) {
        if (rhs == 0 || lhs.isSpecialState() || lhs.getSign() == 0)
            return 0;
        size_t n = lhs.size();
        constexpr size_t STACK_LIMIT = 64;
        uint64_t stack_buf[STACK_LIMIT];
        std::unique_ptr heap_buf;
        uint64_t* q = stack_buf;
        if (n > STACK_LIMIT) {
            heap_buf.reset(new uint64_t[n]);
            q = heap_buf.get();
        }
        return mpn::divmod_1(q, lhs.data(), n, rhs);
    }

    // === int64_t version (decompose the sign and delegate to the uint64_t version) ===

    static uint64_t abs64(int64_t v) {
        return (v >= 0) ? static_cast (v)
                        : static_cast (-(v + 1)) + 1u;
    }

    Int addScalar(const Int& lhs, int64_t rhs) {
        return (rhs >= 0) ? addScalar(lhs, static_cast (rhs))
                          : subScalar(lhs, abs64(rhs));
    }
    Int addScalar(Int&& lhs, int64_t rhs) {
        return (rhs >= 0) ? addScalar(std::move(lhs), static_cast (rhs))
                          : subScalar(std::move(lhs), abs64(rhs));
    }

    Int subScalar(const Int& lhs, int64_t rhs) {
        return (rhs >= 0) ? subScalar(lhs, static_cast (rhs))
                          : addScalar(lhs, abs64(rhs));
    }
    Int subScalar(Int&& lhs, int64_t rhs) {
        return (rhs >= 0) ? subScalar(std::move(lhs), static_cast (rhs))
                          : addScalar(std::move(lhs), abs64(rhs));
    }
    Int rsubScalar(int64_t lhs, const Int& rhs) {
        if (lhs >= 0) return rsubScalar(static_cast (lhs), rhs);
        return -(addScalar(rhs, abs64(lhs)));
    }

    Int mulScalar(const Int& lhs, int64_t rhs) {
        Int result = mulScalar(lhs, abs64(rhs));
        if (rhs < 0 && result.getSign() != 0) result.setSign(-result.getSign());
        return result;
    }
    Int mulScalar(Int&& lhs, int64_t rhs) {
        Int result = mulScalar(std::move(lhs), abs64(rhs));
        if (rhs < 0 && result.getSign() != 0) result.setSign(-result.getSign());
        return result;
    }

    Int divScalar(const Int& lhs, int64_t rhs) {
        Int result = divScalar(lhs, abs64(rhs));
        if (rhs < 0 && result.getSign() != 0) result.setSign(-result.getSign());
        return result;
    }
    Int divScalar(Int&& lhs, int64_t rhs) {
        Int result = divScalar(std::move(lhs), abs64(rhs));
        if (rhs < 0 && result.getSign() != 0) result.setSign(-result.getSign());
        return result;
    }

    Int modScalar(const Int& lhs, int64_t rhs) {
        uint64_t rem = modScalarU(lhs, abs64(rhs));
        if (rem == 0) return Int(0);
        Int result(rem);
        if (lhs.getSign() < 0) result.setSign(-1);
        return result;
    }

    // === Compound assignment (uint64_t) ===

    Int& Int::plusEqScalar(uint64_t rhs) {
        if (isSpecialState()) return *this;
        if (rhs == 0) return *this;
        if (getSign() == 0) { *this = Int(rhs); return *this; }
        IntOps::addWord(*this, rhs);
        return *this;
    }
    Int& Int::minusEqScalar(uint64_t rhs) {
        if (isSpecialState()) return *this;
        if (rhs == 0) return *this;
        if (getSign() == 0) { *this = Int(rhs); setSign(-1); return *this; }
        IntOps::subWord(*this, rhs);
        return *this;
    }
    Int& Int::mulEqScalar(uint64_t rhs) {
        if (isSpecialState()) return *this;
        if (getSign() == 0 || rhs == 0) { *this = Int(0); return *this; }
        if (rhs == 1) return *this;
        IntOps::multiplyWordUnchecked(*this, rhs);
        return *this;
    }
    Int& Int::divEqScalar(uint64_t rhs) {
        if (rhs == 0) { setState(NumericState::NaN); return *this; }
        if (isSpecialState() || getSign() == 0 || rhs == 1) return *this;
        IntOps::divWord(*this, rhs);
        return *this;
    }
    Int& Int::modEqScalar(uint64_t rhs) {
        if (rhs == 0 || isSpecialState() || getSign() == 0) return *this;
        uint64_t rem = *this % rhs;
        int old_sign = getSign();
        *this = Int(rem);
        if (rem != 0 && old_sign < 0) setSign(-1);
        return *this;
    }

    // === Compound assignment (int64_t) ===

    Int& Int::plusEqScalar(int64_t rhs) {
        return (rhs >= 0) ? plusEqScalar(static_cast (rhs))
                          : minusEqScalar(abs64(rhs));
    }
    Int& Int::minusEqScalar(int64_t rhs) {
        return (rhs >= 0) ? minusEqScalar(static_cast (rhs))
                          : plusEqScalar(abs64(rhs));
    }
    Int& Int::mulEqScalar(int64_t rhs) {
        mulEqScalar(abs64(rhs));
        if (rhs < 0 && getSign() != 0) setSign(-getSign());
        return *this;
    }
    Int& Int::divEqScalar(int64_t rhs) {
        divEqScalar(abs64(rhs));
        if (rhs < 0 && getSign() != 0) setSign(-getSign());
        return *this;
    }
    Int& Int::modEqScalar(int64_t rhs) {
        int old_sign = getSign();
        modEqScalar(abs64(rhs));
        if (old_sign < 0 && getSign() > 0) setSign(-1);
        return *this;
    }

    //--------------------------------------------------------------------------
    // Single-word comparison operators
    //--------------------------------------------------------------------------

    std::partial_ordering cmpScalar(const Int& lhs, int64_t rhs) {
        if (lhs.isSpecialState()) [[unlikely]] {
            if (lhs.isNaN()) return std::partial_ordering::unordered;
            if (lhs.isInfinite()) return (lhs.getSign() > 0) ? std::partial_ordering::greater
                                                              : std::partial_ordering::less;
        }
        int rhs_sign = (rhs > 0) ? 1 : (rhs < 0) ? -1 : 0;
        if (lhs.getSign() != rhs_sign) {
            return (lhs.getSign() < rhs_sign) ? std::partial_ordering::less
                                               : std::partial_ordering::greater;
        }
        if (rhs_sign == 0) return std::partial_ordering::equivalent;

        uint64_t abs_rhs = abs64(rhs);
        if (lhs.size() > 1) {
            return (lhs.getSign() > 0) ? std::partial_ordering::greater
                                        : std::partial_ordering::less;
        }
        uint64_t abs_lhs = lhs.word(0);
        if (abs_lhs == abs_rhs) return std::partial_ordering::equivalent;
        bool lhs_abs_greater = (abs_lhs > abs_rhs);
        if (lhs.getSign() < 0) lhs_abs_greater = !lhs_abs_greater;
        return lhs_abs_greater ? std::partial_ordering::greater
                               : std::partial_ordering::less;
    }

    bool eqScalar(const Int& lhs, int64_t rhs) {
        if (lhs.isSpecialState()) return false;
        int rhs_sign = (rhs > 0) ? 1 : (rhs < 0) ? -1 : 0;
        if (lhs.getSign() != rhs_sign) return false;
        if (rhs_sign == 0) return true;
        if (lhs.size() != 1) return false;
        return lhs.word(0) == abs64(rhs);
    }

    std::partial_ordering cmpScalar(const Int& lhs, uint64_t rhs) {
        if (lhs.isSpecialState()) [[unlikely]] {
            if (lhs.isNaN()) return std::partial_ordering::unordered;
            if (lhs.isInfinite()) return (lhs.getSign() > 0) ? std::partial_ordering::greater
                                                              : std::partial_ordering::less;
        }
        if (lhs.getSign() < 0) return std::partial_ordering::less;
        if (lhs.getSign() == 0) {
            return (rhs == 0) ? std::partial_ordering::equivalent
                              : std::partial_ordering::less;
        }
        if (rhs == 0) return std::partial_ordering::greater;
        if (lhs.size() > 1) return std::partial_ordering::greater;
        uint64_t v = lhs.word(0);
        if (v == rhs) return std::partial_ordering::equivalent;
        return (v > rhs) ? std::partial_ordering::greater
                         : std::partial_ordering::less;
    }

    bool eqScalar(const Int& lhs, uint64_t rhs) {
        if (lhs.isSpecialState()) return false;
        if (rhs == 0) return lhs.getSign() == 0;
        if (lhs.getSign() != 1 || lhs.size() != 1) return false;
        return lhs.word(0) == rhs;
    }

} // namespace sangi