#include #include #include #include #include #include #include #include #include #include #include #include #include namespace sangi {

// Forward declarations (used by operator+= / -=)
static int mergeRequested(int req_a, int req_b);
static int mergeEffective(int eff_a, int eff_b);

// Compute 5^n by repeated squaring (O(M(n) log n))
// Since 10^n = 5^n * 2^n, handle 2^n with a shift and only compute 5^n
static Int pow5_fast(size_t n) {
    if (n == 0) return Int(1);
    if (n == 1) return Int(5);
    Int base(5);
    Int result(1);
    size_t exp = n;
    while (exp > 0) {
        if (exp & 1) result = result * base;
        exp >>= 1;
        if (exp > 0) base = base * base;
    }
    return result;
}

// Compute 10^n by repeated squaring (O(M(n) log n))
[[maybe_unused]] static Int pow10_fast(size_t n) {
    if (n == 0) return Int(1);
    if (n == 1) return Int(10);
    Int base(10);
    Int result(1);
    size_t exp = n;
    while (exp > 0) {
        if (exp & 1) result = result * base;
        exp >>= 1;
        if (exp > 0) base = base * base;
    }
    return result;
}



// defaultPrecision() and rounding_mode_ are already defined as inline static thread_local in Float.hpp.
// The variable exists in the caller's TU even across DLL boundaries, so no export is needed.
thread_local int64_t Float::emin_ = Float::EXPONENT_MIN;
thread_local int64_t Float::emax_ = Float::EXPONENT_MAX;
thread_local unsigned Float::exception_flags_ = 0;

void Float::checkExponentBounds() {
    if (is_infinity_ || is_nan_ || mantissa_.isZero()) {
        return;
    }
    if (exponent_ > EXPONENT_MAX) {
        // Exponent overflow → ±∞
        is_infinity_ = true;
        effective_bits_ = 0;
        requested_bits_ = 0;
        mantissa_ = Int(0);
        exponent_ = 0;
    } else if (exponent_ < EXPONENT_MIN) {
        // Exponent underflow → 0
        mantissa_ = Int(0);
        exponent_ = 0;
        is_negative_ = false;
        effective_bits_ = 0;
        requested_bits_ = 0;
    }
}

// Convert decimal digit count to bit count (log₂(10) ≈ 3.32)
int Float::precisionToBits(int precision) {
    // To accurately guarantee N decimal digits, in addition to N * log2(10) bits
    // guard bits are required. Without guard bits the final 1-2 digits become inaccurate due to rounding error.
    // 8 guard bits ≈ 2.4 digits protect the trailing digit.
    constexpr int GUARD = 8;
    return static_cast (std::ceil(precision * 3.32192809488736)) + GUARD;
}

// Convert bit count to decimal digit count
int Float::bitsToPrecision(int bits) {
    if (bits <= 0) return 1;

    // Subtract the guard bits added by precisionToBits before converting to digit count
    // Since precisionToBits(N) = ceil(N * log2(10)) + GUARD,
    // the inverse conversion subtracts GUARD
    constexpr int GUARD = 8;  // same value as in precisionToBits
    int effective_bits = std::max(1, bits - GUARD);
    int precision = static_cast (std::floor(effective_bits / 3.32192809488736));

    return std::max(1, precision);
}

//======================================================================
// Constructors
//======================================================================

Float::Float()
    : mantissa_(0), exponent_(0), is_negative_(false),
      is_infinity_(false), is_nan_(false),
      effective_bits_(INT_MAX), requested_bits_(INT_MAX) {
}

Float::Float(int value)
    : mantissa_(std::abs(value)), exponent_(0), is_negative_(value < 0),
      is_infinity_(false), is_nan_(false),
      effective_bits_(INT_MAX), requested_bits_(INT_MAX) {
    normalize();
}

Float::Float(int64_t value)
    : mantissa_(std::abs(value)), exponent_(0), is_negative_(value < 0),
      is_infinity_(false), is_nan_(false),
      effective_bits_(INT_MAX), requested_bits_(INT_MAX) {
    normalize();
}

