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

// IntBase.cpp
// Basic implementation of the Int class

#include <math/core/mp/Int/IntBase.hpp>
#include <math/core/mp/Int/IntSpecialStates.hpp>
#include <math/core/mp/Int/IntOps.hpp>
#include <math/core/mp/Int/IntIO.hpp>
#include <math/core/numeric_state.hpp>
#include <math/core/common.hpp> // for exception classes such as OverflowError
#include <sstream>
#include <string>
#include <cmath>
#include <iomanip> // std::setfill, std::setw, etc.
#include <algorithm> // std::reverse, std::max
#include <bit> // C++20 std::countl_zero, etc.

namespace sangi {

    // Default constructor
    Int::Int() : m_sign(0), m_state(NumericState::Normal), m_error(NumericError::None) {
        // Initialize to zero - m_words may stay empty
    }

    // Constructor from an int value
    Int::Int(int value) : m_state(NumericState::Normal), m_error(NumericError::None) {
        if (value == 0) {
            m_sign = 0;
            // m_words stays empty
        }
        else if (value > 0) {
            m_sign = 1;
            m_words.push_back(static_cast<uint64_t>(value));
        }
        else {
            m_sign = -1;
            m_words.push_back(static_cast<uint64_t>(-static_cast<int64_t>(value)));
        }
    }

    // Constructor from an int64_t value
    Int::Int(int64_t value) : m_state(NumericState::Normal), m_error(NumericError::None) {
        if (value == 0) {
            m_sign = 0;
            // m_words stays empty
        }
        else if (value > 0) {
            m_sign = 1;
            m_words.push_back(static_cast<uint64_t>(value));
        }
        else if (value == INT64_MIN) {
            // INT64_MIN is a special case
            m_sign = -1;
            m_words.push_back(static_cast<uint64_t>(INT64_MAX) + 1);
        }
        else {
            m_sign = -1;
            m_words.push_back(static_cast<uint64_t>(-value));
        }
    }

    // Constructor from a uint64_t value
    Int::Int(uint64_t value) : m_state(NumericState::Normal), m_error(NumericError::None) {
        if (value == 0) {
            m_sign = 0;
            // m_words stays empty
        }
        else {
            m_sign = 1;
            m_words.push_back(value);
        }
    }

    // Constructor from a string
    Int::Int(std::string_view str, int base) {
        *this = IntIOUtils::fromString(str, base);
    }

    // Copy constructor
    Int::Int(const Int& other)
        : m_words(other.m_words),
        m_sign(other.m_sign),
        m_state(other.m_state),
        m_error(other.m_error) {
    }

    // Move constructor
    Int::Int(Int&& other) noexcept
        : m_words(std::move(other.m_words)),
        m_sign(other.m_sign),
        m_state(other.m_state),
        m_error(other.m_error) {
        other.m_sign = 0;
        other.m_state = NumericState::Normal;
        other.m_error = NumericError::None;
    }

    // Assignment operator
    Int& Int::operator=(const Int& other) {
        if (this != &other) {
            m_words = other.m_words;
            m_sign = other.m_sign;
            m_state = other.m_state;
            m_error = other.m_error;
        }
        return *this;
    }

    // Move assignment operator
    Int& Int::operator=(Int&& other) noexcept {
        if (this != &other) {
            m_words = std::move(other.m_words);
            m_sign = other.m_sign;
            m_state = other.m_state;
            m_error = other.m_error;

            other.m_sign = 0;
            other.m_state = NumericState::Normal;
            other.m_error = NumericError::None;
        }
        return *this;
    }

    // Assignment from an int value
    Int& Int::operator=(int value) {
        if (value == 0) {
            m_words.clear();
            m_sign = 0;
        }
        else if (value > 0) {
            m_words.clear();
            m_words.push_back(static_cast<uint64_t>(value));
            m_sign = 1;
        }
        else {
            m_words.clear();
            m_words.push_back(static_cast<uint64_t>(-static_cast<int64_t>(value)));
            m_sign = -1;
        }
        m_state = NumericState::Normal;
        m_error = NumericError::None;
        return *this;
    }

    // Accessor methods for special states
    Int Int::nan() {
        Int result;
        result.m_state = NumericState::NaN;
        result.m_error = NumericError::ExplicitNaN;
        return result;
    }

    Int Int::infinity() {
        Int result;
        result.m_state = NumericState::PositiveInfinity;
        result.m_sign = 1;
        return result;
    }

    Int Int::negInfinity() {
        Int result;
        result.m_state = NumericState::NegativeInfinity;
        result.m_sign = -1;
        return result;
    }

    // Normalize the internal state
    void Int::normalize() {
        // Remove leading zeros from the word array
        while (!m_words.empty() && m_words.back() == 0) {
            m_words.pop_back();
        }

        // Check whether the value is zero
        if (m_words.empty()) {
            m_sign = 0;
        }

        // No artificial upper limit on the word count. On out-of-memory, std::bad_alloc is thrown.
    }

