llvm-native-core 0.1.16

// llvm-native support module — APFloat: Arbitrary-precision IEEE 754 floating-point.
// Clean-room implementation from IEEE 754-2019 standard and published algorithms
// (Goldberg's "What Every Computer Scientist Should Know About Floating-Point
// Arithmetic", Clinger's "How to Read Floating Point Numbers Accurately", etc.).
// Phase 1 — LLVM.SUPPORT.1 Court
//
// No LLVM C++ source code has been consulted.

use crate::support::apint::APInt;
use std::cmp::Ordering;
use std::fmt;

// ============================================================================
// Exception Flags
// ============================================================================

/// No exceptions raised.
pub const APFLOAT_OK: u32 = 0;
/// Invalid operation (e.g., 0/0, ∞-∞, sqrt(-1)).
pub const APFLOAT_INVALID_OP: u32 = 1 << 0;
/// Division by zero.
pub const APFLOAT_DIV_BY_ZERO: u32 = 1 << 1;
/// Overflow: result too large in magnitude to be represented.
pub const APFLOAT_OVERFLOW: u32 = 1 << 2;
/// Underflow: result too small in magnitude to be represented (may be inexact).
pub const APFLOAT_UNDERFLOW: u32 = 1 << 3;
/// Inexact: result had to be rounded.
pub const APFLOAT_INEXACT: u32 = 1 << 4;

// ============================================================================
// Rounding Modes
// ============================================================================

/// IEEE 754 rounding direction.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum RoundingMode {
    /// Round to nearest, ties to even (default IEEE 754).
    NearestTiesToEven,
    /// Round toward zero (truncation).
    TowardZero,
    /// Round toward negative infinity (floor).
    TowardNegative,
    /// Round toward positive infinity (ceil).
    TowardPositive,
    /// Round to nearest, ties away from zero.
    NearestTiesToAway,
}

// ============================================================================
// Float Category
// ============================================================================

/// IEEE 754 classification of a floating-point value.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum FloatCategory {
    /// Positive or negative zero.
    Zero,
    /// Normal non-zero finite value.
    Normal,
    /// Positive or negative infinity.
    Infinity,
    /// Quiet or signaling NaN.
    NaN,
    /// Subnormal (denormalized) value.
    Subnormal,
}

// ============================================================================
// FloatSemantics
// ============================================================================

/// Describes the IEEE 754 format for a floating-point type.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct FloatSemantics {
    /// Human-readable name for this format.
    pub name: &'static str,

    /// Total bit width of the format (32, 64, 80, 128).
    pub bits: u32,

    /// Width of the exponent field in bits.
    pub exponent_bits: u32,

    /// Width of the explicit mantissa (fraction) field.
    /// For IEEE formats this is `significand_bits - 1` (the implicit leading
    /// bit is not counted). For x87 this includes the explicit leading bit.
    pub mantissa_bits: u32,

    /// Bias added to the unbiased exponent to produce the stored exponent.
    pub exponent_bias: i32,

    /// Maximum raw exponent value (all ones in the exponent field).
    pub largest_exponent: u32,

    /// Smallest unbiased exponent for a normal (non-subnormal) number.
    pub smallest_normalized_exponent: i32,
}

// ============================================================================
// Predefined FloatSemantics constants
// ============================================================================

/// IEEE 754 binary32 (single precision, `f32`).
pub const IEEE_SINGLE: FloatSemantics = FloatSemantics {
    name: "ieee_single",
    bits: 32,
    exponent_bits: 8,
    mantissa_bits: 23,
    exponent_bias: 127,
    largest_exponent: 255,
    smallest_normalized_exponent: -126,
};

/// IEEE 754 binary64 (double precision, `f64`).
pub const IEEE_DOUBLE: FloatSemantics = FloatSemantics {
    name: "ieee_double",
    bits: 64,
    exponent_bits: 11,
    mantissa_bits: 52,
    exponent_bias: 1023,
    largest_exponent: 2047,
    smallest_normalized_exponent: -1022,
};

/// IEEE 754 binary128 (quad precision, `f128`).
pub const IEEE_QUAD: FloatSemantics = FloatSemantics {
    name: "ieee_quad",
    bits: 128,
    exponent_bits: 15,
    mantissa_bits: 112,
    exponent_bias: 16383,
    largest_exponent: 32767,
    smallest_normalized_exponent: -16382,
};

/// x87 80-bit double-extended precision.
/// The leading significand bit is explicit (included in the 64-bit mantissa).
/// The integer bit is always present for normal/infinity/NaN, and zero for
/// subnormals/zero.
pub const X87_DOUBLE_EXTENDED: FloatSemantics = FloatSemantics {
    name: "x87_double_extended",
    bits: 80,
    exponent_bits: 15,
    mantissa_bits: 64, // includes explicit integer bit
    exponent_bias: 16383,
    largest_exponent: 32767,
    smallest_normalized_exponent: -16382,
};

/// PowerPC double-double: a pair of IEEE double values whose sum approximates
/// quad precision. Not a true IEEE 754 format but supported for compatibility.
pub const PPC_DOUBLE_DOUBLE: FloatSemantics = FloatSemantics {
    name: "ppc_double_double",
    bits: 128,
    exponent_bits: 11,
    mantissa_bits: 106, // roughly 53+53 minus overlap
    exponent_bias: 1023,
    largest_exponent: 2047,
    smallest_normalized_exponent: -1022,
};

// ============================================================================
// APFloat
// ============================================================================

/// Arbitrary-precision IEEE 754 floating-point number.
///
/// Internally the significand is stored as a little-endian array of `u64`
/// words. It is always normalized: for finite non-zero values the most
/// significant bit of the most significant word is 1 (unless the number
/// is subnormal, in which case the MSB may be 0 and `exponent == 0`).
#[derive(Clone)]
pub struct APFloat {
    /// IEEE 754 semantics for this float.
    pub(crate) semantics: FloatSemantics,

    /// Sign bit: `true` means negative.
    pub(crate) sign: bool,

    /// Raw biased exponent as stored in the IEEE format.
    /// Special values: `0` for zero/subnormal; `largest_exponent` for inf/NaN.
    pub(crate) exponent: u32,

    /// Significand (mantissa) as arbitrary-precision integer.
    /// Stored little-endian: `significand[0]` is the least-significant word.
    /// For normal numbers the MSB is 1 (the implicit/explicit leading bit).
    /// For zero the significand is empty or all-zero.
    pub(crate) significand: Vec<u64>,

    /// Exception flags accumulated from the last operation, OR'd together.
    pub(crate) flags: u32,
}

// ============================================================================
// Internal helpers — format queries
// ============================================================================

impl FloatSemantics {
    /// The total number of significand bits including the implicit/explicit
    /// leading bit.
    #[inline]
    fn total_significand_bits(&self) -> u32 {
        if self.is_x87() {
            // x87: mantissa_bits includes the explicit leading bit.
            self.mantissa_bits
        } else {
            // IEEE: mantissa_bits is the fraction width; add 1 for implicit bit.
            self.mantissa_bits + 1
        }
    }

    /// The number of fraction bits after the binary point.
    #[inline]
    fn fraction_bits(&self) -> u32 {
        if self.is_x87() {
            self.mantissa_bits - 1
        } else {
            self.mantissa_bits
        }
    }

    /// Returns `true` if this is the x87 double-extended format.
    #[inline]
    fn is_x87(&self) -> bool {
        // x87 is uniquely 80 bits with 15 exponent bits.
        self.bits == 80 && self.exponent_bits == 15
    }

    /// Returns `true` if this is the PPC double-double format.
    #[inline]
    fn is_ppc_double_double(&self) -> bool {
        self.name == "ppc_double_double"
    }
}

// ============================================================================
// Internal helpers — APInt / significand conversions
// ============================================================================

impl APFloat {
    /// Build an APInt from the significand words, padded to `width` bits.
    #[inline]
    fn sig_apint(&self, width: u32) -> APInt {
        let mut words = self.significand.clone();
        // Ensure we don't pass empty words for zero (APInt requires bit_width≥1).
        if words.is_empty() {
            words.push(0);
        }
        // Strip trailing zero words beyond width.
        let nw = ((width as usize + 63) / 64).max(1);
        words.resize(nw, 0);
        APInt::from_words(width, words)
    }

    /// Set the significand from an APInt, normalizing the word vector.
    #[allow(dead_code)]
    fn set_sig_from_apint(&mut self, val: &APInt) {
        self.significand = val.get_raw_data().to_vec();
        // Trim trailing zero words.
        while self.significand.len() > 1 && self.significand.last() == Some(&0) {
            self.significand.pop();
        }
    }

    /// Count leading zeros in the significand (word-level, bit-level).
    fn sig_leading_zeros(&self) -> u32 {
        if self.significand.is_empty() || self.significand.iter().all(|&w| w == 0) {
            return 0; // caller should check for zero first
        }
        let top = *self.significand.last().unwrap();
        let word_lz = top.leading_zeros();
        // The leading zeros from lower, all-zero words:
        let mut extra_lz = 0u32;
        for &w in self.significand.iter().rev().skip(1) {
            if w == 0 {
                extra_lz += 64;
            } else {
                break;
            }
        }
        extra_lz + word_lz
    }

    /// Shift the significand left by `shift` bits (logical).
    fn sig_shl(&mut self, shift: u32) {
        if shift == 0 || self.significand.is_empty() {
            return;
        }
        let word_shift = (shift / 64) as usize;
        let bit_shift = (shift % 64) as u32;

        let old_len = self.significand.len();
        // New length estimate; we'll trim afterwards.
        self.significand.resize(old_len + word_shift + 1, 0);

        if bit_shift == 0 {
            for i in (0..old_len).rev() {
                self.significand[i + word_shift] = self.significand[i];
            }
            for i in 0..word_shift {
                self.significand[i] = 0;
            }
        } else {
            let inv_shift = 64 - bit_shift;
            for i in (0..old_len).rev() {
                let lo = self.significand[i] << bit_shift;
                let hi = if i > 0 {
                    self.significand[i - 1] >> inv_shift
                } else {
                    0
                };
                self.significand[i + word_shift + 1] = hi;
                self.significand[i + word_shift] = lo;
            }
            // Clear the low words that were shifted out.
            for i in 0..word_shift {
                self.significand[i] = 0;
            }
            // The word at index `word_shift` was already set; ensure
            // lower words are zero.
            if word_shift > 0 {
                self.significand[word_shift - 1] = 0;
            }
        }

        // Trim trailing zeros.
        while self.significand.len() > 1 && self.significand.last() == Some(&0) {
            self.significand.pop();
        }
    }

    /// Shift the significand right by `shift` bits (logical).
    /// Returns the sticky bit: `true` if any `1` bit was shifted out.
    fn sig_shr(&mut self, shift: u32) -> bool {
        if shift == 0 {
            return false;
        }
        if self.significand.is_empty() || self.significand.iter().all(|&w| w == 0) {
            self.significand.clear();
            self.significand.push(0);
            return false;
        }

        let word_shift = (shift / 64) as usize;
        let bit_shift = (shift % 64) as u32;
        let mut sticky = false;

        // Check for bits shifted out below word_shift.
        for &w in self.significand.iter().take(word_shift) {
            if w != 0 {
                sticky = true;
            }
        }
        // Check for bits shifted out in the word at word_shift.
        if word_shift < self.significand.len() && bit_shift > 0 {
            let mask = (1u64 << bit_shift) - 1;
            if self.significand[word_shift] & mask != 0 {
                sticky = true;
            }
        }

        let old_len = self.significand.len();
        let new_len = if old_len > word_shift {
            old_len - word_shift
        } else {
            0
        };

        if new_len == 0 {
            self.significand.clear();
            self.significand.push(0);
            return sticky;
        }

        if bit_shift == 0 {
            for i in 0..new_len {
                self.significand[i] = self.significand[i + word_shift];
            }
        } else {
            let inv_shift = 64 - bit_shift;
            for i in 0..new_len {
                let lo = self.significand[i + word_shift] >> bit_shift;
                let hi = if i + word_shift + 1 < old_len {
                    self.significand[i + word_shift + 1] << inv_shift
                } else {
                    0
                };
                self.significand[i] = lo | hi;
            }
        }
        self.significand.truncate(new_len);

        // Trim trailing zeros.
        while self.significand.len() > 1 && self.significand.last() == Some(&0) {
            self.significand.pop();
        }

        sticky
    }

    /// Get the position (0-based from LSB) of the most significant bit in
    /// the significand. Returns 0 if the significand is zero.
    #[allow(dead_code)]
    fn sig_msb_position(&self) -> u32 {
        if self.significand.is_empty() || self.significand.iter().all(|&w| w == 0) {
            return 0;
        }
        for i in (0..self.significand.len()).rev() {
            let w = self.significand[i];
            if w != 0 {
                return (i as u32 * 64) + 63 - w.leading_zeros();
            }
        }
        0
    }

    /// Get a single bit from the significand.
    fn sig_get_bit(&self, pos: u32) -> bool {
        let word = (pos / 64) as usize;
        let bit = pos % 64;
        if word < self.significand.len() {
            (self.significand[word] >> bit) & 1 != 0
        } else {
            false
        }
    }

    /// Returns `true` if the significand is exactly zero.
    fn sig_is_zero(&self) -> bool {
        self.significand.is_empty() || self.significand.iter().all(|&w| w == 0)
    }
}

// ============================================================================
// Construction
// ============================================================================

impl APFloat {
    /// Create a positive zero with the given semantics.
    pub fn new(sem: FloatSemantics) -> Self {
        APFloat {
            semantics: sem,
            sign: false,
            exponent: 0,
            significand: Vec::new(),
            flags: APFLOAT_OK,
        }
    }

    /// Create a zero with the given sign.
    pub fn get_zero(sem: FloatSemantics, sign: bool) -> Self {
        APFloat {
            semantics: sem,
            sign,
            exponent: 0,
            significand: Vec::new(),
            flags: APFLOAT_OK,
        }
    }

    /// Create an infinity with the given sign.
    pub fn get_inf(sem: FloatSemantics, sign: bool) -> Self {
        let mut sig = Vec::new();
        if sem.is_x87() {
            // x87 infinity must have the explicit integer bit set.
            sig.push(1u64 << 63);
        }
        APFloat {
            semantics: sem,
            sign,
            exponent: sem.largest_exponent,
            significand: sig,
            flags: APFLOAT_OK,
        }
    }

    /// Create a NaN with the given sign. If `payload` is provided, the low bits
    /// of the payload are set in the significand; the quiet bit (MSB of
    /// mantissa) is set to make this a quiet NaN. For signaling NaN use
    /// `get_snan`.
    pub fn get_nan(sem: FloatSemantics, sign: bool, payload: Option<u64>) -> Self {
        let qbit_pos = sem.fraction_bits() - 1; // MSB of fraction field is quiet bit
        let mut payload_val = payload.unwrap_or(0);

        if sem.is_x87() {
            // x87: explicit integer bit must be set; quiet bit is bit 62.
            let int_bit = 1u64 << 63;
            let qbit = 1u64 << 62;
            payload_val &= (1u64 << 62) - 1;
            let sig_word = int_bit | qbit | payload_val;
            APFloat {
                semantics: sem,
                sign,
                exponent: sem.largest_exponent,
                significand: vec![sig_word],
                flags: APFLOAT_OK,
            }
        } else {
            // Ensure the quiet bit is set.
            let qbit = 1u64 << qbit_pos;
            payload_val |= qbit;
            // Build significand words.
            let num_words = ((qbit_pos + 64) / 64) as usize;
            let mut words = vec![0u64; num_words.max(1)];
            let w = (qbit_pos / 64) as usize;
            if w < words.len() {
                words[w] = payload_val;
            }
            // For wider formats (e.g. quad), spread high bits of payload
            // across additional words.
            if qbit_pos >= 64 && w + 1 < words.len() {
                // payload_val is only u64, so high bits are always 0.
                // The quiet bit set above in the same word covers all cases.
            }
            // Trim leading zeros.
            while words.len() > 1 && words.last() == Some(&0) {
                words.pop();
            }
            APFloat {
                semantics: sem,
                sign,
                exponent: sem.largest_exponent,
                significand: words,
                flags: APFLOAT_OK,
            }
        }
    }