Float::Float(int64_t mantissa, int64_t exponent, bool is_negative)
    : mantissa_(std::abs(mantissa)), exponent_(exponent), is_negative_(mantissa < 0 ? !is_negative : is_negative),
    is_infinity_(false), is_nan_(false),
    effective_bits_(INT_MAX), requested_bits_(INT_MAX) {
    // Flip the sign if a negative value was passed
    if (mantissa < 0) {
        mantissa_ = Int(-mantissa);
    }
    normalize();
}

Float::Float(const Int& value)
    : mantissa_(abs(value)), exponent_(0), is_negative_(value.isNegative()),
      is_infinity_(false), is_nan_(false),
      effective_bits_(INT_MAX), requested_bits_(INT_MAX) {
    normalize();
}

Float::Float(double value) : is_infinity_(false), is_nan_(false),
    effective_bits_(53), requested_bits_(precisionToBits(defaultPrecision())) {
    // Handle special values
    if (std::isnan(value)) {
        is_nan_ = true;
        is_negative_ = false;
        mantissa_ = Int(0);
        exponent_ = 0;
        effective_bits_ = 0;
        requested_bits_ = 0;
        return;
    }

    is_negative_ = std::signbit(value);
    value = std::abs(value);

    if (std::isinf(value)) {
        is_infinity_ = true;
        mantissa_ = Int(0);
        exponent_ = 0;
        effective_bits_ = 0;
        requested_bits_ = 0;
        return;
    }

    // Handle zero
    if (value == 0.0) {
        mantissa_ = Int(0);
        exponent_ = 0;
        effective_bits_ = INT_MAX;
        requested_bits_ = INT_MAX;
        return;
    }

    // Handle the double directly (without using frexp)
    // Obtain the bit representation of the floating-point number
    uint64_t bits;
    std::memcpy(&bits, &value, sizeof(double));

    // Extraction from IEEE 754 double-precision format
    // sign bit (1bit), exponent (11bits), mantissa (52bits)
    uint64_t exponent_bits = (bits >> 52) & 0x7FF;
    uint64_t mantissa_bits = bits & 0x000FFFFFFFFFFFFF;

    // Biased exponent (1023 in IEEE 754)
    int64_t true_exponent = static_cast (exponent_bits) - 1023;

    // Add the implicit leading bit (1) to the mantissa
    mantissa_ = Int(mantissa_bits | 0x10000000000000); // add 2^52

    // Adjust the exponent (since the mantissa is multiplied by 2^52)
    exponent_ = true_exponent - 52;

    // Determining the precision fields: binary exponent analysis
    // value = full_mantissa × 2^(true_exponent - 52)
    // Remove trailing zeros from the mantissa to get the form odd_part × 2^e_adj
    // If e_adj >= 0 it is an integer → exact (INT_MAX)
    // If e_adj < 0 and |e_adj| <= 20 it is a clean decimal fraction → exact (INT_MAX)
    if (value == 0.0) {
        effective_bits_ = INT_MAX;
        requested_bits_ = INT_MAX;
    } else if (exponent_bits != 0x7FF) {  // other than inf, NaN
        uint64_t m = mantissa_bits | 0x10000000000000;
        int trailing_zeros = 0;
        while ((m & 1) == 0) {
            m >>= 1;
            trailing_zeros++;
        }
        int64_t e_adj = true_exponent - 52 + trailing_zeros;
        if ((e_adj >= 0) || ((-e_adj) <= 20)) {
            effective_bits_ = INT_MAX;
            requested_bits_ = INT_MAX;
        }
        // otherwise: effective_bits_ = 53, requested_bits_ = default (already set in the initializer list)
    }

    normalize();
}

Float::Float(std::string_view str)
    : is_infinity_(false), is_nan_(false),
      effective_bits_(0), requested_bits_(0) {
    initFromDecimalString(str, -1);
}