    size_t Int::trimZeroWords() {
        // Remove zero words on the MSB side (equivalent to Int::normalize())
        while (!m_words.empty() && m_words.back() == 0) {
            m_words.pop_back();
        }
        if (m_words.empty()) {
            m_sign = 0;
            return 0;
        }
        // Remove zero words on the LSB side (SboWords::erase = memmove, no alloc)
        size_t lsb = 0;
        while (lsb < m_words.size() - 1 && m_words[lsb] == 0) {
            ++lsb;
        }
        if (lsb > 0) {
            m_words.erase(m_words.begin(), m_words.begin() + lsb);
        }
        return lsb;
    }

    // String conversion
    std::string Int::toString(int base) const {
        return IntIOUtils::toString(*this, base);
    }

    // Hexadecimal string representation
    std::string Int::toHexString(bool uppercase) const {
        // Handle special states
        std::string specialResult = IntSpecialStates::handleToString(*this);
        if (!specialResult.empty()) {
            return specialResult;
        }

        if (isZero()) {
            return "0";
        }

        std::string result;
        // Generate the hexadecimal representation
        for (auto it = m_words.rbegin(); it != m_words.rend(); ++it) {
            std::stringstream ss;
            ss << std::hex;
            if (uppercase) {
                ss << std::uppercase;
            }
            ss << std::setfill('0');
            // Zero-pad to 16 digits for all but the first value
            if (it != m_words.rbegin()) {
                ss << std::setw(16);
            }
            ss << *it;
            result += ss.str();
        }

        // Remove leading zeros
        result.erase(0, result.find_first_not_of('0'));
        if (result.empty()) {
            return "0";
        }

        // Add the sign
        if (m_sign < 0) {
            result = "-" + result;
        }

        return result;
    }

    // Binary string representation
    std::string Int::toBinaryString() const {
        // Handle special states
        std::string specialResult = IntSpecialStates::handleToString(*this);
        if (!specialResult.empty()) {
            return specialResult;
        }

        if (isZero()) {
            return "0";
        }

        std::string result;
        for (auto it = m_words.rbegin(); it != m_words.rend(); ++it) {
            uint64_t word = *it;
            for (int i = 63; i >= 0; --i) {
                if (word & (1ULL << i)) {
                    result += '1';
                }
                else if (!result.empty() || i == 0) {
                    // Ignore leading zeros, but append if a value has already been added or this is the least-significant bit
                    result += '0';
                }
            }
        }

        // Add the sign
        if (m_sign < 0) {
            result = "-" + result;
        }

        return result;
    }

    // Octal string representation
    std::string Int::toOctalString() const {
        // Handle special states
        std::string specialResult = IntSpecialStates::handleToString(*this);
        if (!specialResult.empty()) {
            return specialResult;
        }

        if (isZero()) {
            return "0";
        }

        Int copy = *this;
        if (copy.m_sign < 0) {
            copy.m_sign = 1; // Handle the sign at the end
        }

        std::string result;
        Int base(8);

        // Convert the value to an octal string
        while (!copy.isZero()) {
            Int remainder;
            copy = IntOps::divmod(copy, base, remainder);
            result.push_back('0' + remainder.toInt());
        }

        // Reverse the result (from least-significant to most-significant digit)
        std::reverse(result.begin(), result.end());

        // Add the sign
        if (m_sign < 0) {
            result = "-" + result;
        }

        return result;
    }

    // Conversion to int
    int Int::toInt() const {
        if (isNaN() || isInfinite()) {
            throw OverflowError("Cannot convert NaN or infinity to int");
        }

        if (isZero()) {
            return 0;
        }

        // Can only convert up to one word
        if (m_words.size() > 1) {
            // INT_MIN special case: a negative value that needs a second word
            if (m_words.size() == 2 && m_words[1] == 0 && m_sign < 0) {
                // This case is permitted
            }
            else {
                throw OverflowError("Int value too large for int conversion");
            }
        }

        uint64_t wordValue = m_words[0];

        // For a positive value
        if (m_sign > 0) {
            if (wordValue > static_cast<uint64_t>(std::numeric_limits<int>::max())) {
                throw OverflowError("Int value too large for int conversion");
            }
            return static_cast<int>(wordValue);
        }

        // For a negative value
        // Special case: INT_MIN (-2147483648) is the negation of 2147483648
        if (wordValue == static_cast<uint64_t>(std::numeric_limits<int>::max()) + 1) {
            return std::numeric_limits<int>::min();
        }

        if (wordValue > static_cast<uint64_t>(std::numeric_limits<int>::max())) {
            throw OverflowError("Int value too small for int conversion");
        }

        return -static_cast<int>(wordValue);
    }