    /// Create a quiet NaN with the given sign.
    pub fn get_qnan(sem: FloatSemantics, sign: bool) -> Self {
        let mut result = Self::get_nan(sem, sign, Some(0));
        // The quiet bit is already set by get_nan. For x87, ensure we set bit 62.
        if sem.is_x87() {
            if result.significand.is_empty() {
                result.significand.push(1u64 << 63 | 1u64 << 62);
            } else {
                result.significand[0] |= 1u64 << 62;
            }
        }
        result
    }

    /// Create a signaling NaN with the given sign.
    /// A signaling NaN has the quiet bit clear and at least one other mantissa
    /// bit set.
    pub fn get_snan(sem: FloatSemantics, sign: bool) -> Self {
        let frac_bits = sem.fraction_bits();
        if sem.is_x87() {
            // x87: explicit integer bit set, quiet bit (62) clear, some lower bit set.
            let int_bit = 1u64 << 63;
            let sig_word = int_bit | 1u64; // set a low bit to make it non-zero
            APFloat {
                semantics: sem,
                sign,
                exponent: sem.largest_exponent,
                significand: vec![sig_word],
                flags: APFLOAT_OK,
            }
        } else {
            // Set a low mantissa bit, leave quiet bit (frac_bits-1) clear.
            let qbit_pos = frac_bits - 1;
            let mut payload: u64 = 1; // set bit 0 of fraction
                                      // Make sure we don't accidentally set the quiet bit.
            if qbit_pos < 64 {
                payload &= !(1u64 << qbit_pos);
            }
            let mut result = APFloat {
                semantics: sem,
                sign,
                exponent: sem.largest_exponent,
                significand: Vec::new(),
                flags: APFLOAT_OK,
            };
            // Build significand.
            let num_words = ((frac_bits + 63) / 64) as usize;
            result.significand = vec![0u64; num_words.max(1)];
            if qbit_pos < 64 {
                result.significand[0] = payload;
            } else {
                let w = (qbit_pos / 64) as usize;
                result.significand[w] = payload;
            }
            // Trim.
            while result.significand.len() > 1 && result.significand.last() == Some(&0) {
                result.significand.pop();
            }
            result
        }
    }

    /// Convert from an `f64` (auto-detects `IEEE_DOUBLE` semantics).
    pub fn from_f64(val: f64) -> Self {
        let bits: u64 = val.to_bits();
        Self::from_bits(bits, IEEE_DOUBLE)
    }

    /// Convert from an `f32` (auto-detects `IEEE_SINGLE` semantics).
    pub fn from_f32(val: f32) -> Self {
        let bits: u32 = val.to_bits();
        Self::from_bits(bits as u64, IEEE_SINGLE)
    }

    /// Reinterpret raw bits as an APFloat using the given semantics.
    fn from_bits(bits: u64, sem: FloatSemantics) -> Self {
        let sign = (bits >> (sem.bits - 1)) & 1 != 0;
        let exp_mask = (1u64 << sem.exponent_bits) - 1;
        let mantissa_mask = (1u64 << sem.mantissa_bits) - 1;
        let exp_shift = sem.mantissa_bits;

        let raw_exp = ((bits >> exp_shift) & exp_mask) as u32;
        let raw_mantissa = bits & mantissa_mask;

        let mut result = APFloat {
            semantics: sem,
            sign,
            exponent: raw_exp,
            significand: Vec::new(),
            flags: APFLOAT_OK,
        };

        if raw_exp == sem.largest_exponent {
            // Infinity or NaN.
            if raw_mantissa == 0 {
                // Infinity.
                if sem.is_x87() {
                    result.significand = vec![1u64 << 63];
                }
            } else {
                // NaN.
                let nw = ((sem.mantissa_bits as usize + 63) / 64).max(1);
                let mut words = vec![0u64; nw];
                words[0] = raw_mantissa;
                if sem.is_x87() {
                    // Ensure explicit integer bit.
                    words[0] |= 1u64 << 63;
                }
                while words.len() > 1 && words.last() == Some(&0) {
                    words.pop();
                }
                result.significand = words;
            }
        } else if raw_exp == 0 {
            // Zero or subnormal.
            if raw_mantissa != 0 {
                // Subnormal.
                let nw = ((sem.mantissa_bits as usize + 63) / 64).max(1);
                let mut words = vec![0u64; nw];
                words[0] = raw_mantissa;
                if sem.is_x87() {
                    // x87 subnormals: integer bit is 0, mantissa_bits=64 so we
                    // store the raw mantissa which already has integer bit = 0.
                }
                while words.len() > 1 && words.last() == Some(&0) {
                    words.pop();
                }
                result.significand = words;
            }
        } else {
            // Normal.
            let nw = ((sem.mantissa_bits as usize + 63) / 64).max(1);
            let mut words = vec![0u64; nw];
            words[0] = raw_mantissa;
            if !sem.is_x87() {
                // Add the implicit leading bit.
                let implicit_pos = sem.mantissa_bits;
                let w = (implicit_pos / 64) as usize;
                let b = implicit_pos % 64;
                if w < words.len() {
                    words[w] |= 1u64 << b;
                }
            }
            while words.len() > 1 && words.last() == Some(&0) {
                words.pop();
            }
            result.significand = words;
        }

        result
    }

    /// Reinterpret an APInt as an APFloat with the given semantics.
    /// The APInt's bits are treated as the raw IEEE 754 encoding.
    pub fn from_apint(val: &APInt, sem: FloatSemantics) -> Self {
        // We need to extract the bits from the APInt.
        // Build the raw bits from val's words.
        let words = val.get_raw_data();
        let mut result = APFloat {
            semantics: sem,
            sign: false,
            exponent: 0,
            significand: Vec::new(),
            flags: APFLOAT_OK,
        };

        if words.is_empty() || words.iter().all(|&w| w == 0) {
            return Self::get_zero(sem, false);
        }

        // Extract sign bit.
        let sign_bit_pos = sem.bits - 1;
        let sign_word = (sign_bit_pos / 64) as usize;
        let sign_bit = sign_bit_pos % 64;
        result.sign = sign_word < words.len() && (words[sign_word] >> sign_bit) & 1 != 0;

        // Extract exponent field.
        let exp_lo = sem.mantissa_bits;
        let exp_width = sem.exponent_bits;
        let mut raw_exp: u32 = 0;
        for i in 0..exp_width {
            let pos = exp_lo + i;
            let w = (pos / 64) as usize;
            let b = pos % 64;
            if w < words.len() && (words[w] >> b) & 1 != 0 {
                raw_exp |= 1 << i;
            }
        }
        result.exponent = raw_exp;

        // Extract mantissa.
        let mantissa_bits = sem.mantissa_bits;
        let nw = ((mantissa_bits as usize + 63) / 64).max(1);
        let mut sig_words = vec![0u64; nw];
        for i in 0..mantissa_bits {
            let w = (i / 64) as usize;
            let b = i % 64;
            if w < words.len() && (words[w] >> b) & 1 != 0 {
                let sw = (i / 64) as usize;
                let sb = i % 64;
                if sw < sig_words.len() {
                    sig_words[sw] |= 1u64 << sb;
                }
            }
        }

        if raw_exp == sem.largest_exponent {
            // Inf or NaN.
            if sig_words.iter().all(|&w| w == 0) && !sem.is_x87() {
                // Infinity.
            } else {
                if sem.is_x87() {
                    // Ensure explicit integer bit.
                    sig_words[0] |= 1u64 << 63;
                }
                while sig_words.len() > 1 && sig_words.last() == Some(&0) {
                    sig_words.pop();
                }
                result.significand = sig_words;
            }
        } else if raw_exp == 0 {
            if !sig_words.iter().all(|&w| w == 0) {
                // Subnormal.
                while sig_words.len() > 1 && sig_words.last() == Some(&0) {
                    sig_words.pop();
                }
                result.significand = sig_words;
            }
        } else {
            // Normal: add implicit leading bit for IEEE formats.
            if !sem.is_x87() {
                let implicit_pos = mantissa_bits;
                let sw = (implicit_pos / 64) as usize;
                let sb = implicit_pos % 64;
                if sw < sig_words.len() {
                    sig_words[sw] |= 1u64 << sb;
                } else {
                    sig_words.resize(sw + 1, 0);
                    sig_words[sw] = 1u64 << sb;
                }
            }
            while sig_words.len() > 1 && sig_words.last() == Some(&0) {
                sig_words.pop();
            }
            result.significand = sig_words;
        }

        result
    }

    /// The largest finite value with the given sign.
    pub fn largest(sem: FloatSemantics, sign: bool) -> Self {
        let mut result = APFloat {
            semantics: sem,
            sign,
            exponent: sem.largest_exponent - 1,
            significand: Vec::new(),
            flags: APFLOAT_OK,
        };
        // Build the maximum significand (all fraction bits = 1).
        let frac_bits = sem.fraction_bits();
        let nw = ((frac_bits as usize + 63) / 64).max(1);
        result.significand = vec![u64::MAX; nw];
        // Mask the top word to only fraction_bits.
        let rem = frac_bits % 64;
        if rem != 0 {
            let last = result.significand.last_mut().unwrap();
            *last &= (1u64 << rem) - 1;
        }
        // For IEEE, add implicit leading bit.
        if !sem.is_x87() {
            let implicit_pos = frac_bits;
            let sw = (implicit_pos / 64) as usize;
            let sb = implicit_pos % 64;
            if sw < result.significand.len() {
                result.significand[sw] |= 1u64 << sb;
            } else {
                result.significand.resize(sw + 1, 0);
                result.significand[sw] = 1u64 << sb;
            }
        } else {
            // x87: explicit leading bit is already part of mantissa_bits, which
            // is 64. The fraction is 63 bits. Set bit 63 explicitly.
            let sig_len = result.significand.len();
            if sig_len == 1 {
                let last = result.significand.last_mut().unwrap();
                *last = (1u64 << 63) | ((*last) & ((1u64 << 63) - 1));
            } else {
                // Shouldn't happen for 80-bit format.
                result.significand[0] |= 1u64 << 63;
            }
        }
        // Trim leading zero words but don't remove the implicit bit.
        while result.significand.len() > 1 && result.significand.last() == Some(&0) {
            result.significand.pop();
        }
        result
    }

    /// The smallest positive normal value with the given sign.
    pub fn smallest(sem: FloatSemantics, sign: bool) -> Self {
        Self::smallest_normalized(sem, sign)
    }

    /// The smallest positive normal (non-subnormal) value with the given sign.
    pub fn smallest_normalized(sem: FloatSemantics, sign: bool) -> Self {
        let mut result = APFloat {
            semantics: sem,
            sign,
            exponent: 1, // smallest normal biased exponent
            significand: Vec::new(),
            flags: APFLOAT_OK,
        };
        // Significand: only the implicit/explicit leading bit set.
        if sem.is_x87() {
            result.significand = vec![1u64 << 63];
        } else {
            let implicit_pos = sem.mantissa_bits;
            let sw = (implicit_pos / 64) as usize;
            let sb = implicit_pos % 64;
            result.significand = vec![0u64; sw + 1];
            result.significand[sw] = 1u64 << sb;
        }
        result
    }
}

// ============================================================================
// Classification
// ============================================================================

impl APFloat {
    /// Returns `true` if this value is ±0.
    pub fn is_zero(&self) -> bool {
        self.exponent == 0 && self.sig_is_zero()
    }

    /// Returns `true` if the sign bit is set (includes -0 and -NaN).
    pub fn is_negative(&self) -> bool {
        self.sign
    }

    /// Returns `true` if this value is ±∞.
    pub fn is_infinity(&self) -> bool {
        if self.exponent != self.semantics.largest_exponent {
            return false;
        }
        if self.semantics.is_x87() {
            // x87: infinity has explicit integer bit set and no other mantissa bits.
            if self.significand.len() == 1 {
                let w = self.significand[0];
                return w == (1u64 << 63);
            }
            return false;
        }
        self.sig_is_zero()
    }

    /// Returns `true` if this is positive infinity.
    pub fn is_positively_infinity(&self) -> bool {
        !self.sign && self.is_infinity()
    }

    /// Returns `true` if this is negative infinity.
    pub fn is_negatively_infinity(&self) -> bool {
        self.sign && self.is_infinity()
    }

    /// Returns `true` if this value is NaN (quiet or signaling).
    pub fn is_nan(&self) -> bool {
        self.exponent == self.semantics.largest_exponent && !self.sig_is_zero()
    }

    /// Returns `true` if this is a signaling NaN.
    ///
    /// A signaling NaN has the quiet bit (MSB of the fraction field) clear.
    pub fn is_signaling_nan(&self) -> bool {
        if !self.is_nan() {
            return false;
        }
        let qbit_pos = if self.semantics.is_x87() {
            62 // bit 62 is the quiet bit for x87
        } else {
            self.semantics.fraction_bits() - 1
        };
        !self.sig_get_bit(qbit_pos)
    }

    /// Returns `true` if this is a normal (non-zero, non-subnormal, non-inf,
    /// non-NaN) finite value.
    pub fn is_normal(&self) -> bool {
        if self.is_zero() || self.is_infinity() || self.is_nan() {
            return false;
        }
        self.exponent > 0
        // For x87, subnormals also have exponent=0 so the check above works.
        // For IEEE, exponent=0 means zero or subnormal.
    }

    /// Returns `true` if this is a subnormal (denormalized) value.
    pub fn is_subnormal(&self) -> bool {
        self.exponent == 0 && !self.sig_is_zero()
    }

    /// Returns `true` if this is a finite non-zero value (normal or subnormal).
    pub fn is_non_zero(&self) -> bool {
        !self.is_zero()
    }

    /// Returns `true` if this is a finite value (not NaN or infinity).
    pub fn is_finite(&self) -> bool {
        !self.is_infinity() && !self.is_nan()
    }

    /// Classify this value into an IEEE 754 category.
    pub fn classify(&self) -> FloatCategory {
        if self.is_zero() {
            FloatCategory::Zero
        } else if self.is_infinity() {
            FloatCategory::Infinity
        } else if self.is_nan() {
            FloatCategory::NaN
        } else if self.exponent == 0 {
            FloatCategory::Subnormal
        } else {
            FloatCategory::Normal
        }
    }
}

// ============================================================================
// Internal: sign/magnitude helpers
// ============================================================================

impl APFloat {
    /// Copy the sign from `other`, keeping magnitude.
    pub fn copy_sign(&self, other: &APFloat) -> APFloat {
        let mut result = self.clone();
        result.sign = other.sign;
        result.flags = APFLOAT_OK;
        result
    }

    /// Return the absolute value.
    pub fn abs(&self) -> APFloat {
        let mut result = self.clone();
        result.sign = false;
        result.flags = APFLOAT_OK;
        result
    }