Float::Float(std::string_view str, int precision)
    : is_infinity_(false), is_nan_(false),
      effective_bits_(0), requested_bits_(0) {
    initFromDecimalString(str, precision);
}

// If requested_precision >= 0, round non-exact values to that decimal precision.
// If requested_precision < 0 (default), use max(string significant digits, defaultPrecision()).
void Float::initFromDecimalString(std::string_view str, int requested_precision) {
    // Check for special values
    std::string s(str);
    std::transform(s.begin(), s.end(), s.begin(), [](unsigned char c) { return std::tolower(c); });

    if (s == "nan" || s == "inf" || s == "+inf" || s == "-inf" || s == "infinity" || s == "+infinity" || s == "-infinity") {
        if (s == "nan") {
            is_nan_ = true;
            is_negative_ = false;
        } else {
            is_infinity_ = true;
            is_negative_ = (s[0] == '-');
        }
        mantissa_ = 0;
        exponent_ = 0;
        return;
    }

    // High-precision decimal string parsing
    // Handle the sign
    size_t pos = 0;
    is_negative_ = false;
    if (pos < str.size() && (str[pos] == '+' || str[pos] == '-')) {
        is_negative_ = (str[pos] == '-');
        pos++;
    }

    // Collect the digits of the integer and fractional parts
    std::string digits;
    int decimal_digits = 0;
    bool found_dot = false;
    int exponent_10 = 0;

    for (; pos < str.size(); pos++) {
        if (str[pos] == '.') {
            found_dot = true;
        } else if (str[pos] == 'e' || str[pos] == 'E') {
            // Exponent part
            pos++;
            exponent_10 = std::stoi(std::string(str.substr(pos)));
            break;
        } else if (str[pos] >= '0' && str[pos] <= '9') {
            digits += str[pos];
            if (found_dot) decimal_digits++;
        } else {
            break;
        }
    }

    if (digits.empty()) {
        // Parse failure
        is_nan_ = true;
        mantissa_ = 0;
        exponent_ = 0;
        is_negative_ = false;
        return;
    }

    // Remove leading zeros
    size_t first_nonzero = digits.find_first_not_of('0');
    if (first_nonzero == std::string::npos) {
        // All zeros
        mantissa_ = Int(0);
        exponent_ = 0;
        is_negative_ = false;
        effective_bits_ = INT_MAX;
        requested_bits_ = INT_MAX;
        return;
    }

    // Compute the precision bit count from the significant digit count
    int significant_digits = static_cast (digits.size() - first_nonzero);
    int string_precision_bits;

    if (requested_precision >= 0) {
        // Two-argument version: round non-exact values to the explicitly specified decimal precision.
        // If the string has more significant digits than that, prefer the string side (do not discard information).
        string_precision_bits = std::max(precisionToBits(requested_precision),
                                         precisionToBits(significant_digits));
    } else {
        // One-argument version: guarantee at least defaultPrecision even for non-exact strings
        // (max(string, default)).
        // Motivation: addressing the problem where effective_bits of Float("0.1") was interpreted
        //       as "1 digit = 12 bits", causing precision degradation in LinAlg use.
        //   - Float("0.1") is mathematically 1/10; the user expects "0.1 at full precision",
        //     not "a 3-digit value"
        //   - to explicitly state a precision constraint, use the two-argument version or setPrecision
        //   - integer/exact strings are INT_MAX, so this branch does not affect them
        string_precision_bits = precisionToBits(significant_digits);
        int default_precision_bits = precisionToBits(defaultPrecision());
        if (default_precision_bits > string_precision_bits) {
            string_precision_bits = default_precision_bits;
        }
    }

    // Convert the decimal digit string to an Int
    Int numerator(digits);

    // value = numerator * 10^(exponent_10 - decimal_digits)
    int net_exp10 = exponent_10 - decimal_digits;

    // "Mathematically exact value" determination:
    //   - integer strings (e.g. "1", "100", "1e3", net_exp10 >= 0) are exact
    //   - pure binary finite fractions (e.g. "0.5", "0.25", net_exp10 < 0 where 5^D divides
    //     numerator) are also exact
    // When judged exact, set effective_bits_ = INT_MAX to preserve arithmetic precision.
    // This ensures that in the typical case of using Float("1") as input to LinAlg etc., later
    // operations (= 1/5 etc.) are executed at default precision.
    bool is_exact = false;

    if (net_exp10 == 0) {
        // value = numerator (integer) — always exact
        mantissa_ = numerator;
        exponent_ = 0;
        is_exact = true;
        normalize();
    } else if (net_exp10 > 0) {
        // value = numerator * 10^net_exp10 = numerator * 5^e * 2^e — integer and exact
        Int five_pow = pow5_fast(static_cast (net_exp10));
        mantissa_ = numerator * five_pow;
        exponent_ = net_exp10; // 2^net_exp10
        is_exact = true;
        normalize();
    }

    if (is_exact) {
        effective_bits_ = INT_MAX;
        requested_bits_ = INT_MAX;
    } else if (net_exp10 < 0) {
        // value = numerator / 10^D = numerator / (5^D * 2^D)
        // Divide by 5^D, and handle 2^D via exponent adjustment
        size_t D = static_cast (-net_exp10);
        Int five_pow = pow5_fast(D);

        // Pure binary finite-fraction test: exact if numerator is divisible by 5^D.
        // e.g. "0.5" (5/10 = 1/2), "0.25" (25/100 = 1/4), "1.5" (15/10 = 3/2)
        // e.g. "0.1" (1/10) is not exact, since 5^1=5 does not divide 1
        Int remainder_check = numerator % five_pow;
        if (remainder_check.isZero()) {
            // Perfect binary finite fraction → exact
            mantissa_ = numerator / five_pow;
            exponent_ = -static_cast (D);
            effective_bits_ = INT_MAX;
            requested_bits_ = INT_MAX;
            normalize();
        } else {
            int num_bits = static_cast (numerator.bitLength());
            int den_bits = static_cast (five_pow.bitLength());
            // Set so the quotient secures string_precision_bits + guard bits
            // quotient_bits ≈ num_bits + target_bits - den_bits ≥ string_precision_bits + 20
            int target_bits = std::max(den_bits - num_bits, 0) + string_precision_bits + 20;

            Int scaled = numerator << target_bits;
            Int remainder;
            Int quotient = IntOps::divmod(scaled, five_pow, remainder);

            mantissa_ = std::move(quotient);
            exponent_ = -static_cast (target_bits) - static_cast (D);
            effective_bits_ = string_precision_bits;
            requested_bits_ = string_precision_bits;
            normalize();
        }
    }
}