    // Conversion to a 64-bit integer
    int64_t Int::toInt64() const {
        if (isNaN() || isInfinite()) {
            throw OverflowError("Cannot convert NaN or infinity to int64_t");
        }

        if (isZero()) {
            return 0;
        }

        if (m_words.size() > 1 ||
            (m_sign > 0 && m_words[0] > static_cast<uint64_t>(std::numeric_limits<int64_t>::max())) ||
            (m_sign < 0 && m_words[0] > static_cast<uint64_t>(std::numeric_limits<int64_t>::max()) + 1)) {
            throw OverflowError("Int value too large for int64_t conversion");
        }

        int64_t value = static_cast<int64_t>(m_words[0]);
        if (m_sign < 0) {
            if (value < 0 && value != INT64_MIN) {
                throw OverflowError("Int value too small for int64_t conversion");
            }
            if (value == INT64_MIN) {
                return INT64_MIN;
            }
            value = -value;
        }
        return value;
    }

    // Conversion to an unsigned 64-bit integer
    uint64_t Int::toUInt64() const {
        if (isNaN() || isInfinite() || m_sign < 0) {
            throw OverflowError("Cannot convert NaN, infinity, or negative value to uint64_t");
        }

        if (isZero()) {
            return 0;
        }

        if (m_words.size() > 1) {
            throw OverflowError("Int value too large for uint64_t conversion");
        }

        return m_words[0];
    }

    // Conversion to size_t
    size_t Int::toSizeT() const {
        if (isNaN() || isInfinite() || m_sign < 0) {
            throw OverflowError("Cannot convert NaN, infinity, or negative value to size_t");
        }

        if (isZero()) {
            return 0;
        }

        if (m_words.size() > 1 || (sizeof(size_t) < sizeof(uint64_t) &&
            m_words[0] > static_cast<uint64_t>(std::numeric_limits<size_t>::max()))) {
            throw OverflowError("Int value too large for size_t conversion");
        }

        return static_cast<size_t>(m_words[0]);
    }

    // Conversion to a double-precision floating-point number
    double Int::toDouble() const {
        if (isNaN()) {
            return std::numeric_limits<double>::quiet_NaN();
        }

        if (m_state == NumericState::PositiveInfinity) {
            return std::numeric_limits<double>::infinity();
        }

        if (m_state == NumericState::NegativeInfinity) {
            return -std::numeric_limits<double>::infinity();
        }

        if (isZero()) {
            return 0.0;
        }

        double result = 0.0;
        double factor = 1.0;

        for (size_t i = 0; i < m_words.size(); ++i) {
            result += static_cast<double>(m_words[i]) * factor;
            factor *= static_cast<double>(UINT64_MAX) + 1.0;
        }

        if (m_sign < 0) {
            result = -result;
        }

        return result;
    }

    double Int::toDouble2Exp(int64_t* exp) const {
        if (isZero() || isNaN() || isInfinite()) {
            *exp = 0;
            if (isNaN()) return std::numeric_limits<double>::quiet_NaN();
            if (m_state == NumericState::PositiveInfinity) return std::numeric_limits<double>::infinity();
            if (m_state == NumericState::NegativeInfinity) return -std::numeric_limits<double>::infinity();
            return 0.0;
        }

        size_t bl = bitLength();
        *exp = static_cast<int64_t>(bl);

        // Take the top 53 bits and build a double in [0.5, 1.0)
        size_t n = m_words.size();
        uint64_t top = m_words[n - 1];
        int top_bits = 64 - std::countl_zero(top);

        // Collect the top 53 bits to build mantissa ∈ [2^52, 2^53)
        uint64_t mantissa;
        if (top_bits >= 53) {
            mantissa = top >> (top_bits - 53);
        } else {
            int need = 53 - top_bits;
            mantissa = top << need;  // First left-shift to widen to 53 bits
            if (n >= 2 && need <= 64) {
                uint64_t w1 = m_words[n - 2];
                mantissa |= (w1 >> (64 - need));
            }
            // n==1 and top_bits < 53: lower bits are filled with 0 (correct)
        }

        // d = mantissa / 2^53 ∈ [0.5, 1.0)
        double d = static_cast<double>(mantissa) * (1.0 / static_cast<double>(1ULL << 53));
        // d ∈ [0.25, 0.5) → fix: double it and decrement exp by 1
        // In fact: mantissa is a 53-bit value with MSB=1, so mantissa ∈ [2^52, 2^53)
        // d = mantissa / 2^53 ∈ [0.5, 1.0) ✓

        if (m_sign < 0) d = -d;
        return d;
    }