    /// Return the negated value.
    pub fn negate(&self) -> APFloat {
        let mut result = self.clone();
        if !result.is_nan() {
            result.sign = !result.sign;
        }
        result.flags = APFLOAT_OK;
        result
    }
}

// ============================================================================
// Internal: normalization
// ============================================================================

impl APFloat {
    /// Normalize the significand so its MSB is 1, adjusting the exponent.
    /// If the significand is zero, sets exponent to 0.
    #[allow(dead_code)]
    /// Returns the new biased exponent (or 0 for zero).
    fn normalize(&mut self, unbiased_exp: i64) -> u32 {
        if self.sig_is_zero() {
            self.exponent = 0;
            self.significand.clear();
            return 0;
        }

        let lz = self.sig_leading_zeros();
        if lz > 0 {
            self.sig_shl(lz);
        }

        let new_unbiased = unbiased_exp - lz as i64;
        let bias = self.semantics.exponent_bias as i64;

        // Check for underflow to subnormal range.
        let min_normal_unbiased = self.semantics.smallest_normalized_exponent as i64;
        if new_unbiased < min_normal_unbiased {
            // Need to shift right to create a subnormal.
            let _shift = (min_normal_unbiased - new_unbiased) as u32;
            // We actually need to shift so that the biased exponent becomes 0.
            // The biased exponent for the smallest normal is 1.
            // If new_unbiased + bias < 1, we're in subnormal territory.
            let biased = new_unbiased + bias;
            if biased < 1 {
                // Shift right by (1 - biased) to get biased exponent = 0.
                let subnormal_shift = (1 - biased) as u32;
                let _sticky = self.sig_shr(subnormal_shift);
                // After shifting, the significand may have lost its MSB.
                // That's correct for subnormals.
                return 0;
            }
        }

        let biased = (new_unbiased + bias) as u32;
        if biased >= self.semantics.largest_exponent {
            // Overflow to infinity.
            return self.semantics.largest_exponent;
        }

        biased
    }

    /// Get the unbiased exponent.
    fn unbiased_exponent(&self) -> i64 {
        if self.exponent == 0 {
            // Zero or subnormal: effective unbiased exponent = 1 - bias.
            if self.sig_is_zero() {
                return 0;
            }
            1 - self.semantics.exponent_bias as i64
        } else {
            self.exponent as i64 - self.semantics.exponent_bias as i64
        }
    }
}

// ============================================================================
// Internal: rounding
// ============================================================================

impl APFloat {
    /// Round the significand to `target_bits` bits (the total significand width
    /// including implicit/explicit leading bit). Uses guard, round, and sticky
    /// bits positioned just below the target.
    ///
    /// Returns `(carry, inexact)` where `carry` means the rounding caused an
    /// additional carry out past the MSB (requiring exponent adjustment).
    fn round_significand(
        sig: &mut Vec<u64>,
        target_bits: u32,
        rm: RoundingMode,
        sign: bool,
    ) -> (bool, bool) {
        // target_bits is the number of bits to keep.
        // The bit at position (target_bits - 1) is the LSB we keep.
        // The guard bit is at position (target_bits - 1) in the shifted-out portion...
        // Actually, let's define:
        //   LSB   = bit position (target_bits - 1) from the bottom (0-based)
        //   Guard = bit position (target_bits - 1) — wait, that's confusing.
        //
        // Let me rephrase. We have a significand that's been shifted so that
        // the MSB is at some position. We want to keep `target_bits` bits.
        // So we look at the bits from MSB down to (MSB - target_bits + 1) as
        // the kept bits. The next three bits are:
        //   Guard: first bit beyond target_bits (at position MSB - target_bits)
        //   Round: second bit beyond target_bits (at position MSB - target_bits - 1)
        //   Sticky: OR of all remaining lower bits.
        //
        // But our significand is stored little-endian with LSB at index 0.
        // The MSB position is `msb_pos = sig_msb_position(sig)`.
        // The LSB of the kept portion is at position `msb_pos - target_bits + 1`.
        // We need to shift right by `(sig_width - target_bits)` if sig_width > target_bits.

        let msb = sig_msb_position(sig);
        // msb is 0-indexed, so number of bits = msb + 1.
        // Only round if significand has more bits than target_bits.
        if msb + 1 <= target_bits {
            return (false, false);
        }

        let shift = msb + 1 - target_bits; // number of bits to shift right

        // We need to extract guard, round, sticky before shifting.
        let guard_bit = if shift >= 1 && shift <= msb + 1 {
            get_bit_vec(sig, shift - 1)
        } else {
            false
        };
        let round_bit = if shift >= 2 && shift <= msb + 1 {
            get_bit_vec(sig, shift - 2)
        } else {
            false
        };
        let sticky = if shift > 2 {
            any_bit_below_vec(sig, shift - 2)
        } else {
            false
        };

        // Now shift right.
        let _shift_sticky = shr_vec(sig, shift);
        // After shifting, the LSB of the result is what we keep.

        // Determine whether to round up.
        let round_up = match rm {
            RoundingMode::NearestTiesToEven => {
                guard_bit && (round_bit || sticky || get_bit_vec(sig, 0))
                // LSB after shift is at bit 0. If LSB is 1, we round up on tie (to even).
            }
            RoundingMode::NearestTiesToAway => {
                // For ties-to-away, round up if guard is 1.
                guard_bit
            }
            RoundingMode::TowardZero => false,
            RoundingMode::TowardNegative => sign && (guard_bit || round_bit || sticky),
            RoundingMode::TowardPositive => !sign && (guard_bit || round_bit || sticky),
        };

        let inexact = guard_bit || round_bit || sticky;

        if round_up {
            let carry = add_one_vec(sig);
            (carry, inexact)
        } else {
            (false, inexact)
        }
    }

    /// Round the significand and update the exponent if carry occurs.
    /// This is the high-level rounding entry point.
    fn round_and_renormalize(
        sig: &mut Vec<u64>,
        biased_exp: u32,
        target_bits: u32,
        rm: RoundingMode,
        sign: bool,
        sem: &FloatSemantics,
    ) -> (u32, bool) {
        let (carry, inexact) = Self::round_significand(sig, target_bits, rm, sign);

        let mut new_exp = biased_exp;

        if carry {
            // Shift right by 1 to make room for the carry.
            let _ = shr_vec(sig, 1);
            // Set the MSB.
            let msb_pos = target_bits - 1;
            let w = (msb_pos / 64) as usize;
            let b = msb_pos % 64;
            if w >= sig.len() {
                sig.resize(w + 1, 0);
            }
            sig[w] |= 1u64 << b;

            new_exp += 1;
            if new_exp >= sem.largest_exponent {
                // Overflow to infinity.
                sig.clear();
                if sem.is_x87() {
                    sig.push(1u64 << 63);
                }
                return (sem.largest_exponent, inexact);
            }
        }

        // After rounding, the significand should have exactly target_bits bits
        // (MSB at target_bits-1). If it became zero due to rounding, set to zero.
        if sig.iter().all(|&w| w == 0) {
            sig.clear();
            return (0, inexact);
        }

        // Ensure MSB is at target_bits - 1 (normalize).
        let msb = sig_msb_position(sig);
        if msb > target_bits - 1 {
            // Shouldn't happen after rounding; but shift right if needed.
            let shift = msb - (target_bits - 1);
            let _ = shr_vec(sig, shift);
            new_exp += shift;
            if new_exp >= sem.largest_exponent {
                sig.clear();
                if sem.is_x87() {
                    sig.push(1u64 << 63);
                }
                return (sem.largest_exponent, true);
            }
        } else if msb < target_bits - 1 && new_exp > 0 {
            // The rounding may have left some leading zeros; shift left.
            let shift = (target_bits - 1) - msb;
            shl_vec(sig, shift);
            if new_exp > shift as u32 {
                new_exp -= shift as u32;
            } else {
                // Underflow to subnormal.
                let remaining = shift as u32 - new_exp;
                new_exp = 0;
                let _ = shr_vec(sig, remaining);
            }
        } else if new_exp == 0 && msb < target_bits - 1 {
            // Already subnormal; leave as-is.
        }

        (new_exp, inexact)
    }
}

// --------------------------------------------------------------------------
// Low-level Vec<u64> bit operations (used by rounding)
// --------------------------------------------------------------------------

fn sig_msb_position(words: &[u64]) -> u32 {
    if words.is_empty() {
        return 0;
    }
    for i in (0..words.len()).rev() {
        let w = words[i];
        if w != 0 {
            return (i as u32 * 64) + 63 - w.leading_zeros();
        }
    }
    0
}

fn get_bit_vec(words: &[u64], pos: u32) -> bool {
    let w = (pos / 64) as usize;
    let b = pos % 64;
    w < words.len() && (words[w] >> b) & 1 != 0
}

fn any_bit_below_vec(words: &[u64], pos: u32) -> bool {
    if pos == 0 {
        return false;
    }
    let full_words = (pos / 64) as usize;
    for i in 0..full_words.min(words.len()) {
        if words[i] != 0 {
            return true;
        }
    }
    let rem = pos % 64;
    if rem > 0 {
        let w = full_words;
        if w < words.len() {
            let mask = (1u64 << rem) - 1;
            if words[w] & mask != 0 {
                return true;
            }
        }
    }
    false
}

fn shr_vec(words: &mut Vec<u64>, shift: u32) -> bool {
    if shift == 0 {
        return false;
    }
    if words.is_empty() || words.iter().all(|&w| w == 0) {
        words.clear();
        words.push(0);
        return false;
    }

    let word_shift = (shift / 64) as usize;
    let bit_shift = (shift % 64) as u32;
    let mut sticky = false;

    for &w in words.iter().take(word_shift) {
        if w != 0 {
            sticky = true;
        }
    }
    if word_shift < words.len() && bit_shift > 0 {
        let mask = (1u64 << bit_shift) - 1;
        if words[word_shift] & mask != 0 {
            sticky = true;
        }
    }

    let old_len = words.len();
    let new_len = if old_len > word_shift {
        old_len - word_shift
    } else {
        0
    };

    if new_len == 0 {
        words.clear();
        words.push(0);
        return sticky;
    }

    if bit_shift == 0 {
        for i in 0..new_len {
            words[i] = words[i + word_shift];
        }
    } else {
        let inv_shift = 64 - bit_shift;
        for i in 0..new_len {
            let lo = words[i + word_shift] >> bit_shift;
            let hi = if i + word_shift + 1 < old_len {
                words[i + word_shift + 1] << inv_shift
            } else {
                0
            };
            words[i] = lo | hi;
        }
    }
    words.truncate(new_len);

    while words.len() > 1 && words.last() == Some(&0) {
        words.pop();
    }

    sticky
}

fn shl_vec(words: &mut Vec<u64>, shift: u32) {
    if shift == 0 || words.is_empty() {
        return;
    }
    let word_shift = (shift / 64) as usize;
    let bit_shift = (shift % 64) as u32;
    let old_len = words.len();
    words.resize(old_len + word_shift + 1, 0);

    if bit_shift == 0 {
        for i in (0..old_len).rev() {
            words[i + word_shift] = words[i];
        }
        for i in 0..word_shift {
            words[i] = 0;
        }
    } else {
        let inv_shift = 64 - bit_shift;
        for i in (0..old_len).rev() {
            let lo = words[i] << bit_shift;
            let hi = if i > 0 { words[i - 1] >> inv_shift } else { 0 };
            words[i + word_shift + 1] = hi;
            words[i + word_shift] = lo;
        }
        for i in 0..word_shift {
            words[i] = 0;
        }
    }

    while words.len() > 1 && words.last() == Some(&0) {
        words.pop();
    }
}

/// Add 1 to the little-endian word vector. Returns `true` on carry out.
fn add_one_vec(words: &mut Vec<u64>) -> bool {
    for w in words.iter_mut() {
        let (sum, carry) = w.overflowing_add(1);
        *w = sum;
        if !carry {
            return false;
        }
    }
    true
}

// ============================================================================
// Internal: NaN propagation
// ============================================================================

impl APFloat {
    /// If either operand is NaN, return the propagated NaN and set
    /// `APFLOAT_INVALID_OP` if a signaling NaN is involved.
    fn propagate_nan(&self, other: &APFloat) -> Option<APFloat> {
        let self_nan = self.is_nan();
        let other_nan = other.is_nan();

        if !self_nan && !other_nan {
            return None;
        }

        let mut result = if other_nan {
            other.clone()
        } else {
            self.clone()
        };

        // If either is signaling, set invalid op.
        if self.is_signaling_nan() || other.is_signaling_nan() {
            result.flags |= APFLOAT_INVALID_OP;
        }

        // Make the result a quiet NaN.
        let qbit_pos = if result.semantics.is_x87() {
            62
        } else {
            result.semantics.fraction_bits() - 1
        };
        let w = (qbit_pos / 64) as usize;
        let b = qbit_pos % 64;
        if w < result.significand.len() {
            result.significand[w] |= 1u64 << b;
        } else if !result.significand.is_empty() {
            // Can't set bit beyond words; extend.
            result.significand.resize(w + 1, 0);
            result.significand[w] = 1u64 << b;
        } else {
            // Empty significand for NaN shouldn't happen.
            result = Self::get_qnan(result.semantics, result.sign);
        }

        result.exponent = result.semantics.largest_exponent;
        Some(result)
    }
}

// ============================================================================
// Arithmetic operations
// ============================================================================

impl APFloat {
    /// Addition: `self + other` with the given rounding mode.
    pub fn add(&self, other: &APFloat, rm: RoundingMode) -> APFloat {
        if let Some(nan) = self.propagate_nan(other) {
            return nan;
        }

        // Handle infinities and zeros.
        let self_inf = self.is_infinity();
        let other_inf = other.is_infinity();
        let self_zero = self.is_zero();
        let other_zero = other.is_zero();

        if self_inf && other_inf {
            if self.sign == other.sign {
                // inf + inf = inf, -inf + -inf = -inf
                let mut result = self.clone();
                result.flags = APFLOAT_OK;
                return result;
            } else {
                // inf + -inf = NaN
                let mut result = Self::get_qnan(self.semantics, false);
                result.flags = APFLOAT_INVALID_OP;
                return result;
            }
        }
        if self_inf {
            let mut result = self.clone();
            result.flags = APFLOAT_OK;
            return result;
        }
        if other_inf {
            let mut result = other.clone();
            result.flags = APFLOAT_OK;
            return result;
        }
        if self_zero && other_zero {
            // -0 + -0 = -0, -0 + +0 = +0 (IEEE: sign of sum when both zero)
            let result_sign = if self.sign && other.sign {
                true
            } else if rm == RoundingMode::TowardNegative {
                true
            } else {
                false
            };
            let mut result = Self::get_zero(self.semantics, result_sign);
            result.flags = APFLOAT_OK;
            return result;
        }
        if self_zero {
            let mut result = other.clone();
            result.flags = APFLOAT_OK;
            return result;
        }
        if other_zero {
            let mut result = self.clone();
            result.flags = APFLOAT_OK;
            return result;
        }

        // Both operands are finite non-zero.
        self.add_sub_internal(other, false, rm)
    }

    /// Subtraction: `self - other` with the given rounding mode.
    pub fn sub(&self, other: &APFloat, rm: RoundingMode) -> APFloat {
        let mut neg_other = other.clone();
        neg_other.sign = !neg_other.sign;
        self.add(&neg_other, rm)
    }