Float::Float(const Int& mantissa, int64_t exponent, bool is_negative)
    : mantissa_(mantissa), exponent_(exponent), is_negative_(is_negative),
      is_infinity_(false), is_nan_(false),
      effective_bits_(INT_MAX), requested_bits_(INT_MAX) {
    normalize();
}

Float::Float(Int&& mantissa, int64_t exponent, bool is_negative)
    : mantissa_(std::move(mantissa)), exponent_(exponent), is_negative_(is_negative),
      is_infinity_(false), is_nan_(false),
      effective_bits_(INT_MAX), requested_bits_(INT_MAX) {
    normalize();
}

//======================================================================
// Assignment operators
//======================================================================

Float& Float::operator=(int64_t value) {
    mantissa_ = Int(std::abs(value));
    exponent_ = 0;
    is_negative_ = (value < 0);
    is_infinity_ = false;
    is_nan_ = false;
    effective_bits_ = INT_MAX;
    requested_bits_ = INT_MAX;
    normalize();
    return *this;
}

Float& Float::operator=(double value) {
    *this = Float(value);
    return *this;
}

//======================================================================
// Operator implementations
//======================================================================

Float& Float::operator+=(const Float& rhs) {
    // Delegate special values to the legacy path
    if (isNaN() || rhs.isNaN()) [[unlikely]] { *this = Float::nan(); return *this; }
    if (isInfinity() || rhs.isInfinity()) [[unlikely]] {
        *this = static_cast (*this) + rhs; return *this;
    }
    if (isZero()) [[unlikely]] { *this = rhs; return *this; }
    if (rhs.isZero()) [[unlikely]] return *this;

    // Read the precision fields first (before this is overwritten)
    int eff = mergeEffective(effective_bits_, rhs.effective_bits_);
    int req = mergeRequested(requested_bits_, rhs.requested_bits_);

    if (is_negative_ == rhs.is_negative_) {
        // Same sign: mantissa addition (no cancellation)
        addUnsignedInPlace(rhs);
        effective_bits_ = eff;
        requested_bits_ = req;
    } else {
        // Opposite signs: legacy path, since cancellation computation is complex
        int64_t mag_this = static_cast (mantissa_.bitLength()) + exponent_;
        int64_t mag_rhs = static_cast (rhs.mantissa_.bitLength()) + rhs.exponent_;
        int64_t max_mag = std::max(mag_this, mag_rhs);
        bool this_neg = is_negative_;

        subtractUnsignedInPlace(rhs);

        if (isZero()) {
            if (eff >= INT_MAX) {
                effective_bits_ = INT_MAX;
                requested_bits_ = INT_MAX;
            } else {
                effective_bits_ = 0;
                requested_bits_ = req;
            }
        } else {
            is_negative_ = is_negative_ ? rhs.isNegative() : this_neg;
            int64_t result_mag = static_cast (mantissa_.bitLength()) + exponent_;
            int lost_bits = static_cast (std::max(int64_t(0), max_mag - result_mag));
            if (eff >= INT_MAX) {
                effective_bits_ = INT_MAX;
            } else {
                effective_bits_ = std::max(1, eff - lost_bits);
            }
            requested_bits_ = req;
        }
    }
    return *this;
}