    bool Int::isDivisibleBy2Exp(size_t b) const {
        if (b == 0 || isZero()) return true;
        if (isNaN() || isInfinite()) return false;

        size_t full_words = b / 64;
        size_t rem_bits = b % 64;
        size_t n = m_words.size();

        // Whether the lower full_words words are all 0
        size_t check = std::min(full_words, n);
        for (size_t i = 0; i < check; ++i) {
            if (m_words[i] != 0) return false;
        }
        if (full_words >= n) return true;  // If b > bitLength, all bits are within the zero range

        // Mask check on the remaining bits
        if (rem_bits > 0) {
            uint64_t mask = (1ULL << rem_bits) - 1;
            if ((m_words[full_words] & mask) != 0) return false;
        }
        return true;
    }

    bool Int::isCongruent2Exp(const Int& c, size_t b) const {
        if (b == 0) return true;
        if (isNaN() || isInfinite() || c.isNaN() || c.isInfinite()) return false;

        // Whether the lower b bits match (check the lower bits of the difference, accounting for sign)
        // Same sign: directly compare the lower b bits of the magnitudes
        // Different sign: whether the lower b bits of (this - c) are 0
        // Simple implementation: take the difference and use isDivisibleBy2Exp
        // However this incurs an O(n) subtraction. Optimize as needed.
        Int diff = *this - c;
        return diff.isDivisibleBy2Exp(b);
    }

    // Compute the bit length
    size_t Int::bitLength() const {
        if (isZero() || isNaN() || isInfinite()) {
            return 0;
        }

        if (m_words.empty()) {
            return 0;
        }

        // Compute the number of bits in the most-significant word
        uint64_t msw = m_words.back();

        // Use the C++20 function
#if defined(__cpp_lib_bitops) && __cpp_lib_bitops >= 201907L
        int msw_bits = 64 - std::countl_zero(msw);
#else
    // Fallback for older compilers
        int msw_bits = 0;
        for (int i = 63; i >= 0; --i) {
            if ((msw & (1ULL << i)) != 0) {
                msw_bits = i + 1;
                break;
            }
        }
#endif

        // Total number of bits
        return (m_words.size() - 1) * 64 + msw_bits;
    }

    // Get the bit at a specific position
    bool Int::getBit(size_t position) const {
        if (isZero() || isNaN() || isInfinite()) {
            return false;
        }

        size_t word_index = position / 64;
        size_t bit_index = position % 64;

        if (word_index >= m_words.size()) {
            return false;
        }

        return (m_words[word_index] & (1ULL << bit_index)) != 0;
    }

    // Set a bit
    void Int::setBit(size_t position) {
        if (isSpecialState()) return;
        if (m_sign < 0) {
            // Negative number: OR operation with two's-complement semantics
            *this |= (Int(1) << static_cast<int>(position));
            return;
        }
        size_t word_index = position / 64;
        size_t bit_index = position % 64;
        if (word_index >= m_words.size()) {
            m_words.resize(word_index + 1, 0);
        }
        m_words[word_index] |= (1ULL << bit_index);
        if (m_sign == 0) m_sign = 1;
    }

    // Clear a bit
    void Int::clearBit(size_t position) {
        if (isSpecialState()) return;
        if (m_sign < 0) {
            // Negative number: AND operation with two's-complement semantics
            *this &= ~(Int(1) << static_cast<int>(position));
            return;
        }
        size_t word_index = position / 64;
        size_t bit_index = position % 64;
        if (word_index >= m_words.size()) return;  // Already 0
        m_words[word_index] &= ~(1ULL << bit_index);
        normalize();
    }

    // Flip a bit
    void Int::complementBit(size_t position) {
        if (isSpecialState()) return;
        if (m_sign < 0) {
            // Negative number: XOR operation with two's-complement semantics
            // Special states already handled at entry + mask is always Normal from the Int(1)<<n construction → Unchecked
            Int mask = Int(1) << static_cast<int>(position);
            IntOps::bitwiseXorUnchecked(*this, mask);
            return;
        }
        size_t word_index = position / 64;
        size_t bit_index = position % 64;
        if (word_index >= m_words.size()) {
            m_words.resize(word_index + 1, 0);
        }
        m_words[word_index] ^= (1ULL << bit_index);
        if (m_sign == 0) m_sign = 1;
        normalize();
    }