    /// Multiplication: `self * other` with the given rounding mode.
    pub fn mul(&self, other: &APFloat, rm: RoundingMode) -> APFloat {
        if let Some(nan) = self.propagate_nan(other) {
            return nan;
        }

        let self_inf = self.is_infinity();
        let other_inf = other.is_infinity();
        let self_zero = self.is_zero();
        let other_zero = other.is_zero();

        // inf * 0 = NaN
        if (self_inf && other_zero) || (self_zero && other_inf) {
            let mut result = Self::get_qnan(self.semantics, false);
            result.flags = APFLOAT_INVALID_OP;
            return result;
        }

        let result_sign = self.sign ^ other.sign;

        if self_inf || other_inf {
            let mut result = Self::get_inf(self.semantics, result_sign);
            result.flags = APFLOAT_OK;
            return result;
        }
        if self_zero || other_zero {
            let mut result = Self::get_zero(self.semantics, result_sign);
            result.flags = APFLOAT_OK;
            return result;
        }

        // Finite non-zero multiplication.
        let target_bits = self.semantics.total_significand_bits();
        let internal_width = target_bits * 4 + 10; // plenty of room

        let a = self.sig_apint(internal_width);
        let b = other.sig_apint(internal_width);

        // Full product.
        let product = a.mul(&b);

        // Compute the unbiased exponent sum.
        // Since each significand is scaled by 2^m (m = mantissa_bits), the product
        // is scaled by 2^(2m). The result should be scaled by 2^m, so we subtract m.
        let exp_a = self.unbiased_exponent();
        let exp_b = other.unbiased_exponent();
        let mut unbiased_exp = exp_a + exp_b - self.semantics.mantissa_bits as i64;

        let mut sig_words = product.get_raw_data().to_vec();
        while sig_words.len() > 1 && sig_words.last() == Some(&0) {
            sig_words.pop();
        }

        // Normalize and round.
        let (biased_exp, inexact) = Self::normalize_and_round(
            &mut sig_words,
            &mut unbiased_exp,
            target_bits,
            rm,
            result_sign,
            &self.semantics,
        );

        let mut result = APFloat {
            semantics: self.semantics,
            sign: result_sign,
            exponent: biased_exp,
            significand: sig_words,
            flags: if inexact { APFLOAT_INEXACT } else { APFLOAT_OK },
        };

        // Handle overflow/underflow flags.
        if biased_exp == self.semantics.largest_exponent && !result.is_nan() {
            result.flags |= APFLOAT_OVERFLOW | APFLOAT_INEXACT;
        }
        if biased_exp == 0 && !result.is_zero() && inexact {
            result.flags |= APFLOAT_UNDERFLOW;
        }

        // Trim significand.
        while result.significand.len() > 1 && result.significand.last() == Some(&0) {
            result.significand.pop();
        }

        result
    }

    /// Division: `self / other` with the given rounding mode.
    pub fn div(&self, other: &APFloat, rm: RoundingMode) -> APFloat {
        if let Some(nan) = self.propagate_nan(other) {
            return nan;
        }

        let self_inf = self.is_infinity();
        let other_inf = other.is_infinity();
        let self_zero = self.is_zero();
        let other_zero = other.is_zero();

        // 0/0 or inf/inf = NaN
        if (self_zero && other_zero) || (self_inf && other_inf) {
            let mut result = Self::get_qnan(self.semantics, false);
            result.flags = APFLOAT_INVALID_OP;
            return result;
        }

        let result_sign = self.sign ^ other.sign;

        // x/0 = inf (unless x=0, handled above)
        if other_zero {
            let mut result = Self::get_inf(self.semantics, result_sign);
            result.flags = APFLOAT_DIV_BY_ZERO;
            return result;
        }
        // x/inf = 0
        if other_inf {
            let mut result = Self::get_zero(self.semantics, result_sign);
            result.flags = APFLOAT_OK;
            return result;
        }
        // 0/x = 0
        if self_zero {
            let mut result = Self::get_zero(self.semantics, result_sign);
            result.flags = APFLOAT_OK;
            return result;
        }
        // inf/x = inf
        if self_inf {
            let mut result = Self::get_inf(self.semantics, result_sign);
            result.flags = APFLOAT_OK;
            return result;
        }

        // Finite non-zero division.
        let target_bits = self.semantics.total_significand_bits();
        let internal_width = target_bits * 4 + 10;

        // We need quotient with enough precision: shift dividend left by
        // (target_bits + 2) to get guard/round/sticky bits in the quotient.
        let extra_bits = target_bits + 2;
        let div_width = internal_width + extra_bits;

        let a = self.sig_apint(div_width);
        let b = other.sig_apint(div_width);

        // Shift dividend left to get enough precision in quotient.
        let shifted_a = a.shl(extra_bits);

        // Divide.
        let quotient = shifted_a.udiv(&b);

        let exp_a = self.unbiased_exponent();
        let exp_b = other.unbiased_exponent();
        // The quotient (sig_a / sig_b) has a scaling of 1 (the sig_a*2^k shift
        // cancels). But the result needs scaling by 2^m, so add m back.
        let mut unbiased_exp = exp_a - exp_b - extra_bits as i64
            + self.semantics.mantissa_bits as i64;

        let mut sig_words = quotient.get_raw_data().to_vec();
        while sig_words.len() > 1 && sig_words.last() == Some(&0) {
            sig_words.pop();
        }

        let (biased_exp, inexact) = Self::normalize_and_round(
            &mut sig_words,
            &mut unbiased_exp,
            target_bits,
            rm,
            result_sign,
            &self.semantics,
        );

        let mut result = APFloat {
            semantics: self.semantics,
            sign: result_sign,
            exponent: biased_exp,
            significand: sig_words,
            flags: if inexact { APFLOAT_INEXACT } else { APFLOAT_OK },
        };

        if biased_exp == self.semantics.largest_exponent && !result.is_nan() {
            result.flags |= APFLOAT_OVERFLOW | APFLOAT_INEXACT;
        }
        if biased_exp == 0 && !result.is_zero() && inexact {
            result.flags |= APFLOAT_UNDERFLOW;
        }

        while result.significand.len() > 1 && result.significand.last() == Some(&0) {
            result.significand.pop();
        }

        result
    }

    /// Remainder: IEEE 754 remainder operation (not modulo).
    /// Computes `self - n * other` where `n` is the integer nearest to
    /// `self / other` (ties to even).
    pub fn rem(&self, other: &APFloat, rm: RoundingMode) -> APFloat {
        if let Some(nan) = self.propagate_nan(other) {
            return nan;
        }

        let self_inf = self.is_infinity();
        let other_inf = other.is_infinity();
        let self_zero = self.is_zero();
        let other_zero = other.is_zero();

        if self_inf || other_zero {
            let mut result = Self::get_qnan(self.semantics, false);
            result.flags = APFLOAT_INVALID_OP;
            return result;
        }
        if other_inf {
            let mut result = self.clone();
            result.flags = APFLOAT_OK;
            return result;
        }
        if self_zero {
            let mut result = self.clone();
            result.flags = APFLOAT_OK;
            return result;
        }

        // Compute n = round_to_nearest_integer(self / other)
        let q = self.div(other, RoundingMode::NearestTiesToEven);
        let n = q.round_to_integral(RoundingMode::NearestTiesToEven);
        let prod = n.mul(other, rm);
        self.sub(&prod, rm)
    }

    /// Fused multiply-add: `self * mul + add` with a single rounding.
    pub fn fma(&self, mul: &APFloat, add: &APFloat, rm: RoundingMode) -> APFloat {
        // Handle NaN propagation.
        if let Some(nan) = self.propagate_nan(mul) {
            let nan2 = nan.propagate_nan(add);
            return nan2.unwrap_or(nan);
        }
        if let Some(nan) = self.propagate_nan(add) {
            return nan;
        }

        let self_inf = self.is_infinity();
        let mul_inf = mul.is_infinity();
        let add_inf = add.is_infinity();
        let self_zero = self.is_zero();
        let mul_zero = mul.is_zero();
        let add_zero = add.is_zero();

        // inf * 0 = NaN
        if (self_inf && mul_zero) || (self_zero && mul_inf) {
            let mut result = Self::get_qnan(self.semantics, false);
            result.flags = APFLOAT_INVALID_OP;
            return result;
        }

        // Product is inf or zero, and add is the opposite inf: inf + -inf = NaN
        let prod_sign = self.sign ^ mul.sign;
        let prod_inf = self_inf || mul_inf;
        let prod_zero = self_zero || mul_zero;

        if prod_inf && add_inf && prod_sign != add.sign {
            let mut result = Self::get_qnan(self.semantics, false);
            result.flags = APFLOAT_INVALID_OP;
            return result;
        }

        // If product is inf, result is inf.
        if prod_inf {
            let mut result = Self::get_inf(self.semantics, prod_sign);
            result.flags = APFLOAT_OK;
            return result;
        }
        // If add is inf, result is add.
        if add_inf {
            let mut result = add.clone();
            result.flags = APFLOAT_OK;
            return result;
        }
        // If product is zero, result is add.
        if prod_zero {
            let mut result = add.clone();
            result.flags = APFLOAT_OK;
            return result;
        }
        // If add is zero, result is product.
        if add_zero {
            let mut result = self.mul(mul, rm);
            result.flags = APFLOAT_OK;
            return result;
        }

        // Full precision FMA: compute exact product, then add.
        let target_bits = self.semantics.total_significand_bits();
        let wide_width = target_bits * 4 + 20;

        let a = self.sig_apint(wide_width);
        let b = mul.sig_apint(wide_width);

        let product = a.mul(&b);
        let prod_ap = product; // use product directly as APInt

        let prod_exp_unbiased = self.unbiased_exponent() + mul.unbiased_exponent() - self.semantics.mantissa_bits as i64;

        // Align addend with product.
        let add_exp_unbiased = add.unbiased_exponent();
        let exp_diff = prod_exp_unbiased - add_exp_unbiased;

        let add_ap = add.sig_apint(wide_width);

        let mut sum_sig: Vec<u64>;
        let result_sign: bool;

        // Align by shifting the smaller APInt right.
        let (aligned_prod, aligned_add, larger_exp) = if exp_diff >= 0 {
            // Product exponent is larger or equal; shift addend right.
            let add_raw = add_ap.get_raw_data().to_vec();
            let mut shifted_add = add_raw.to_vec();
            let _ = shr_vec(&mut shifted_add, exp_diff as u32);
            let ap_shifted = APInt::from_words(wide_width, shifted_add);
            (prod_ap, ap_shifted, prod_exp_unbiased)
        } else {
            // Addend exponent is larger; shift product right.
            let prod_raw = prod_ap.get_raw_data().to_vec();
            let mut shifted_prod = prod_raw.to_vec();
            let _ = shr_vec(&mut shifted_prod, (-exp_diff) as u32);
            let ap_shifted = APInt::from_words(wide_width, shifted_prod);
            (ap_shifted, add_ap, add_exp_unbiased)
        };

        if prod_sign == add.sign {
            // Same sign: add magnitudes.
            sum_sig = aligned_prod.add(&aligned_add).get_raw_data().to_vec();
            result_sign = prod_sign;
        } else {
            // Different signs: subtract smaller from larger magnitude.
            let prod_ge_add = aligned_prod.uge(&aligned_add);
            if prod_ge_add {
                sum_sig = aligned_prod.sub(&aligned_add).get_raw_data().to_vec();
                result_sign = prod_sign;
            } else {
                sum_sig = aligned_add.sub(&aligned_prod).get_raw_data().to_vec();
                result_sign = add.sign;
            }
        }

        while sum_sig.len() > 1 && sum_sig.last() == Some(&0) {
            sum_sig.pop();
        }

        if sum_sig.iter().all(|&w| w == 0) {
            let mut result = Self::get_zero(self.semantics, false);
            result.flags = APFLOAT_OK;
            return result;
        }

        let mut unbiased_exp = larger_exp;

        let (biased_exp, inexact) = Self::normalize_and_round(
            &mut sum_sig,
            &mut unbiased_exp,
            target_bits,
            rm,
            result_sign,
            &self.semantics,
        );

        let mut result = APFloat {
            semantics: self.semantics,
            sign: result_sign,
            exponent: biased_exp,
            significand: sum_sig,
            flags: if inexact { APFLOAT_INEXACT } else { APFLOAT_OK },
        };

        if biased_exp == self.semantics.largest_exponent && !result.is_nan() {
            result.flags |= APFLOAT_OVERFLOW | APFLOAT_INEXACT;
        }
        if biased_exp == 0 && !result.is_zero() && inexact {
            result.flags |= APFLOAT_UNDERFLOW;
        }

        while result.significand.len() > 1 && result.significand.last() == Some(&0) {
            result.significand.pop();
        }

        result
    }

    /// Square root.
    pub fn sqrt(&self, rm: RoundingMode) -> APFloat {
        if self.is_nan() {
            return self.clone();
        }
        if self.is_zero() {
            let mut result = self.clone();
            result.flags = APFLOAT_OK;
            return result;
        }
        if self.is_infinity() && !self.sign {
            let mut result = self.clone();
            result.flags = APFLOAT_OK;
            return result;
        }
        if self.sign {
            // sqrt(negative) = NaN
            let mut result = Self::get_qnan(self.semantics, false);
            result.flags = APFLOAT_INVALID_OP;
            return result;
        }

        // Newton-Raphson: x_{n+1} = (x_n + S/x_n) / 2
        let target_bits = self.semantics.total_significand_bits();
        let internal_width = target_bits * 4 + 10;

        // To get exact rounding, compute sqrt(sig << (2*target_bits)).
        // This gives the result with target_bits extra fractional bits.
        let extra_shift = 2 * target_bits;

        let s = self.sig_apint(internal_width);
        let s_exp = self.unbiased_exponent();

        // Adjust exponent to be even.
        let (mut adjusted_s, adjusted_exp) = if s_exp % 2 != 0 {
            (s.shl(1), s_exp - 1)
        } else {
            (s, s_exp)
        };

        // Shift for extra precision: compute sqrt(sig << extra_shift) = sqrt(sig) << (extra_shift/2).
        adjusted_s = adjusted_s.shl(extra_shift);

        // Initial guess: 2^(bit_len/2).
        let s_bits = adjusted_s.count_leading_zeros();
        let bit_len = internal_width - s_bits;
        let mut x = APInt::from_u64(internal_width, 1);
        x = x.shl((bit_len + 1) / 2);

        // Newton iterations (15 should be enough for 53-bit precision).
        for _ in 0..15 {
            let q = adjusted_s.udiv(&x);
            let sum = x.add(&q);
            let x_next = sum.lshr(1);
            if x_next.eq(&x) {
                break;
            }
            x = x_next;
        }

        let mut sig_words = x.get_raw_data().to_vec();
        while sig_words.len() > 1 && sig_words.last() == Some(&0) {
            sig_words.pop();
        }

        // Result unbiased exponent:
        // sqrt(sig << extra_shift) * 2^((E - m)/2) / 2^(extra_shift/2)
        // = sqrt(sig) * 2^(extra_shift/2) * 2^((E - m)/2) / 2^(extra_shift/2)
        // = sqrt(sig) * 2^((E - m)/2)
        //
        // To represent as result_sig * 2^(ub - m):
        // result_sig * 2^(ub - m) = sqrt(sig) * 2^((E - m)/2)
        // With result_sig = sqrt(sig) << (extra_shift/2) (approximately),
        // ub = (E + m)/2.
        //
        // But we computed x = floor(sqrt(sig << extra_shift)).
        // x ≈ sqrt(sig) << (extra_shift/2).
        // The normalize_and_round function will align this to the target precision.
        // The initial exponent accounts for the extra_shift scaling:
        // x = sqrt(sig << extra_shift) ≈ sqrt(sig) << (extra_shift/2)
        // Value = x * 2^(ub - m) = sqrt(sig) * 2^((E - m)/2)
        // x contributes 2^(target_bits) extra, so subtract that from ub.
        let mut unbiased_exp = (adjusted_exp + self.semantics.mantissa_bits as i64) / 2
            - target_bits as i64;

        let (biased_exp, inexact) = Self::normalize_and_round(
            &mut sig_words,
            &mut unbiased_exp,
            target_bits,
            rm,
            false,
            &self.semantics,
        );

        let mut result = APFloat {
            semantics: self.semantics,
            sign: false,
            exponent: biased_exp,
            significand: sig_words,
            flags: if inexact { APFLOAT_INEXACT } else { APFLOAT_OK },
        };

        while result.significand.len() > 1 && result.significand.last() == Some(&0) {
            result.significand.pop();
        }

        result
    }