Float& Float::operator-=(const Float& rhs) {
    // Delegate special values to the legacy path
    if (isNaN() || rhs.isNaN()) [[unlikely]] { *this = Float::nan(); return *this; }
    if (isInfinity() || rhs.isInfinity()) [[unlikely]] {
        *this = static_cast (*this) - rhs; return *this;
    }
    if (isZero()) [[unlikely]] {
        *this = rhs;
        is_negative_ = !is_negative_;
        return *this;
    }
    if (rhs.isZero()) [[unlikely]] return *this;

    int eff = mergeEffective(effective_bits_, rhs.effective_bits_);
    int req = mergeRequested(requested_bits_, rhs.requested_bits_);

    if (is_negative_ != rhs.is_negative_) {
        // Opposite signs: (+a)-(-b)=a+b → mantissa addition (no cancellation)
        addUnsignedInPlace(rhs);
        effective_bits_ = eff;
        requested_bits_ = req;
    } else {
        // Same sign: mantissa subtraction (with cancellation)
        int64_t mag_this = static_cast (mantissa_.bitLength()) + exponent_;
        int64_t mag_rhs = static_cast (rhs.mantissa_.bitLength()) + rhs.exponent_;
        int64_t max_mag = std::max(mag_this, mag_rhs);
        bool this_neg = is_negative_;

        subtractUnsignedInPlace(rhs);

        if (isZero()) {
            if (eff >= INT_MAX) {
                effective_bits_ = INT_MAX;
                requested_bits_ = INT_MAX;
            } else {
                effective_bits_ = 0;
                requested_bits_ = req;
            }
        } else {
            // subtractUnsigned: is_negative_=true means |this| < |rhs|
            // operator-: same sign with |this|>=|rhs| → sign=this_neg, |this|<|rhs| → sign=!this_neg
            is_negative_ = is_negative_ ? !this_neg : this_neg;
            int64_t result_mag = static_cast (mantissa_.bitLength()) + exponent_;
            int lost_bits = static_cast (std::max(int64_t(0), max_mag - result_mag));
            if (eff >= INT_MAX) {
                effective_bits_ = INT_MAX;
            } else {
                effective_bits_ = std::max(1, eff - lost_bits);
            }
            requested_bits_ = req;
        }
    }
    return *this;
}