    // scanBit1: position of the first set bit at or after startPos
    size_t Int::scanBit1(size_t startPos) const {
        if (isSpecialState()) return SIZE_MAX;

        if (m_sign < 0) {
            // Negative number: two's complement = ~(|this|-1)
            // scanBit1 = find the first 0-bit in (|this|-1)
            Int mag_minus_1 = abs(*this) - 1;
            size_t word_index = startPos / 64;
            size_t bit_index = startPos % 64;
            // Find a 0 bit within the words of (|this|-1)
            if (word_index < mag_minus_1.m_words.size()) {
                uint64_t word = ~mag_minus_1.m_words[word_index];
                word >>= bit_index;
                if (word != 0) {
                    return startPos + static_cast<size_t>(std::countr_zero(word));
                }
                for (size_t i = word_index + 1; i < mag_minus_1.m_words.size(); i++) {
                    uint64_t w = ~mag_minus_1.m_words[i];
                    if (w != 0) {
                        return i * 64 + static_cast<size_t>(std::countr_zero(w));
                    }
                }
            }
            // Outside storage everything is 0 → complemented, everything is 1
            size_t beyond = mag_minus_1.m_words.size() * 64;
            return (startPos >= beyond) ? startPos : beyond;
        }

        if (isZero()) return SIZE_MAX;

        size_t word_index = startPos / 64;
        size_t bit_index = startPos % 64;
        if (word_index >= m_words.size()) return SIZE_MAX;

        // First word (partial)
        uint64_t word = m_words[word_index] >> bit_index;
        if (word != 0) {
            return startPos + static_cast<size_t>(std::countr_zero(word));
        }
        // Subsequent words
        for (size_t i = word_index + 1; i < m_words.size(); i++) {
            if (m_words[i] != 0) {
                return i * 64 + static_cast<size_t>(std::countr_zero(m_words[i]));
            }
        }
        return SIZE_MAX;
    }

    // scanBit0: position of the first clear bit at or after startPos
    size_t Int::scanBit0(size_t startPos) const {
        if (isSpecialState()) return SIZE_MAX;

        if (m_sign < 0) {
            // Negative number: two's complement = ~(|this|-1)
            // scanBit0 = find the first 1-bit in (|this|-1)
            Int mag_minus_1 = abs(*this) - 1;
            size_t word_index = startPos / 64;
            size_t bit_index = startPos % 64;
            if (word_index < mag_minus_1.m_words.size()) {
                uint64_t word = mag_minus_1.m_words[word_index] >> bit_index;
                if (word != 0) {
                    return startPos + static_cast<size_t>(std::countr_zero(word));
                }
                for (size_t i = word_index + 1; i < mag_minus_1.m_words.size(); i++) {
                    if (mag_minus_1.m_words[i] != 0) {
                        return i * 64 + static_cast<size_t>(std::countr_zero(mag_minus_1.m_words[i]));
                    }
                }
            }
            // Outside storage everything is 0 → complemented, everything is 1 → no 0 bit
            return SIZE_MAX;
        }

        if (isZero()) return startPos;  // All bits 0

        size_t word_index = startPos / 64;
        size_t bit_index = startPos % 64;

        if (word_index < m_words.size()) {
            // Find a 0 bit in the first word (partial)
            uint64_t word = ~m_words[word_index];
            word >>= bit_index;
            if (word != 0) {
                return startPos + static_cast<size_t>(std::countr_zero(word));
            }
            // Subsequent words
            for (size_t i = word_index + 1; i < m_words.size(); i++) {
                if (m_words[i] != UINT64_MAX) {
                    return i * 64 + static_cast<size_t>(std::countr_zero(~m_words[i]));
                }
            }
        }
        // Outside storage everything is 0
        size_t beyond = m_words.size() * 64;
        return (startPos >= beyond) ? startPos : beyond;
    }

    // hammingDistance: number of differing bits between two integers (compatible with GMP mpz_hamdist)
    size_t Int::hammingDistance(const Int& a, const Int& b) {
        // Special states
        if (a.isSpecialState() || b.isSpecialState()) return SIZE_MAX;
        // If the signs differ, return infinity (same as GMP)
        if ((a.m_sign < 0) != (b.m_sign < 0)) return SIZE_MAX;

        if (a.m_sign >= 0 && b.m_sign >= 0) {
            // Both non-negative: popcount of the XOR
            size_t max_size = std::max(a.m_words.size(), b.m_words.size());
            size_t count = 0;
            for (size_t i = 0; i < max_size; i++) {
                uint64_t wa = (i < a.m_words.size()) ? a.m_words[i] : 0;
                uint64_t wb = (i < b.m_words.size()) ? b.m_words[i] : 0;
                count += static_cast<size_t>(std::popcount(wa ^ wb));
            }
            return count;
        }

        // Both negative: two's complement = ~(|x|-1)
        // hamming(-a, -b) = popcount( ~(|a|-1) XOR ~(|b|-1) ) = popcount( (|a|-1) XOR (|b|-1) )
        Int a_mag_m1 = abs(a) - 1;
        Int b_mag_m1 = abs(b) - 1;
        size_t max_size = std::max(a_mag_m1.m_words.size(), b_mag_m1.m_words.size());
        size_t count = 0;
        for (size_t i = 0; i < max_size; i++) {
            uint64_t wa = (i < a_mag_m1.m_words.size()) ? a_mag_m1.m_words[i] : 0;
            uint64_t wb = (i < b_mag_m1.m_words.size()) ? b_mag_m1.m_words[i] : 0;
            count += static_cast<size_t>(std::popcount(wa ^ wb));
        }
        return count;
    }