    // ------------------------------------------------------------------
    // Internal: addition/subtraction of finite non-zero values
    // ------------------------------------------------------------------

    fn add_sub_internal(&self, other: &APFloat, _is_sub: bool, rm: RoundingMode) -> APFloat {
        let target_bits = self.semantics.total_significand_bits();
        let internal_width = target_bits * 4 + 10;

        let exp_a = self.unbiased_exponent();
        let exp_b = other.unbiased_exponent();

        // Determine which has larger magnitude.
        let (larger, smaller, larger_exp) = if exp_a > exp_b
            || (exp_a == exp_b
                && self
                    .sig_apint(internal_width)
                    .uge(&other.sig_apint(internal_width)))
        {
            (self, other, exp_a)
        } else {
            (other, self, exp_b)
        };

        let exp_diff = larger_exp - smaller.unbiased_exponent();

        let larger_sig = larger.sig_apint(internal_width);
        let smaller_sig = smaller.sig_apint(internal_width);

        // Align: shift smaller significand right.
        let mut smaller_words = smaller_sig.get_raw_data().to_vec();
        let _sticky = shr_vec(&mut smaller_words, exp_diff as u32);
        let aligned_smaller = APInt::from_words(internal_width, smaller_words);

        let result_sign: bool;
        let sum_words: Vec<u64>;

        let _eff_sub = larger.sign != smaller.sign;
        // is_sub parameter: for addition, effective op is (sign_a == sign_b)?
        // Actually, the caller handles signs via `sub` which negates other.
        // Here we just do sign-magnitude add/sub based on signs.
        if larger.sign == smaller.sign {
            // Same sign: add magnitudes.
            sum_words = larger_sig.add(&aligned_smaller).get_raw_data().to_vec();
            result_sign = larger.sign;
        } else {
            // Different signs: subtract smaller magnitude from larger.
            let diff = larger_sig.sub(&aligned_smaller);
            sum_words = diff.get_raw_data().to_vec();
            result_sign = larger.sign;
        }

        let mut sig = sum_words;
        while sig.len() > 1 && sig.last() == Some(&0) {
            sig.pop();
        }

        if sig.iter().all(|&w| w == 0) {
            let mut result = Self::get_zero(self.semantics, false);
            result.flags = APFLOAT_OK;
            return result;
        }

        let mut larger_unbiased_exp = larger_exp;

        let (biased_exp, inexact) = Self::normalize_and_round(
            &mut sig,
            &mut larger_unbiased_exp,
            target_bits,
            rm,
            result_sign,
            &self.semantics,
        );

        let mut result = APFloat {
            semantics: self.semantics,
            sign: result_sign,
            exponent: biased_exp,
            significand: sig,
            flags: if inexact { APFLOAT_INEXACT } else { APFLOAT_OK },
        };

        if biased_exp == self.semantics.largest_exponent && !result.is_nan() {
            result.flags |= APFLOAT_OVERFLOW | APFLOAT_INEXACT;
        }
        if biased_exp == 0 && !result.is_zero() && inexact {
            result.flags |= APFLOAT_UNDERFLOW;
        }

        while result.significand.len() > 1 && result.significand.last() == Some(&0) {
            result.significand.pop();
        }

        result
    }

    /// Normalize significand words and round to target_bits.
    /// `unbiased_exp` is updated in place.
    /// Returns (biased_exponent, inexact).
    fn normalize_and_round(
        sig_words: &mut Vec<u64>,
        unbiased_exp: &mut i64,
        target_bits: u32,
        rm: RoundingMode,
        sign: bool,
        sem: &FloatSemantics,
    ) -> (u32, bool) {
        // Early return for zero.
        if sig_words.is_empty() || sig_words.iter().all(|&w| w == 0) {
            return (0, false);
        }

        // Normalize: shift the significand so its MSB is at target_bits - 1.
        // If the MSB is above target_bits-1 (e.g., after addition carry),
        // shift right and increase exponent. If below, shift left and decrease.
        let msb = sig_msb_position(sig_words);
        let target_msb = target_bits.saturating_sub(1);
        if msb > target_msb {
            // Shift right and increase exponent.
            let shift = msb - target_msb;
            let _ = shr_vec(sig_words, shift);
            *unbiased_exp += shift as i64;
        } else if msb < target_msb && *unbiased_exp > sem.smallest_normalized_exponent as i64 - 1 {
            // Shift left and decrease exponent, but don't go subnormal.
            let max_shift = (*unbiased_exp - (sem.smallest_normalized_exponent as i64 - 1)) as u32;
            let needed_shift = (target_msb - msb).min(max_shift);
            if needed_shift > 0 {
                shl_vec(sig_words, needed_shift);
                *unbiased_exp -= needed_shift as i64;
            }
        }

        // Compute biased exponent and then round.
        let bias = sem.exponent_bias as i64;
        let biased_raw = *unbiased_exp + bias;

        let (biased_exp, inexact) = if biased_raw <= 0 {
            // Subnormal: shift right to make biased exponent = 0, then round.
            let shift = (1 - biased_raw) as u32;
            let _ = shr_vec(sig_words, shift);
            Self::round_and_renormalize(sig_words, 0, target_bits, rm, sign, sem)
        } else if biased_raw >= sem.largest_exponent as i64 {
            // Overflow to infinity.
            sig_words.clear();
            if sem.is_x87() {
                sig_words.push(1u64 << 63);
            }
            (sem.largest_exponent, true)
        } else {
            // Normal: round with the correct biased exponent.
            Self::round_and_renormalize(sig_words, biased_raw as u32, target_bits, rm, sign, sem)
        };

        (biased_exp, inexact)
    }

    /// Round this value to an integral value (integer-valued float).
    fn round_to_integral(&self, rm: RoundingMode) -> APFloat {
        if self.is_zero() || self.is_infinity() || self.is_nan() {
            return self.clone();
        }
        if self.exponent as i64 - self.semantics.exponent_bias as i64
            >= self.semantics.fraction_bits() as i64
        {
            // Already integral.
            return self.clone();
        }
        if self.unbiased_exponent() < 0 {
            // |self| < 1; result is ±0 or ±1 depending on rounding.
            match rm {
                RoundingMode::TowardNegative => {
                    if self.sign {
                        return Self::from_f64(-1.0).convert(self.semantics, rm);
                    } else {
                        return Self::get_zero(self.semantics, false);
                    }
                }
                RoundingMode::TowardPositive => {
                    if self.sign {
                        return Self::get_zero(self.semantics, true);
                    } else {
                        return Self::from_f64(1.0).convert(self.semantics, rm);
                    }
                }
                RoundingMode::NearestTiesToAway => {
                    let half = Self::from_f64(0.5).convert(self.semantics, rm);
                    let abs_self = self.abs();
                    if abs_self.compare(&half) == Some(Ordering::Greater) {
                        let one = Self::from_f64(1.0).convert(self.semantics, rm);
                        return if self.sign { one.negate() } else { one };
                    } else {
                        return Self::get_zero(self.semantics, self.sign);
                    }
                }
                _ => {
                    return Self::get_zero(self.semantics, self.sign);
                }
            }
        }

        // The fractional bits are bits below position (fraction_bits - unbiased_exp).
        let frac_pos = (self.semantics.fraction_bits() as i64 - self.unbiased_exponent()) as u32;
        // We need to mask off frac_pos low bits.
        let mut result = self.clone();
        let mask_shift = frac_pos;
        let _sticky = shr_vec(&mut result.significand, mask_shift);
        shl_vec(&mut result.significand, mask_shift);

        // The sticky bit tells us if we need to round.
        let _sticky = false; // We already have the rounding info
                             // For simplicity, re-round.
        let target_bits = self.semantics.total_significand_bits();
        let (new_exp, _inexact) = Self::round_and_renormalize(
            &mut result.significand,
            result.exponent,
            target_bits,
            rm,
            result.sign,
            &self.semantics,
        );
        result.exponent = new_exp;
        result.flags = APFLOAT_OK;
        result
    }
}

// --------------------------------------------------------------------------
// Helpers for normalize_and_round
// --------------------------------------------------------------------------

#[allow(dead_code)]
fn count_leading_zeros_sig(words: &[u64]) -> u32 {
    if words.is_empty() || words.iter().all(|&w| w == 0) {
        return 0;
    }
    let mut lz = 0u32;
    for i in (0..words.len()).rev() {
        let w = words[i];
        if w == 0 {
            lz += 64;
        } else {
            lz += w.leading_zeros();
            break;
        }
    }
    lz
}

// ============================================================================
// Comparison
// ============================================================================

impl APFloat {
    /// Compare two APFloat values according to IEEE 754 total ordering.
    ///
    /// Returns `None` if either operand is NaN. Otherwise returns the ordering
    /// of the two values, where -0 and +0 compare equal.
    pub fn compare(&self, other: &APFloat) -> Option<Ordering> {
        if self.is_nan() || other.is_nan() {
            return None;
        }

        let self_zero = self.is_zero();
        let other_zero = other.is_zero();

        if self_zero && other_zero {
            return Some(Ordering::Equal);
        }

        // Handle sign differences.
        if self.sign != other.sign {
            if self.sign {
                // self is negative, other is positive.
                return Some(Ordering::Less);
            } else {
                return Some(Ordering::Greater);
            }
        }

        // Same sign.
        let negate_result = self.sign; // if both negative, order is reversed.

        if self_zero {
            return Some(if negate_result {
                Ordering::Greater
            } else {
                Ordering::Less
            });
        }
        if other_zero {
            return Some(if negate_result {
                Ordering::Less
            } else {
                Ordering::Greater
            });
        }

        let self_inf = self.is_infinity();
        let other_inf = other.is_infinity();
        if self_inf && other_inf {
            return Some(Ordering::Equal);
        }
        if self_inf {
            return Some(if negate_result {
                Ordering::Less
            } else {
                Ordering::Greater
            });
        }
        if other_inf {
            return Some(if negate_result {
                Ordering::Greater
            } else {
                Ordering::Less
            });
        }

        // Both finite non-zero, same sign.
        let exp_a = self.unbiased_exponent();
        let exp_b = other.unbiased_exponent();

        let ord = if exp_a != exp_b {
            if exp_a > exp_b {
                Ordering::Greater
            } else {
                Ordering::Less
            }
        } else {
            // Compare significands.
            let internal_width = self.semantics.total_significand_bits() * 2;
            let sig_a = self.sig_apint(internal_width);
            let sig_b = other.sig_apint(internal_width);
            if sig_a.ugt(&sig_b) {
                Ordering::Greater
            } else if sig_a.ult(&sig_b) {
                Ordering::Less
            } else {
                Ordering::Equal
            }
        };

        if negate_result {
            Some(ord.reverse())
        } else {
            Some(ord)
        }
    }

    /// Exact bit-level identity (including NaN payload and sign).
    pub fn bitwise_identical(&self, other: &APFloat) -> bool {
        self.sign == other.sign
            && self.exponent == other.exponent
            && self.significand == other.significand
    }
}

// ============================================================================
// Conversion
// ============================================================================

impl APFloat {
    /// Convert to `f64` with the given rounding mode.
    pub fn to_f64(&self, rm: RoundingMode) -> f64 {
        // Convert to IEEE_DOUBLE semantics first, then extract bits.
        if self.semantics == IEEE_DOUBLE {
            self.to_native_f64()
        } else {
            let converted = self.convert(IEEE_DOUBLE, rm);
            converted.to_native_f64()
        }
    }

    /// Convert to `f32` with the given rounding mode.
    pub fn to_f32(&self, rm: RoundingMode) -> f32 {
        if self.semantics == IEEE_SINGLE {
            self.to_native_f32()
        } else {
            let converted = self.convert(IEEE_SINGLE, rm);
            converted.to_native_f32()
        }
    }

    /// Fast path: extract f64 bits from APFloat in IEEE_DOUBLE semantics.
    fn to_native_f64(&self) -> f64 {
        let mut bits: u64 = 0;

        // Sign.
        if self.sign {
            bits |= 1u64 << 63;
        }

        // Exponent (11 bits).
        let exp = (self.exponent as u64) & 0x7FF;
        bits |= exp << 52;

        // Mantissa (52 bits).
        if self.exponent == IEEE_DOUBLE.largest_exponent {
            // Inf or NaN: raw mantissa (without implicit bit).
            if !self.sig_is_zero() {
                let mut mantissa: u64 = 0;
                for (i, &w) in self.significand.iter().enumerate() {
                    mantissa |= w << (i * 64);
                }
                mantissa &= (1u64 << 52) - 1;
                bits |= mantissa;
            }
        } else if self.exponent == 0 {
            // Zero or subnormal: raw mantissa.
            if !self.sig_is_zero() {
                let mut mantissa: u64 = 0;
                for (i, &w) in self.significand.iter().enumerate() {
                    mantissa |= w << (i * 64);
                }
                mantissa &= (1u64 << 52) - 1;
                bits |= mantissa;
            }
        } else {
            // Normal: remove implicit leading bit.
            let mut mantissa: u64 = 0;
            for (i, &w) in self.significand.iter().enumerate() {
                mantissa |= w << (i * 64);
            }
            // Clear the implicit bit.
            mantissa &= (1u64 << 52) - 1;
            bits |= mantissa;
        }

        f64::from_bits(bits)
    }

    /// Fast path: extract f32 bits from APFloat in IEEE_SINGLE semantics.
    fn to_native_f32(&self) -> f32 {
        let mut bits: u32 = 0;

        if self.sign {
            bits |= 1u32 << 31;
        }

        let exp = (self.exponent & 0xFF) as u32;
        bits |= exp << 23;

        if self.exponent == IEEE_SINGLE.largest_exponent {
            if !self.sig_is_zero() {
                let mantissa = self.significand.first().copied().unwrap_or(0) as u32;
                bits |= mantissa & ((1u32 << 23) - 1);
            }
        } else if self.exponent == 0 {
            if !self.sig_is_zero() {
                let mantissa = self.significand.first().copied().unwrap_or(0) as u32;
                bits |= mantissa & ((1u32 << 23) - 1);
            }
        } else {
            let mantissa = self.significand.first().copied().unwrap_or(0) as u32;
            bits |= mantissa & ((1u32 << 23) - 1);
        }

        f32::from_bits(bits)
    }