Float& Float::operator*=(const Float& rhs) {
    *this = std::move(*this) * rhs;
    return *this;
}

Float& Float::operator/=(const Float& rhs) {
    *this = std::move(*this) / rhs;
    return *this;
}

//======================================================================
// Method implementations
//======================================================================

bool Float::isZero() const {
    return !is_infinity_ && !is_nan_ && mantissa_.isZero();
}

int Float::precision() const {
    if (isZero() || isInfinity() || isNaN()) {
        return defaultPrecision();
    }
    int eff = std::min(effective_bits_, requested_bits_);
    if (eff >= INT_MAX) {
        // exact value (integer-derived) → defaultPrecision compatible
        return defaultPrecision();
    }
    return Float::bitsToPrecision(eff);
}

Float& Float::setPrecision(int precision) {
    if (isZero() || isInfinity() || isNaN()) {
        return *this;
    }

    int precision_bits = precisionToBits(precision);
    int current_bits = static_cast (mantissa_.bitLength());

    // Update requested_bits_
    requested_bits_ = precision_bits;

    if (current_bits < precision_bits) {
        // If precision is insufficient, left-shift the mantissa to pad
        // The value does not change: (mantissa << shift) * 2^(exponent - shift) = mantissa * 2^exponent
        // Do not change effective_bits_ (padding does not increase reliability)
        int shift = precision_bits - current_bits;
        mantissa_ <<= shift;
        exponent_ -= shift;
        checkExponentBounds();
    }

    round(precision_bits, roundingMode());
    return *this;
}

Float& Float::setResultPrecision(int precision) {
    setPrecision(precision);
    if (!isZero() && !isInfinity() && !isNaN()) {
        effective_bits_ = precisionToBits(precision);
    }
    return *this;
}

void Float::truncateToApprox(int precision) {
    if (isZero() || isInfinity() || isNaN()) [[unlikely]] return;

    int target_bits = precisionToBits(precision);
    requested_bits_ = target_bits;
    int bl = static_cast (mantissa_.bitLength());

    if (bl < target_bits) {
        // Padding: extend the mantissa as in setPrecision
        int shift = target_bits - bl;
        mantissa_ <<= shift;
        exponent_ -= shift;
        checkExponentBounds();
        return;
    }

    if (bl <= target_bits) return;

    // Truncate the low part word-by-word (no bit shift or rounding)
    // floor division: keep up to 63 bits of excess → absorbed by guard bits
    int words_to_drop = (bl - target_bits) / 64;
    if (words_to_drop <= 0) return;

    int shift = words_to_drop * 64;
    mantissa_ >>= shift;  // bit_shift==0 → erase only (O(1) memmove)
    exponent_ += shift;
}

void Float::normalize() {
    if (is_infinity_ || is_nan_ || mantissa_.isZero()) [[unlikely]] {
        return;
    }

    // Remove zero words at both ends in-place (no alloc)
    size_t lsb_removed = mantissa_.trimZeroWords();
    if (lsb_removed > 0) {
        exponent_ += static_cast (64 * lsb_removed);
    }

    if (mantissa_.isZero()) [[unlikely]] {
        exponent_ = 0;
        is_negative_ = false;
    }

    checkExponentBounds();
}