    // Stream output
    std::ostream& operator<<(std::ostream& os, const Int& value) {
        os << value.toString();
        return os;
    }

    // Stream input
    std::istream& operator>>(std::istream& is, Int& value) {
        std::string str;
        is >> str;
        value = Int(str);
        return is;
    }

    // Bit-counting functions
    size_t Int::countLeadingZeros() const {
        if (isZero() || isNaN() || isInfinite()) {
            return 64;
        }
        if (m_words.empty()) {
            return 64;
        }
        // Number of zero bits in the most-significant word
        uint64_t msw = m_words.back();

        // Use the C++20 function
#if defined(__cpp_lib_bitops) && __cpp_lib_bitops >= 201907L
        size_t zeros = std::countl_zero(msw);
#else
    // Fallback for older compilers
        size_t zeros = 0;
        for (int i = 63; i >= 0; --i) {
            if ((msw & (1ULL << i)) == 0) {
                zeros++;
            }
            else {
                break;
            }
        }
#endif

        // Unused bits in the most-significant word + all bits of the lower words
        return zeros + static_cast<size_t>(64 * (m_words.size() - 1));
    }

    size_t Int::countTrailingZeros() const {
        if (isZero() || isNaN() || isInfinite()) {
            return 64;  // Special cases unchanged
        }

        if (m_words.empty()) {
            return 64;  // Empty case also unchanged
        }

        // Examine words starting from the least-significant
        for (size_t i = 0; i < m_words.size(); ++i) {
            if (m_words[i] != 0) {
                uint64_t word = m_words[i];

                // Use the C++20 function
#if defined(__cpp_lib_bitops) && __cpp_lib_bitops >= 201907L
                return std::countr_zero(word) + i * 64;
#else
            // Fallback for older compilers - fixed version
                size_t zeros = 0;
                for (size_t j = 0; j < 64; ++j) {
                    if ((word & (1ULL << j)) == 0) {
                        zeros++;
                    }
                    else {
                        // Stop once the first 1 bit is found
                        return zeros + i * 64;
                    }
                }
                // If this word is all zeros (normally unreachable)
                return 64 + i * 64;
#endif
            }
        }

        // If all words are zero
        return 64 * m_words.size();
    }

    size_t Int::popcount() const {
        if (isSpecialState()) return SIZE_MAX;
        if (m_sign < 0) return SIZE_MAX;  // Negative number: infinitely many 1s in two's complement
        if (m_sign == 0) return 0;

        size_t count = 0;
        for (size_t i = 0; i < m_words.size(); ++i) {
            count += static_cast<size_t>(std::popcount(m_words[i]));
        }
        return count;
    }

    // Internal absolute-value comparison function
    int Int::compareAbsoluteValues(const Int& lhs, const Int& rhs) {
        // Compare by word count
        if (lhs.m_words.size() < rhs.m_words.size()) {
            return -1;
        }
        if (lhs.m_words.size() > rhs.m_words.size()) {
            return 1;
        }

        // For the same word count, compare from the most-significant word
        for (size_t i = lhs.m_words.size(); i-- > 0; ) {
            if (lhs.m_words[i] < rhs.m_words[i]) {
                return -1;
            }
            if (lhs.m_words[i] > rhs.m_words[i]) {
                return 1;
            }
        }

        // If all are equal
        return 0;
    }

    // Unary minus operator
    // Special-state path of operator-() (NaN/Infinity)
    // Slow path called from the inline version
    Int Int::negate_slow() const {
        if (isNaN()) return *this;
        if (m_state == NumericState::PositiveInfinity) return NegativeInfinity();
        if (m_state == NumericState::NegativeInfinity) return PositiveInfinity();
        // Should not reach here
        return *this;
    }

    // Left-shift assignment operator
    Int& Int::operator<<=(int shift) {
        // Handle special states
        if (isNaN() || isInfinite() || isZero()) {
            return *this;
        }

        // Handle a negative shift as a right shift
        if (shift < 0) {
            return (*this >>= -shift);
        }

        // A zero shift does nothing
        if (shift == 0) {
            return *this;
        }

        // Shift across word boundaries
        size_t word_shift = shift / 64;
        int bit_shift = shift % 64;

        // Expand the size in preparation
        size_t old_size = m_words.size();
        if (bit_shift != 0 && (m_words.empty() || m_words.back() >> (64 - bit_shift) != 0)) {
            // When the bit shift requires an additional word
            m_words.resize(old_size + word_shift + 1);
        }
        else {
            // Otherwise
            m_words.resize(old_size + word_shift);
        }

        // Word-boundary shift
        if (word_shift > 0) {
            // Process from the most-significant word to avoid overwriting data
            for (size_t i = old_size; i-- > 0; ) {
                m_words[i + word_shift] = m_words[i];
            }
            // Zero-clear the lower words
            for (size_t i = 0; i < word_shift; ++i) {
                m_words[i] = 0;
            }
        }

        // Bitwise shift
        if (bit_shift > 0) {
            // Process from the last word
            uint64_t carry = 0;
            for (size_t i = word_shift; i < m_words.size(); ++i) {
                uint64_t new_carry = m_words[i] >> (64 - bit_shift);
                m_words[i] = (m_words[i] << bit_shift) | carry;
                carry = new_carry;
            }
        }

        // Normalize (remove unnecessary leading zeros)
        normalize();
        return *this;
    }