    /// Bit-level representation as an APInt.
    pub fn to_apint(&self) -> APInt {
        let sem = &self.semantics;
        let bit_width = sem.bits;

        // Build the IEEE 754 bit pattern.
        let mut result = APInt::get_zero(bit_width);

        // Sign bit.
        if self.sign {
            result.set_bit(bit_width - 1);
        }

        // Exponent field.
        for i in 0..sem.exponent_bits {
            if (self.exponent >> i) & 1 != 0 {
                result.set_bit(sem.mantissa_bits + i);
            }
        }

        // Mantissa field.
        let frac_bits = sem.fraction_bits();
        for i in 0..frac_bits {
            if self.sig_get_bit(i) {
                result.set_bit(i);
            }
        }

        // For x87, the explicit integer bit is bit 63 of the mantissa field,
        // which is included in mantissa_bits=64. So bit 63 is part of the
        // fraction_bits=63 range? Actually fraction_bits = mantissa_bits-1 = 63.
        // The explicit integer bit is at position 63, which is NOT in the
        // fraction_bits range (0..63). So we need to set it separately.
        if sem.is_x87() && !self.is_zero() {
            // For x87, the explicit integer bit is at position 63.
            if self.is_infinity() || self.is_nan() || self.exponent > 0 {
                result.set_bit(63);
            }
            // For subnormals/zero, bit 63 is 0.
        }

        result
    }

    /// Convert to a string in C99 hex-float format with the given precision
    /// and padding.
    pub fn to_string(&self, precision: usize, max_padding: usize) -> String {
        if self.is_nan() {
            return "nan".to_string();
        }
        if self.is_infinity() {
            return if self.sign { "-inf" } else { "inf" }.to_string();
        }
        if self.is_zero() {
            return if self.sign { "-0.0" } else { "0.0" }.to_string();
        }

        // Use hex-float format: [-]0xh.hhhhp±d
        let sign_str = if self.sign { "-" } else { "" };
        let unbiased = self.unbiased_exponent();

        // Extract the hex digits from the significand.
        let frac_bits = self.semantics.fraction_bits();
        // We'll generate enough hex digits.
        let hex_digits = (precision.max(1) + 3) / 4;
        let mut hex_str = String::new();

        // Get the integer part (1 or 0 for subnormals).
        let int_bit = if self.exponent > 0 || self.semantics.is_x87() {
            1u8
        } else {
            0u8
        };

        // Extract fraction bits into hex digits.
        let mut frac_accum: u64 = 0;
        let mut bits_collected = 0u32;
        for i in (0..frac_bits).rev() {
            frac_accum = (frac_accum << 1) | (if self.sig_get_bit(i) { 1 } else { 0 });
            bits_collected += 1;
            if bits_collected == 4 {
                hex_str.push(char::from_digit(frac_accum as u32, 16).unwrap_or('0'));
                frac_accum = 0;
                bits_collected = 0;
            }
        }
        if bits_collected > 0 {
            frac_accum <<= 4 - bits_collected;
            hex_str.push(char::from_digit(frac_accum as u32, 16).unwrap_or('0'));
        }

        // Pad or truncate.
        if hex_str.len() < hex_digits {
            let pad = hex_digits - hex_str.len();
            for _ in 0..pad.min(max_padding) {
                hex_str.push('0');
            }
        } else if hex_str.len() > hex_digits {
            hex_str.truncate(hex_digits);
        }

        format!("{}0x{}.{}p{}", sign_str, int_bit, hex_str, unbiased)
    }

    /// Convert this value to another floating-point semantics, rounding as
    /// specified.
    pub fn convert(&self, new_sem: FloatSemantics, rm: RoundingMode) -> APFloat {
        if self.semantics == new_sem {
            return self.clone();
        }

        // Handle special values.
        if self.is_nan() {
            // Convert NaN payload. For simplicity, create a new qNaN.
            let mut result = Self::get_qnan(new_sem, self.sign);
            result.flags = self.flags;
            return result;
        }
        if self.is_infinity() {
            let mut result = Self::get_inf(new_sem, self.sign);
            result.flags = APFLOAT_OK;
            return result;
        }
        if self.is_zero() {
            let mut result = Self::get_zero(new_sem, self.sign);
            result.flags = APFLOAT_OK;
            return result;
        }

        // PPC double-double handling.
        if new_sem.is_ppc_double_double() {
            return self.convert_to_ppc_double_double(rm);
        }
        if self.semantics.is_ppc_double_double() {
            return self.convert_from_ppc_double_double(new_sem, rm);
        }

        // General conversion: unpack, round to new precision, repack.
        let target_bits = new_sem.total_significand_bits();
        let _internal_width = target_bits.max(self.semantics.total_significand_bits()) * 4 + 10;

        let mut sig_words = self.significand.clone();
        while sig_words.len() > 1 && sig_words.last() == Some(&0) {
            sig_words.pop();
        }

        let unbiased = self.unbiased_exponent();

        // Adjust the significand to the new format's precision.
        // First normalize to a canonical position.
        let old_target = self.semantics.total_significand_bits();
        let new_target = new_sem.total_significand_bits();

        let mut cur_unbiased = unbiased;

        // Adjust significand to match the new format's precision.
        // The unbiased exponent is format-independent (value = sig * 2^(unbiased - m)),
        // but the significand must be scaled to the new mantissa width.
        // sig_new = sig_old * 2^(m_new - m_old), unbiased unchanged.
        if new_target > old_target {
            let extend_shift = new_target - old_target;
            shl_vec(&mut sig_words, extend_shift);
            // unbiased is unchanged — the shift left is compensated by the
            // larger m in the new format: value = (sig<<k) * 2^(unbiased - (m+k))
            // = sig * 2^k * 2^(unbiased - m - k) = sig * 2^(unbiased - m).
        }

        if new_target < old_target {
            // Round significand to new_target bits.
            let (carry, _inexact) =
                Self::round_significand(&mut sig_words, new_target, rm, self.sign);
            if carry {
                // Shift right and increase exponent.
                let _ = shr_vec(&mut sig_words, 1);
                let set_pos = new_target - 1;
                let w = (set_pos / 64) as usize;
                let b = set_pos % 64;
                if w >= sig_words.len() {
                    sig_words.resize(w + 1, 0);
                }
                sig_words[w] |= 1u64 << b;
                cur_unbiased += 1;
            }
        }

        // Now truncate significand to new_target bits.
        let trunc_shift = if sig_msb_position(&sig_words) >= new_target {
            sig_msb_position(&sig_words) - new_target + 1
        } else {
            0
        };
        if trunc_shift > 0 {
            let _ = shr_vec(&mut sig_words, trunc_shift);
            cur_unbiased += trunc_shift as i64;
        }

        // Re-bias.
        let bias = new_sem.exponent_bias as i64;
        let biased = cur_unbiased + bias;

        let (final_exp, final_sig) = if biased <= 0 {
            // Subnormal in new format.
            let shift = (1 - biased) as u32;
            let _ = shr_vec(&mut sig_words, shift);
            (0u32, sig_words)
        } else if biased >= new_sem.largest_exponent as i64 {
            // Overflow.
            let mut inf_sig = Vec::new();
            if new_sem.is_x87() {
                inf_sig.push(1u64 << 63);
            }
            (new_sem.largest_exponent, inf_sig)
        } else {
            (biased as u32, sig_words)
        };

        // Final rounding pass on the actual significand.
        // Skip rounding for infinity (already exact).
        let mut result_sig = final_sig;
        let (final_exp, _inexact) = if final_exp == new_sem.largest_exponent && result_sig.is_empty() {
            // Infinity — no rounding needed.
            (final_exp, false)
        } else {
            Self::round_and_renormalize(
                &mut result_sig,
                final_exp,
                new_target,
                rm,
                self.sign,
                &new_sem,
            )
        };
        while result_sig.len() > 1 && result_sig.last() == Some(&0) {
            result_sig.pop();
        }

        let mut result = APFloat {
            semantics: new_sem,
            sign: self.sign,
            exponent: final_exp,
            significand: result_sig,
            flags: APFLOAT_OK,
        };

        if final_exp == new_sem.largest_exponent && !result.is_nan() {
            result.flags |= APFLOAT_OVERFLOW | APFLOAT_INEXACT;
        }

        result
    }

    /// Convert to PPC double-double format (pair of IEEE doubles).
    fn convert_to_ppc_double_double(&self, rm: RoundingMode) -> APFloat {
        // Round to the nearest double for the high part.
        let high = self.convert(IEEE_DOUBLE, rm);
        // Subtract high part from self (exact if possible).
        let high_as_self_sem = high.convert(self.semantics, rm);
        let low_exact = self.sub(&high_as_self_sem, rm);
        // Round the remainder to double.
        let low = low_exact.convert(IEEE_DOUBLE, rm);

        // Build the double-double: we store the two doubles as a pair.
        // For now, return high; the double-double handling in to_apint
        // will reconstruct.
        // Actually, we need a proper double-double representation.
        // Let's store it as: the significand encodes both doubles.
        // For simplicity, just return high with PPC semantics flags.
        let mut result = high.clone();
        result.semantics = PPC_DOUBLE_DOUBLE;
        result.flags = high.flags | low.flags;
        result
    }

    /// Convert from PPC double-double to another format.
    fn convert_from_ppc_double_double(&self, new_sem: FloatSemantics, rm: RoundingMode) -> APFloat {
        // This is a stub: proper double-double conversion requires decoding
        // the pair of doubles. For now, convert the high double.
        // In a real implementation, we'd add the high and low doubles.
        let high_bits = self.significand.first().copied().unwrap_or(0);
        let high_f64 = f64::from_bits(high_bits);
        let high = APFloat::from_f64(high_f64);
        high.convert(new_sem, rm)
    }
}

// ============================================================================
// String parsing
// ============================================================================

impl APFloat {
    /// Parse a string into an APFloat.
    ///
    /// Supports:
    /// - Decimal: `"1.5"`, `"-0.5"`, `"1.5e10"`, `"1.5e-3"`
    /// - Hex float: `"0x1.5p3"`, `"0x1.5p-3"`
    /// - Special: `"inf"`, `"infinity"`, `"nan"`, `"nan(123)"`, `"-inf"`, etc.
    pub fn from_string(s: &str, sem: FloatSemantics) -> Result<Self, String> {
        let s = s.trim().to_lowercase();
        if s.is_empty() {
            return Err("empty string".to_string());
        }

        // Check for specials.
        let (negative, rest) = if s.starts_with('-') {
            (true, &s[1..])
        } else if s.starts_with('+') {
            (false, &s[1..])
        } else {
            (false, s.as_str())
        };

        if rest == "inf" || rest == "infinity" {
            return Ok(Self::get_inf(sem, negative));
        }
        if rest.starts_with("nan") {
            let payload = if rest.starts_with("nan(") && rest.ends_with(')') {
                let inner = &rest[4..rest.len() - 1];
                if inner.is_empty() {
                    None
                } else {
                    // Parse payload as integer (decimal or hex).
                    let (radix, digits) = if inner.starts_with("0x") || inner.starts_with("0X") {
                        (16, &inner[2..])
                    } else {
                        (10, inner)
                    };
                    u64::from_str_radix(digits, radix).ok()
                }
            } else {
                None
            };
            let mut result = Self::get_nan(sem, negative, payload);
            // Check if signaling (payload without quiet bit? — for now always qNaN).
            result.flags = APFLOAT_OK;
            return Ok(result);
        }

        // Determine if this is a hex float.
        if rest.starts_with("0x") || rest.starts_with("0X") {
            return Self::from_hex_string(rest, negative, sem);
        }

        // Decimal float.
        Self::from_decimal_string(rest, negative, sem)
    }

    /// Parse a hex-float string: `0xh.hhhhp±d`
    fn from_hex_string(s: &str, negative: bool, sem: FloatSemantics) -> Result<Self, String> {
        let s = s.to_lowercase();
        // Strip "0x" prefix.
        let s = if s.starts_with("0x") {
            &s[2..]
        } else {
            s.as_str()
        };

        // Split into mantissa and exponent parts.
        let (mantissa_str, exp_str) = if let Some(pos) = s.find('p') {
            (&s[..pos], &s[pos + 1..])
        } else {
            return Err("hex float missing exponent (p)".to_string());
        };

        // Parse exponent.
        let exp: i64 = exp_str
            .parse()
            .map_err(|_| format!("invalid hex exponent: {}", exp_str))?;

        // Parse mantissa: "h.hhh" where integer part and fractional part are hex.
        let (int_part, frac_part) = if let Some(pos) = mantissa_str.find('.') {
            (&mantissa_str[..pos], &mantissa_str[pos + 1..])
        } else {
            (mantissa_str, "")
        };

        // Build the exact value as integer * 2^(4*exp - 4*frac_len).
        // Combine int_part and frac_part as a hex integer.
        let combined = if frac_part.is_empty() {
            int_part.to_string()
        } else {
            format!("{}{}", int_part, frac_part)
        };

        if combined.is_empty() {
            return Ok(Self::get_zero(sem, negative));
        }

        let int_val = APInt::from_string(&combined, 16)
            .map_err(|e| format!("invalid hex mantissa: {}", e))?;

        // The value is: (int_val / 16^frac_len) * 2^exp
        // int_val has frac_len extra hex digits, so divide by 16^frac_len = 2^(4*frac_len).
        let frac_len = frac_part.len() as i64;
        let binary_exp = exp - 4 * frac_len;

        // Build the APFloat from the integer and binary exponent.
        Self::from_apint_and_exp(&int_val, binary_exp, negative, sem)
    }

    /// Parse a decimal float string.
    fn from_decimal_string(s: &str, negative: bool, sem: FloatSemantics) -> Result<Self, String> {
        // Split into mantissa and optional exponent (e.g., "1.5e10").
        let (mantissa_str, exp_str) = if let Some(pos) = s.find('e') {
            (&s[..pos], Some(&s[pos + 1..]))
        } else {
            (s, None)
        };

        let dec_exp: i64 = if let Some(es) = exp_str {
            es.parse()
                .map_err(|_| format!("invalid exponent: {}", es))?
        } else {
            0
        };

        // Split mantissa into integer and fractional parts.
        let (int_part, frac_part) = if let Some(pos) = mantissa_str.find('.') {
            (&mantissa_str[..pos], &mantissa_str[pos + 1..])
        } else {
            (mantissa_str, "")
        };

        let digits = format!("{}{}", int_part, frac_part);
        if digits.is_empty() || digits.chars().all(|c| c == '0') {
            return Ok(Self::get_zero(sem, negative));
        }

        // Strip leading zeros.
        let digits = digits.trim_start_matches('0');
        if digits.is_empty() {
            return Ok(Self::get_zero(sem, negative));
        }

        let frac_len = frac_part.len() as i64;
        let effective_exp = dec_exp - frac_len;

        // Parse the digit string as an integer.
        let int_val = APInt::from_string(digits, 10)
            .map_err(|e| format!("invalid decimal mantissa: {}", e))?;

        // The value is: int_val * 10^effective_exp
        Self::from_decimal_apint_and_exp(&int_val, effective_exp, negative, sem)
    }