void Float::round(int precision_bits, RoundingMode mode) {
    if (isZero() || isInfinity() || isNaN()) {
        return;
    }

    int bit_length = static_cast (mantissa_.bitLength());
    if (bit_length <= precision_bits) {
        return;  // already at or below the required precision
    }

    int shift = bit_length - precision_bits;

    // ── Extract guard and sticky bits directly from the raw limbs ──
    // Avoid Int copies and mask generation; make the rounding decision with O(n) reads only.
    const uint64_t* data = mantissa_.data();
    size_t n = mantissa_.size();

    // guard bit = bit (shift - 1): the bit just below the precision boundary
    bool guard_bit = false;
    if (shift >= 1) {
        size_t gw = static_cast (shift - 1) / 64;
        unsigned gb = static_cast ((shift - 1) % 64);
        if (gw < n) {
            guard_bit = ((data[gw] >> gb) & 1) != 0;
        }
    }

    // sticky bits = bits 0..(shift-2): true if any of them is non-zero
    bool sticky = false;
    if (shift >= 2) {
        size_t full_words = static_cast (shift - 1) / 64;
        for (size_t i = 0; i < full_words && i < n; ++i) {
            if (data[i] != 0) { sticky = true; break; }
        }
        if (!sticky && full_words < n) {
            unsigned partial = static_cast ((shift - 1) % 64);
            if (partial > 0) {
                uint64_t mask = (1ULL << partial) - 1;
                if (data[full_words] & mask) sticky = true;
            }
        }
    }

    // Remove the low bits
    mantissa_ >>= shift;
    exponent_ += shift;

    // Determine the rounding direction
    bool any_discarded = guard_bit || sticky;
    bool round_up = false;

    switch (mode) {
        case RoundingMode::ToNearest:
            // Round to even: round up when guard=1 and (sticky=1 or result LSB=1)
            if (guard_bit) {
                round_up = sticky || ((mantissa_.word(0) & 1) != 0);
            }
            break;

        case RoundingMode::TowardZero:
            break;

        case RoundingMode::TowardPositive:
            round_up = !is_negative_ && any_discarded;
            break;

        case RoundingMode::TowardNegative:
            round_up = is_negative_ && any_discarded;
            break;

        case RoundingMode::AwayFromZero:
            round_up = any_discarded;
            break;
    }

    if (round_up) {
        IntOps::addDelta(mantissa_, 1);
        if (static_cast (mantissa_.bitLength()) > precision_bits) {
            mantissa_ >>= 1;
            exponent_ += 1;
        }
    }

    // Since we truncated during rounding, limit effective_bits_ to at most precision_bits
    if (effective_bits_ > precision_bits && effective_bits_ < INT_MAX) {
        effective_bits_ = precision_bits;
    }

    // MIN policy: rounding clamps effective_bits_ to at most precision_bits
    if constexpr (PRECISION_POLICY == PrecisionPolicy::MIN_PROPAGATION) {
        if (effective_bits_ < INT_MAX) {
            effective_bits_ = std::min(effective_bits_, precision_bits);
        }
    }

    // Range check the exponent
    checkExponentBounds();
}