    // Right-shift assignment operator
    Int& Int::operator>>=(int shift) {
        // Handle special states
        if (isNaN() || isInfinite() || isZero()) {
            return *this;
        }

        // Handle a negative shift as a left shift
        if (shift < 0) {
            return (*this <<= -shift);
        }

        // A zero shift does nothing
        if (shift == 0) {
            return *this;
        }

        // Shift across word boundaries
        size_t word_shift = shift / 64;
        if (word_shift >= m_words.size()) {
            // All bits are shifted out
            m_words.clear();
            m_sign = 0;
            return *this;
        }

        // Bit shift within words
        int bit_shift = shift % 64;
        if (bit_shift == 0) {
            // Word-only shift
            m_words.erase(m_words.begin(), m_words.begin() + word_shift);
        }
        else {
            // Word + bitwise shift
            for (size_t i = 0; i < m_words.size() - word_shift; ++i) {
                m_words[i] = (i + word_shift < m_words.size()) ?
                    ((m_words[i + word_shift] >> bit_shift) |
                        ((i + word_shift + 1 < m_words.size()) ? (m_words[i + word_shift + 1] << (64 - bit_shift)) : 0)) : 0;
            }
            m_words.resize(m_words.size() - word_shift);
        }

        // Check whether it is zero
        normalize();
        return *this;
    }

    Int& Int::operator+=(const Int& other) {
        if (isSpecialState() || other.isSpecialState()) {
            *this = IntSpecialStates::handleAddition(*this, other);
            return *this;
        }
        if (other.getSign() == 0) return *this;
        if (getSign() == 0) { *this = other; return *this; }

        if (m_sign == other.m_sign) {
            IntOps::addAbsolute(*this, other);
        } else {
            if (IntOps::compareAbsLess(*this, other)) {
                // ★ |*this| < |other|: expand the buffer, copy other, then subtract in place
                size_t on = other.m_words.size();
                size_t tn = m_words.size();
                // Temporarily save the data of *this
                constexpr size_t BUF_LIMIT = 64;
                uint64_t stack_buf[BUF_LIMIT];
                uint64_t* tbuf = (tn <= BUF_LIMIT) ? stack_buf : new uint64_t[tn];
                std::copy_n(m_words.data(), tn, tbuf);
                // Copy other into *this
                m_words.resize(on);
                std::copy_n(other.m_words.data(), on, m_words.data());
                // in-place: *this -= old_this
                mpn::sub(m_words.data(), m_words.data(), on, tbuf, tn);
                if (tbuf != stack_buf) delete[] tbuf;
                // Normalize
                while (on > 1 && m_words.data()[on - 1] == 0) on--;
                m_words.resize(on);
                m_sign = other.m_sign;
                if (on == 1 && m_words.data()[0] == 0) m_sign = 0;
            } else {
                IntOps::subtractAbsoluteInPlace(*this, other);
            }
        }
        return *this;
    }

    Int& Int::operator-=(const Int& other) {
        if (isSpecialState() || other.isSpecialState()) {
            *this = IntSpecialStates::handleSubtraction(*this, other);
            return *this;
        }
        if (other.getSign() == 0) return *this;
        if (getSign() == 0) {
            *this = other;
            if (m_sign != 0) m_sign = -m_sign;
            return *this;
        }

        int rhsEffSign = -other.m_sign;
        if (m_sign == rhsEffSign) {
            // Same sign: add the magnitudes (e.g. 5-(-3)=8)
            IntOps::addAbsolute(*this, other);
        } else {
            // Different sign: subtract the magnitudes (e.g. 5-3=2)
            if (IntOps::compareAbsLess(*this, other)) {
                Int result;
                IntOps::subtractAbsolute(other, *this, result);
                result.setSign(rhsEffSign);
                *this = std::move(result);
            } else {
                IntOps::subtractAbsoluteInPlace(*this, other);
            }
        }
        return *this;
    }

    Int& Int::operator*=(const Int& other) {
        *this = *this * other;
        return *this;
    }

    Int& Int::operator/=(const Int& other) {
        *this = *this / other;
        return *this;
    }