    /// Build an APFloat from an APInt mantissa and a binary exponent.
    fn from_apint_and_exp(
        mantissa: &APInt,
        binary_exp: i64,
        negative: bool,
        sem: FloatSemantics,
    ) -> Result<Self, String> {
        // The value is mantissa * 2^binary_exp.
        // We need to normalize and round.
        let words = mantissa.get_raw_data();
        let mut sig = words.to_vec();
        while sig.len() > 1 && sig.last() == Some(&0) {
            sig.pop();
        }

        if sig.iter().all(|&w| w == 0) {
            return Ok(Self::get_zero(sem, negative));
        }

        let target_bits = sem.total_significand_bits();
        let msb = sig_msb_position(&sig);
        let _sig_bits = msb + 1;

        let mut unbiased_exp = binary_exp + sem.mantissa_bits as i64;

        // Normalize significand (MSB at top of the vector).
        // Actually it's already normalized by construction.

        let (biased_exp, _inexact) = Self::normalize_and_round(
            &mut sig,
            &mut unbiased_exp,
            target_bits,
            RoundingMode::NearestTiesToEven,
            negative,
            &sem,
        );

        let mut result = APFloat {
            semantics: sem,
            sign: negative,
            exponent: biased_exp,
            significand: sig,
            flags: APFLOAT_OK,
        };

        while result.significand.len() > 1 && result.significand.last() == Some(&0) {
            result.significand.pop();
        }

        Ok(result)
    }

    /// Build an APFloat from a decimal integer and exponent of 10.
    fn from_decimal_apint_and_exp(
        mantissa: &APInt,
        dec_exp: i64,
        negative: bool,
        sem: FloatSemantics,
    ) -> Result<Self, String> {
        // We need to compute mantissa * 10^dec_exp in binary.
        // This is the core of Clinger's algorithm.
        //
        // Strategy: convert to a binary fraction. Use the fact that
        // 10^k = 2^k * 5^k, so we can multiply/divide by powers of 5
        // and adjust the binary exponent accordingly.

        let target_bits = sem.total_significand_bits();
        let internal_width = target_bits * 4 + 20;

        let mut sig = mantissa.get_raw_data().to_vec();
        while sig.len() > 1 && sig.last() == Some(&0) {
            sig.pop();
        }

        if sig.iter().all(|&w| w == 0) {
            return Ok(Self::get_zero(sem, negative));
        }

        // Compute the binary value.
        // Value = mantissa * 10^dec_exp
        //        = mantissa * 2^dec_exp * 5^dec_exp

        // We'll work with the mantissa as a big integer and adjust a binary
        // exponent separately.
        let mut binary_exp: i64 = 0;

        // Convert sig to APInt with wide width.
        let mut big = APInt::from_words(internal_width, sig);

        if dec_exp > 0 {
            // Multiply by 5^dec_exp.
            let five = APInt::from_u64(internal_width, 5);
            let mut remaining = dec_exp;
            // Use exponentiation by squaring.
            let mut pow5 = five.clone();
            while remaining > 0 {
                if remaining & 1 != 0 {
                    big = big.mul(&pow5);
                    // Normalize (shift right if too wide) — APInt multiplication is
                    // wrapping for fixed width. We need full precision.
                    // Recreate with wider width if needed.
                }
                pow5 = pow5.mul(&pow5);
                remaining >>= 1;
            }
            binary_exp += dec_exp; // multiply by 2^dec_exp
        } else if dec_exp < 0 {
            // Divide by 5^(-dec_exp).
            let pos_exp = (-dec_exp) as u64;
            let five = APInt::from_u64(internal_width, 5);
            let mut remaining = pos_exp;
            let mut pow5 = five.clone();
            // Compute 5^pos_exp.
            let mut pow5_val = APInt::from_u64(internal_width, 1);
            while remaining > 0 {
                if remaining & 1 != 0 {
                    pow5_val = pow5_val.mul(&pow5);
                }
                pow5 = pow5.mul(&pow5);
                remaining >>= 1;
            }
            // We need to maintain enough precision; shift the dividend left.
            // big * 2^N / pow5_val, where N is enough bits for exact rounding.
            let extra_bits = target_bits + 10;
            // Create a wider dividend.
            let mut wide_dividend = big.zext(internal_width + extra_bits + 64);
            wide_dividend = wide_dividend.shl(extra_bits);
            let divisor = pow5_val.zext(internal_width + extra_bits + 64);
            wide_dividend = wide_dividend.udiv(&divisor);
            big = wide_dividend;
            binary_exp += dec_exp - extra_bits as i64;
        }

        // Now we have: value = big * 2^binary_exp
        let mut sig = big.get_raw_data().to_vec();
        while sig.len() > 1 && sig.last() == Some(&0) {
            sig.pop();
        }

        if sig.iter().all(|&w| w == 0) {
            return Ok(Self::get_zero(sem, negative));
        }

        let _msb = sig_msb_position(&sig);
        let mut unbiased_exp = binary_exp + sem.mantissa_bits as i64;

        let (biased_exp, _inexact) = Self::normalize_and_round(
            &mut sig,
            &mut unbiased_exp,
            target_bits,
            RoundingMode::NearestTiesToEven,
            negative,
            &sem,
        );

        let mut result = APFloat {
            semantics: sem,
            sign: negative,
            exponent: biased_exp,
            significand: sig,
            flags: APFLOAT_OK,
        };

        while result.significand.len() > 1 && result.significand.last() == Some(&0) {
            result.significand.pop();
        }

        Ok(result)
    }
}

// ============================================================================
// Trait implementations
// ============================================================================

impl fmt::Debug for APFloat {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(
            f,
            "APFloat({}, {})",
            self.semantics.name,
            self.to_string(6, 0)
        )
    }
}

impl fmt::Display for APFloat {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "{}", self.to_string(6, 0))
    }
}

impl PartialEq for APFloat {
    fn eq(&self, other: &APFloat) -> bool {
        self.compare(other) == Some(Ordering::Equal)
    }
}

// ============================================================================
// Tests
// ============================================================================

#[cfg(test)]
mod tests {
    use super::*;

    // ------------------------------------------------------------------
    // Construction tests
    // ------------------------------------------------------------------

    #[test]
    fn test_new_zero() {
        let f = APFloat::new(IEEE_DOUBLE);
        assert!(f.is_zero());
        assert!(!f.sign);
        assert_eq!(f.classify(), FloatCategory::Zero);
    }

    #[test]
    fn test_get_zero() {
        let pos = APFloat::get_zero(IEEE_DOUBLE, false);
        let neg = APFloat::get_zero(IEEE_DOUBLE, true);
        assert!(pos.is_zero());
        assert!(neg.is_zero());
        assert!(!pos.is_negative());
        assert!(neg.is_negative());
        // -0 == +0 per IEEE 754
        assert_eq!(pos.compare(&neg), Some(Ordering::Equal));
    }

    #[test]
    fn test_get_inf() {
        let pos_inf = APFloat::get_inf(IEEE_DOUBLE, false);
        let neg_inf = APFloat::get_inf(IEEE_DOUBLE, true);
        assert!(pos_inf.is_infinity());
        assert!(neg_inf.is_infinity());
        assert!(!pos_inf.is_negative());
        assert!(neg_inf.is_negative());
        assert!(pos_inf.is_positively_infinity());
        assert!(neg_inf.is_negatively_infinity());
    }

    #[test]
    fn test_get_nan() {
        let nan = APFloat::get_nan(IEEE_DOUBLE, false, None);
        assert!(nan.is_nan());
        assert!(!nan.is_signaling_nan());
        let snan = APFloat::get_snan(IEEE_DOUBLE, false);
        assert!(snan.is_nan());
        assert!(snan.is_signaling_nan());
    }

    #[test]
    fn test_get_qnan() {
        let qnan = APFloat::get_qnan(IEEE_DOUBLE, true);
        assert!(qnan.is_nan());
        assert!(!qnan.is_signaling_nan());
        assert!(qnan.is_negative());
    }

    #[test]
    fn test_get_snan() {
        let snan = APFloat::get_snan(IEEE_DOUBLE, true);
        assert!(snan.is_nan());
        assert!(snan.is_signaling_nan());
        assert!(snan.is_negative());
    }

    #[test]
    fn test_get_snan_single() {
        let snan = APFloat::get_snan(IEEE_SINGLE, false);
        assert!(snan.is_nan());
        assert!(snan.is_signaling_nan());
    }

    #[test]
    fn test_largest() {
        let largest = APFloat::largest(IEEE_DOUBLE, false);
        assert!(largest.is_normal());
        assert!(!largest.is_negative());
        // Convert to f64 and check.
        let val = largest.to_f64(RoundingMode::NearestTiesToEven);
        assert_eq!(val, f64::MAX);
    }

    #[test]
    fn test_smallest_normalized() {
        let smallest = APFloat::smallest_normalized(IEEE_DOUBLE, false);
        assert!(smallest.is_normal());
        let val = smallest.to_f64(RoundingMode::NearestTiesToEven);
        assert_eq!(val, f64::MIN_POSITIVE);
    }

    // ------------------------------------------------------------------
    // f32/f64 round-trip tests
    // ------------------------------------------------------------------

    #[test]
    fn test_f64_roundtrip_zero() {
        let f = APFloat::from_f64(0.0);
        assert!(f.is_zero());
        assert_eq!(f.to_f64(RoundingMode::NearestTiesToEven), 0.0);
    }

    #[test]
    fn test_f64_roundtrip_neg_zero() {
        let f = APFloat::from_f64(-0.0);
        assert!(f.is_zero());
        assert!(f.is_negative());
        assert_eq!(f.to_f64(RoundingMode::NearestTiesToEven), -0.0);
    }

    #[test]
    fn test_f64_roundtrip_one() {
        let f = APFloat::from_f64(1.0);
        assert!(!f.is_zero());
        assert_eq!(f.to_f64(RoundingMode::NearestTiesToEven), 1.0);
    }

    #[test]
    fn test_f64_roundtrip_neg_one() {
        let f = APFloat::from_f64(-1.0);
        assert_eq!(f.to_f64(RoundingMode::NearestTiesToEven), -1.0);
    }

    #[test]
    fn test_f64_roundtrip_pi() {
        let f = APFloat::from_f64(std::f64::consts::PI);
        assert_eq!(
            f.to_f64(RoundingMode::NearestTiesToEven),
            std::f64::consts::PI
        );
    }

    #[test]
    fn test_f64_roundtrip_inf() {
        let f = APFloat::from_f64(f64::INFINITY);
        assert!(f.is_infinity());
        assert!(!f.is_negative());
        assert_eq!(f.to_f64(RoundingMode::NearestTiesToEven), f64::INFINITY);
    }

    #[test]
    fn test_f64_roundtrip_neg_inf() {
        let f = APFloat::from_f64(f64::NEG_INFINITY);
        assert!(f.is_infinity());
        assert!(f.is_negative());
        assert_eq!(f.to_f64(RoundingMode::NearestTiesToEven), f64::NEG_INFINITY);
    }

    #[test]
    fn test_f64_roundtrip_nan() {
        let f = APFloat::from_f64(f64::NAN);
        assert!(f.is_nan());
        let result = f.to_f64(RoundingMode::NearestTiesToEven);
        assert!(result.is_nan());
    }

    #[test]
    fn test_f32_roundtrip_zero() {
        let f = APFloat::from_f32(0.0f32);
        assert!(f.is_zero());
        assert_eq!(f.to_f32(RoundingMode::NearestTiesToEven), 0.0f32);
    }

    #[test]
    fn test_f32_roundtrip_one() {
        let f = APFloat::from_f32(1.0f32);
        assert_eq!(f.to_f32(RoundingMode::NearestTiesToEven), 1.0f32);
    }

    #[test]
    fn test_f32_roundtrip_inf() {
        let f = APFloat::from_f32(f32::INFINITY);
        assert!(f.is_infinity());
        assert_eq!(f.to_f32(RoundingMode::NearestTiesToEven), f32::INFINITY);
    }

    #[test]
    fn test_f64_roundtrip_subnormal() {
        let val = f64::from_bits(1); // smallest positive subnormal
        let f = APFloat::from_f64(val);
        assert!(f.is_subnormal());
        assert_eq!(f.to_f64(RoundingMode::NearestTiesToEven), val);
    }

    #[test]
    fn test_f64_roundtrip_max() {
        let f = APFloat::from_f64(f64::MAX);
        assert_eq!(f.to_f64(RoundingMode::NearestTiesToEven), f64::MAX);
    }

    #[test]
    fn test_f64_roundtrip_min_positive() {
        let f = APFloat::from_f64(f64::MIN_POSITIVE);
        assert!(f.is_normal());
        assert_eq!(f.to_f64(RoundingMode::NearestTiesToEven), f64::MIN_POSITIVE);
    }

    // ------------------------------------------------------------------
    // Arithmetic tests
    // ------------------------------------------------------------------

    #[test]
    fn test_add_basic() {
        let a = APFloat::from_f64(1.5);
        let b = APFloat::from_f64(2.5);
        let sum = a.add(&b, RoundingMode::NearestTiesToEven);
        assert_eq!(sum.to_f64(RoundingMode::NearestTiesToEven), 4.0);
    }

    #[test]
    fn test_add_negatives() {
        let a = APFloat::from_f64(-1.5);
        let b = APFloat::from_f64(-2.5);
        let sum = a.add(&b, RoundingMode::NearestTiesToEven);
        assert_eq!(sum.to_f64(RoundingMode::NearestTiesToEven), -4.0);
    }

    #[test]
    fn test_add_mixed_signs() {
        let a = APFloat::from_f64(3.0);
        let b = APFloat::from_f64(-1.0);
        let sum = a.add(&b, RoundingMode::NearestTiesToEven);
        assert_eq!(sum.to_f64(RoundingMode::NearestTiesToEven), 2.0);
    }

    #[test]
    fn test_sub_basic() {
        let a = APFloat::from_f64(5.0);
        let b = APFloat::from_f64(3.0);
        let diff = a.sub(&b, RoundingMode::NearestTiesToEven);
        assert_eq!(diff.to_f64(RoundingMode::NearestTiesToEven), 2.0);
    }

    #[test]
    fn test_mul_basic() {
        let a = APFloat::from_f64(2.0);
        let b = APFloat::from_f64(3.0);
        let prod = a.mul(&b, RoundingMode::NearestTiesToEven);
        assert_eq!(prod.to_f64(RoundingMode::NearestTiesToEven), 6.0);
    }

    #[test]
    fn test_mul_negatives() {
        let a = APFloat::from_f64(-2.0);
        let b = APFloat::from_f64(3.0);
        let prod = a.mul(&b, RoundingMode::NearestTiesToEven);
        assert_eq!(prod.to_f64(RoundingMode::NearestTiesToEven), -6.0);
    }

    #[test]
    fn test_div_basic() {
        let a = APFloat::from_f64(6.0);
        let b = APFloat::from_f64(2.0);
        let quot = a.div(&b, RoundingMode::NearestTiesToEven);
        assert_eq!(quot.to_f64(RoundingMode::NearestTiesToEven), 3.0);
    }

    #[test]
    fn test_div_negatives() {
        let a = APFloat::from_f64(-6.0);
        let b = APFloat::from_f64(2.0);
        let quot = a.div(&b, RoundingMode::NearestTiesToEven);
        assert_eq!(quot.to_f64(RoundingMode::NearestTiesToEven), -3.0);
    }

    // ------------------------------------------------------------------
    // Special cases
    // ------------------------------------------------------------------

    #[test]
    fn test_add_inf_inf_same_sign() {
        let a = APFloat::get_inf(IEEE_DOUBLE, false);
        let b = APFloat::get_inf(IEEE_DOUBLE, false);
        let sum = a.add(&b, RoundingMode::NearestTiesToEven);
        assert!(sum.is_infinity());
        assert!(!sum.is_negative());
    }

    #[test]
    fn test_add_inf_neg_inf() {
        let a = APFloat::get_inf(IEEE_DOUBLE, false);
        let b = APFloat::get_inf(IEEE_DOUBLE, true);
        let sum = a.add(&b, RoundingMode::NearestTiesToEven);
        assert!(sum.is_nan());
        assert!(sum.flags & APFLOAT_INVALID_OP != 0);
    }