void Float::subnormalize(RoundingMode mode) {
    // Do nothing for special values and zero
    if (isZero() || isInfinity() || isNaN()) return;

    // Do nothing if within the normal range
    // IEEE 754: the minimum exponent of normal numbers is emin; with p mantissa bits,
    // subnormal condition: exponent < emin + p - 1 (= the leading mantissa bit does not reach emin)
    int bit_length = static_cast (mantissa_.bitLength());
    // effective exponent = exponent_ + bit_length - 1 (position of MSB)
    int64_t msb_exp = exponent_ + static_cast (bit_length) - 1;

    if (msb_exp >= emin_) return;  // normal number

    // Shift amount: shift the mantissa right to align with emin
    int64_t shift = emin_ - msb_exp;

    if (shift >= bit_length) {
        // Complete underflow → zero
        mantissa_ = Int(0);
        exponent_ = 0;
        is_negative_ = false;
        effective_bits_ = 0;
        requested_bits_ = 0;
        raiseException(FE_UNDERFLOW | FE_INEXACT);
        return;
    }

    // Extract guard/sticky bits for rounding
    int shift_int = static_cast (shift);
    bool guard_bit = mantissa_.getBit(shift_int - 1);
    bool sticky = false;
    for (int i = 0; i < shift_int - 1; ++i) {
        if (mantissa_.getBit(i)) { sticky = true; break; }
    }
    bool any_discarded = guard_bit || sticky;

    // Perform the shift
    mantissa_ >>= shift_int;
    exponent_ += shift_int;  // = emin - (bit_length - 1) + shift = emin

    // Rounding decision
    bool round_up = false;
    switch (mode) {
        case RoundingMode::ToNearest:
            if (guard_bit) {
                round_up = sticky || ((mantissa_.word(0) & 1) != 0);
            }
            break;
        case RoundingMode::TowardZero:
            break;
        case RoundingMode::TowardPositive:
            round_up = !is_negative_ && any_discarded;
            break;
        case RoundingMode::TowardNegative:
            round_up = is_negative_ && any_discarded;
            break;
        case RoundingMode::AwayFromZero:
            round_up = any_discarded;
            break;
    }

    if (round_up) {
        IntOps::addDelta(mantissa_, 1);
    }

    if (mantissa_.isZero()) {
        exponent_ = 0;
        is_negative_ = false;
        effective_bits_ = 0;
        requested_bits_ = 0;
    }

    if (any_discarded) {
        raiseException(FE_UNDERFLOW | FE_INEXACT);
    }
}

std::string Float::toString(int precision) const {
    if (isNaN()) {
        return "NaN";
    }
    
    if (isInfinity()) {
        return is_negative_ ? "-Infinity" : "Infinity";
    }
    
    if (isZero()) {
        return "0.0";
    }

    // By default use the current precision
    if (precision < 0) {
        precision = this->precision();
    }

    // Decide the display format based on the size of the exponent
    if (exponent_ >= 0 && exponent_ < precision) {
        return toDecimalString(precision);
    } else {
        return toScientificString(precision);
    }
}

std::string Float::toDecimalString(int precision) const {
    if (isNaN()) return "NaN";
    if (isInfinity()) return is_negative_ ? "-Infinity" : "Infinity";
    if (isZero()) return "0.0";

    if (precision < 0) {
        precision = std::max(1, this->precision());
    }

    // value = mantissa * 2^exponent  (mantissa >= 0)
    // |value| * 10^D = mantissa * 10^D * 2^exponent
    //                = mantissa * (2^D * 5^D) * 2^exponent
    //                = mantissa * 5^D * 2^(D + exponent)
    //
    // shift = D + exponent
    //   shift >= 0: N = (mantissa * 5^D) << shift    (exact)
    //   shift <  0: N = (mantissa * 5^D) >> (-shift)  (truncation)
    //
    // Convert N = floor(|value| * 10^D) to a decimal string and insert the decimal point.

    int D = precision;

    // Guard bits: to prevent the trailing digit from becoming inaccurate due to
    // rounding in the binary→decimal conversion (truncation of the right shift), compute with extra bits added to the mantissa.
    constexpr int GUARD_BITS = 20;  // 20 bits ≈ 6 digits of guard
    Int extended_mantissa = mantissa_ << GUARD_BITS;
    int64_t adjusted_exponent = exponent_ - GUARD_BITS;

    Int pow5 = pow(Int(5), static_cast (D));
    Int product = extended_mantissa * pow5;

    int64_t shift = static_cast (D) + adjusted_exponent;

    Int N;
    if (shift >= 0) {
        N = product << static_cast (shift);
    } else {
        int rshift = static_cast (-shift);
        N = product >> rshift;
        // Round half up: carry if the MSB of the discarded bit string is 1
        if (rshift >= 1 && product.getBit(rshift - 1)) {
            N += 1;
        }
    }

    // N = round(|value| * 10^D)
    std::string digits = N.toString();

    // Compute the number of integer-part digits independently from the magnitude of |value|
    // (the digit count of N can be larger than expected due to binary rounding)
    int intDigits;
    {
        double val = std::abs(toDouble());
        if (val >= 1.0)
            intDigits = static_cast