    Int& Int::operator%=(const Int& other) {
        *this = *this % other;
        return *this;
    }

    Int& Int::operator&=(const Int& other) {
        *this = *this & other;
        return *this;
    }

    Int& Int::operator|=(const Int& other) {
        *this = *this | other;
        return *this;
    }

    Int& Int::operator^=(const Int& other) {
        *this = pow(*this, other);
        return *this;
    }

    // setState implementation
    void Int::setState(NumericState state, NumericError error) {
        m_state = state;
        m_error = error;

        // Overwrite the sign only for special states (no branch on the Normal path)
        if (state != NumericState::Normal) [[unlikely]] {
            if (state == NumericState::NaN) {
                m_sign = 0;
            }
            else if (state == NumericState::PositiveInfinity) {
                m_sign = 1;
            }
            else if (state == NumericState::NegativeInfinity) {
                m_sign = -1;
            }
        }
    }

    // Factory methods for special states (NaN, etc.)
    Int Int::NaN() {
        Int result;
        result.m_state = NumericState::NaN;
        result.m_error = NumericError::ExplicitNaN;
        return result;
    }

    Int Int::PositiveInfinity() {
        Int result;
        result.m_state = NumericState::PositiveInfinity;
        result.m_sign = 1;
        return result;
    }

    Int Int::NegativeInfinity() {
        Int result;
        result.m_state = NumericState::NegativeInfinity;
        result.m_sign = -1;
        return result;
    }

    // Internal method
    void Int::setNaN(NumericError error) {
        m_state = NumericState::NaN;
        m_error = error;
        m_sign = 0; // The sign of NaN is normally treated as zero
    }

    // Implementation of digitCount
    size_t Int::digitCount(int base) const {
        // Validate the base
        if (base < 2 || base > 36) {
            return 0;  // Invalid base
        }
        // For special states
        if (isSpecialState()) return 0;
        // For 0, the digit count is 1
        if (isZero()) return 1;

        // Take the absolute value
        Int abs_val = abs(*this);

        // Optimize for power-of-two bases
        if (base == 2) {
            return abs_val.bitLength();
        } else if (base == 16) {
            return (abs_val.bitLength() + 3) / 4;  // 4 bits = 1 digit
        } else if (base == 8) {
            return (abs_val.bitLength() + 2) / 3;  // 3 bits ≈ 1 digit
        }

        // General base: convert to a string and count the digits
        std::string str = abs_val.toString(base);
        return str.length();
    }

    // Implementation of sizeInBase (compatible with GMP mpz_sizeinbase)
    size_t Int::sizeInBase(int base) const {
        // Validate the base
        if (base < 2 || base > 62) return 0;
        // For special states
        if (isSpecialState()) return 0;
        // For 0, the digit count is 1
        if (isZero()) return 1;

        size_t bits = bitLength();

        // For power-of-two bases: an exact result
        if ((base & (base - 1)) == 0) {
            int bits_per_digit = 0;
            int b = base;
            while (b > 1) { b >>= 1; bits_per_digit++; }
            return (bits + bits_per_digit - 1) / bits_per_digit;
        }

        // General base: estimate via a double-precision logarithm
        // The result is either exact or larger by 1 (same semantics as GMP)
        double log_base = std::log(static_cast<double>(base));
        size_t result = static_cast<size_t>(bits * std::log(2.0) / log_base) + 1;
        return result;
    }

    // Implementation of ilog2
    size_t Int::ilog2() const {
        // 0 for special states, zero, or negative values
        if (isSpecialState() || !isPositive()) return 0;
        // floor(log2(n)) = bitLength - 1
        return bitLength() - 1;
    }

    // Implementation of ilog10
    size_t Int::ilog10() const {
        // 0 for special states, zero, or negative values
        if (isSpecialState() || !isPositive()) return 0;
        // digitCount(10) returns the exact number of decimal digits
        // floor(log10(n)) = digitCount - 1
        return digitCount(10) - 1;
    }

    // Implementation of isDivisible
    bool Int::isDivisible(const Int& divisor) const {
        // Handle special states
        if (isSpecialState() || divisor.isSpecialState()) return false;
        // Testing divisibility by 0 is meaningless
        if (divisor.isZero()) return false;
        // 0 is divisible by any number
        if (isZero()) return true;
        // Decide via this % divisor == 0
        return (*this % divisor).isZero();
    }

    // Implementation of isCongruent
    bool Int::isCongruent(const Int& other, const Int& modulus) const {
        // Handle special states
        if (isSpecialState() || other.isSpecialState() || modulus.isSpecialState()) return false;
        // If modulus is 0, the test is undefined
        if (modulus.isZero()) return false;
        // Decide via (this - other) % modulus == 0
        return (*this - other).isDivisible(modulus);
    }

} // namespace sangi