    #[test]
    fn test_mul_inf_zero() {
        let a = APFloat::get_inf(IEEE_DOUBLE, false);
        let b = APFloat::get_zero(IEEE_DOUBLE, false);
        let prod = a.mul(&b, RoundingMode::NearestTiesToEven);
        assert!(prod.is_nan());
        assert!(prod.flags & APFLOAT_INVALID_OP != 0);
    }

    #[test]
    fn test_div_zero_by_zero() {
        let a = APFloat::get_zero(IEEE_DOUBLE, false);
        let b = APFloat::get_zero(IEEE_DOUBLE, false);
        let quot = a.div(&b, RoundingMode::NearestTiesToEven);
        assert!(quot.is_nan());
        assert!(quot.flags & APFLOAT_INVALID_OP != 0);
    }

    #[test]
    fn test_div_by_zero() {
        let a = APFloat::from_f64(1.0);
        let b = APFloat::get_zero(IEEE_DOUBLE, false);
        let quot = a.div(&b, RoundingMode::NearestTiesToEven);
        assert!(quot.is_infinity());
        assert!(quot.flags & APFLOAT_DIV_BY_ZERO != 0);
    }

    #[test]
    fn test_sqrt_negative() {
        let a = APFloat::from_f64(-4.0);
        let r = a.sqrt(RoundingMode::NearestTiesToEven);
        assert!(r.is_nan());
        assert!(r.flags & APFLOAT_INVALID_OP != 0);
    }

    #[test]
    fn test_sqrt_zero() {
        let a = APFloat::from_f64(0.0);
        let r = a.sqrt(RoundingMode::NearestTiesToEven);
        assert!(r.is_zero());
    }

    #[test]
    fn test_sqrt_four() {
        let a = APFloat::from_f64(4.0);
        let r = a.sqrt(RoundingMode::NearestTiesToEven);
        assert_eq!(r.to_f64(RoundingMode::NearestTiesToEven), 2.0);
    }

    // ------------------------------------------------------------------
    // NaN propagation
    // ------------------------------------------------------------------

    #[test]
    fn test_nan_add_propagation() {
        let nan = APFloat::get_qnan(IEEE_DOUBLE, false);
        let val = APFloat::from_f64(1.0);
        let sum = nan.add(&val, RoundingMode::NearestTiesToEven);
        assert!(sum.is_nan());
    }

    #[test]
    fn test_snan_add_invalid() {
        let snan = APFloat::get_snan(IEEE_DOUBLE, false);
        let val = APFloat::from_f64(1.0);
        let sum = snan.add(&val, RoundingMode::NearestTiesToEven);
        assert!(sum.is_nan());
        assert!(sum.flags & APFLOAT_INVALID_OP != 0);
    }

    // ------------------------------------------------------------------
    // Comparison edge cases
    // ------------------------------------------------------------------

    #[test]
    fn test_compare_nan_returns_none() {
        let nan = APFloat::get_qnan(IEEE_DOUBLE, false);
        let val = APFloat::from_f64(1.0);
        assert_eq!(nan.compare(&val), None);
        assert_eq!(val.compare(&nan), None);
        assert_eq!(nan.compare(&nan), None);
    }

    #[test]
    fn test_compare_zero_equal() {
        let pos = APFloat::get_zero(IEEE_DOUBLE, false);
        let neg = APFloat::get_zero(IEEE_DOUBLE, true);
        assert_eq!(pos.compare(&neg), Some(Ordering::Equal));
    }

    #[test]
    fn test_compare_neg_less_than_pos() {
        let neg = APFloat::from_f64(-1.0);
        let pos = APFloat::from_f64(1.0);
        assert_eq!(neg.compare(&pos), Some(Ordering::Less));
        assert_eq!(pos.compare(&neg), Some(Ordering::Greater));
    }

    #[test]
    fn test_compare_infinity() {
        let inf = APFloat::get_inf(IEEE_DOUBLE, false);
        let max = APFloat::from_f64(f64::MAX);
        assert_eq!(max.compare(&inf), Some(Ordering::Less));
        assert_eq!(inf.compare(&max), Some(Ordering::Greater));
    }

    #[test]
    fn test_bitwise_identical() {
        let a = APFloat::from_f64(1.0);
        let b = APFloat::from_f64(1.0);
        assert!(a.bitwise_identical(&b));
    }

    #[test]
    fn test_bitwise_not_identical() {
        let a = APFloat::from_f64(1.0);
        let b = APFloat::from_f64(2.0);
        assert!(!a.bitwise_identical(&b));
    }

    // ------------------------------------------------------------------
    // String parsing
    // ------------------------------------------------------------------

    #[test]
    fn test_parse_decimal_simple() {
        let f = APFloat::from_string("1.5", IEEE_DOUBLE).unwrap();
        assert_eq!(f.to_f64(RoundingMode::NearestTiesToEven), 1.5);
    }

    #[test]
    fn test_parse_decimal_negative() {
        let f = APFloat::from_string("-0.5", IEEE_DOUBLE).unwrap();
        assert_eq!(f.to_f64(RoundingMode::NearestTiesToEven), -0.5);
    }

    #[test]
    fn test_parse_decimal_exponent() {
        let f = APFloat::from_string("1.5e10", IEEE_DOUBLE).unwrap();
        assert_eq!(f.to_f64(RoundingMode::NearestTiesToEven), 1.5e10);
    }

    #[test]
    fn test_parse_decimal_negative_exponent() {
        let f = APFloat::from_string("1.5e-3", IEEE_DOUBLE).unwrap();
        assert!((f.to_f64(RoundingMode::NearestTiesToEven) - 0.0015).abs() < 1e-10);
    }

    #[test]
    fn test_parse_hex_float() {
        let f = APFloat::from_string("0x1.5p3", IEEE_DOUBLE).unwrap();
        // 0x1.5 = 1 + 5/16 = 1.3125, * 2^3 = 10.5
        assert_eq!(f.to_f64(RoundingMode::NearestTiesToEven), 10.5);
    }

    #[test]
    fn test_parse_hex_float_negative_exp() {
        let f = APFloat::from_string("0x1.0p-2", IEEE_DOUBLE).unwrap();
        assert_eq!(f.to_f64(RoundingMode::NearestTiesToEven), 0.25);
    }

    #[test]
    fn test_parse_inf() {
        let f = APFloat::from_string("inf", IEEE_DOUBLE).unwrap();
        assert!(f.is_infinity());
        assert!(!f.is_negative());
    }

    #[test]
    fn test_parse_neg_inf() {
        let f = APFloat::from_string("-inf", IEEE_DOUBLE).unwrap();
        assert!(f.is_infinity());
        assert!(f.is_negative());
    }

    #[test]
    fn test_parse_nan() {
        let f = APFloat::from_string("nan", IEEE_DOUBLE).unwrap();
        assert!(f.is_nan());
    }

    #[test]
    fn test_parse_nan_payload() {
        let f = APFloat::from_string("nan(42)", IEEE_DOUBLE).unwrap();
        assert!(f.is_nan());
    }

    #[test]
    fn test_parse_neg_nan() {
        let f = APFloat::from_string("-nan", IEEE_DOUBLE).unwrap();
        assert!(f.is_nan());
        assert!(f.is_negative());
    }

    // ------------------------------------------------------------------
    // Conversion between semantics
    // ------------------------------------------------------------------

    #[test]
    fn test_convert_double_to_single() {
        let d = APFloat::from_f64(std::f64::consts::PI);
        let s = d.convert(IEEE_SINGLE, RoundingMode::NearestTiesToEven);
        let f32_val = s.to_f32(RoundingMode::NearestTiesToEven);
        assert!((f32_val as f64 - std::f64::consts::PI).abs() < 1e-6);
    }

    #[test]
    fn test_convert_single_to_double() {
        let s = APFloat::from_f32(1.5f32);
        let d = s.convert(IEEE_DOUBLE, RoundingMode::NearestTiesToEven);
        assert_eq!(d.to_f64(RoundingMode::NearestTiesToEven), 1.5);
    }

    #[test]
    fn test_convert_overflow_to_inf() {
        let large = APFloat::from_f64(1e300);
        let s = large.convert(IEEE_SINGLE, RoundingMode::NearestTiesToEven);
        assert!(s.is_infinity());
    }

    // ------------------------------------------------------------------
    // Classification tests
    // ------------------------------------------------------------------

    #[test]
    fn test_classify_normal() {
        let f = APFloat::from_f64(1.0);
        assert_eq!(f.classify(), FloatCategory::Normal);
        assert!(f.is_normal());
        assert!(f.is_finite());
        assert!(f.is_non_zero());
    }

    #[test]
    fn test_classify_subnormal() {
        let f = APFloat::from_f64(f64::from_bits(1));
        assert_eq!(f.classify(), FloatCategory::Subnormal);
        assert!(f.is_subnormal());
    }

    #[test]
    fn test_classify_zero() {
        let f = APFloat::from_f64(0.0);
        assert_eq!(f.classify(), FloatCategory::Zero);
        assert!(!f.is_non_zero());
    }

    #[test]
    fn test_classify_inf() {
        let f = APFloat::get_inf(IEEE_DOUBLE, false);
        assert_eq!(f.classify(), FloatCategory::Infinity);
    }

    #[test]
    fn test_classify_nan() {
        let f = APFloat::get_qnan(IEEE_DOUBLE, false);
        assert_eq!(f.classify(), FloatCategory::NaN);
    }

    // ------------------------------------------------------------------
    // FMA tests
    // ------------------------------------------------------------------

    #[test]
    fn test_fma_basic() {
        let a = APFloat::from_f64(2.0);
        let b = APFloat::from_f64(3.0);
        let c = APFloat::from_f64(4.0);
        let r = a.fma(&b, &c, RoundingMode::NearestTiesToEven);
        assert_eq!(r.to_f64(RoundingMode::NearestTiesToEven), 10.0);
    }

    #[test]
    fn test_fma_inf_zero() {
        let a = APFloat::get_inf(IEEE_DOUBLE, false);
        let b = APFloat::get_zero(IEEE_DOUBLE, false);
        let c = APFloat::from_f64(1.0);
        let r = a.fma(&b, &c, RoundingMode::NearestTiesToEven);
        assert!(r.is_nan());
        assert!(r.flags & APFLOAT_INVALID_OP != 0);
    }

    // ------------------------------------------------------------------
    // sqrt tests
    // ------------------------------------------------------------------

    #[test]
    fn test_sqrt_two() {
        let a = APFloat::from_f64(2.0);
        let r = a.sqrt(RoundingMode::NearestTiesToEven);
        let diff = (r.to_f64(RoundingMode::NearestTiesToEven) - 2.0f64.sqrt()).abs();
        assert!(diff < 1e-15);
    }

    #[test]
    fn test_sqrt_inf() {
        let a = APFloat::get_inf(IEEE_DOUBLE, false);
        let r = a.sqrt(RoundingMode::NearestTiesToEven);
        assert!(r.is_infinity());
        assert!(!r.is_negative());
    }

    // ------------------------------------------------------------------
    // Negate / abs / copy_sign tests
    // ------------------------------------------------------------------

    #[test]
    fn test_negate() {
        let a = APFloat::from_f64(1.0);
        let b = a.negate();
        assert_eq!(b.to_f64(RoundingMode::NearestTiesToEven), -1.0);
    }

    #[test]
    fn test_abs() {
        let a = APFloat::from_f64(-1.0);
        let b = a.abs();
        assert_eq!(b.to_f64(RoundingMode::NearestTiesToEven), 1.0);
    }

    #[test]
    fn test_copy_sign() {
        let a = APFloat::from_f64(1.0);
        let b = APFloat::from_f64(-2.0);
        let c = a.copy_sign(&b);
        assert_eq!(c.to_f64(RoundingMode::NearestTiesToEven), -1.0);
    }

    // ------------------------------------------------------------------
    // from_apint / to_apint tests
    // ------------------------------------------------------------------

    #[test]
    fn test_from_apint_zero() {
        let bits = APInt::from_u64(64, 0);
        let f = APFloat::from_apint(&bits, IEEE_DOUBLE);
        assert!(f.is_zero());
    }

    #[test]
    fn test_from_apint_one() {
        // IEEE double for 1.0: 0x3FF0000000000000
        let bits = APInt::from_u64(64, 0x3FF0000000000000u64);
        let f = APFloat::from_apint(&bits, IEEE_DOUBLE);
        assert_eq!(f.to_f64(RoundingMode::NearestTiesToEven), 1.0);
    }

    #[test]
    fn test_to_apint_roundtrip() {
        let f = APFloat::from_f64(1.5);
        let bits = f.to_apint();
        let f2 = APFloat::from_apint(&bits, IEEE_DOUBLE);
        assert!(f.bitwise_identical(&f2));
    }

    // ------------------------------------------------------------------
    // x87 tests
    // ------------------------------------------------------------------

    #[test]
    fn test_x87_zero() {
        let f = APFloat::get_zero(X87_DOUBLE_EXTENDED, false);
        assert!(f.is_zero());
    }

    #[test]
    fn test_x87_inf() {
        let f = APFloat::get_inf(X87_DOUBLE_EXTENDED, false);
        assert!(f.is_infinity());
        assert!(!f.is_negative());
    }

    #[test]
    fn test_x87_nan() {
        let f = APFloat::get_qnan(X87_DOUBLE_EXTENDED, false);
        assert!(f.is_nan());
    }

    // ------------------------------------------------------------------
    // Rounding mode tests
    // ------------------------------------------------------------------

    #[test]
    fn test_round_toward_zero() {
        let a = APFloat::from_f64(1.7);
        let b = APFloat::from_f64(1.0);
        let diff = a.sub(&b, RoundingMode::TowardZero);
        assert_eq!(diff.to_f64(RoundingMode::NearestTiesToEven), 0.7);
    }

    // ------------------------------------------------------------------
    // Exception flag tests
    // ------------------------------------------------------------------

    #[test]
    fn test_flags_clear_on_new_op() {
        let a = APFloat::get_inf(IEEE_DOUBLE, false);
        let b = APFloat::get_inf(IEEE_DOUBLE, true);
        let r = a.add(&b, RoundingMode::NearestTiesToEven);
        assert!(r.flags & APFLOAT_INVALID_OP != 0);
        // A new operation should have clean flags (except those from its own
        // computation).
        let c = APFloat::from_f64(1.0);
        assert_eq!(c.flags, APFLOAT_OK);
    }

    // ------------------------------------------------------------------
    // Display / Debug
    // ------------------------------------------------------------------

    #[test]
    fn test_display() {
        let f = APFloat::from_f64(1.0);
        let s = format!("{}", f);
        assert!(!s.is_empty());
    }

    #[test]
    fn test_debug() {
        let f = APFloat::from_f64(1.0);
        let s = format!("{:?}", f);
        assert!(s.contains("APFloat"));
    }

    // ------------------------------------------------------------------
    // PPC double-double basic tests
    // ------------------------------------------------------------------

    #[test]
    fn test_ppc_double_double_zero() {
        let f = APFloat::get_zero(PPC_DOUBLE_DOUBLE, false);
        assert!(f.is_zero());
    }

    #[test]
    fn test_ppc_double_double_inf() {
        let f = APFloat::get_inf(PPC_DOUBLE_DOUBLE, false);
        assert!(f.is_infinity());
    }
}