la_stack/

interval.rs

#![forbid(unsafe_code)]

//! Outward-rounded intervals, fixed-size determinants, and solution enclosures.
//!
//! See `REFERENCES.md` \[17\] for `FastTwoSum`, \[9-11\] for the binary64 arithmetic
//! model, and \[12\] for the Leibniz determinant identity. The column-subset
//! evaluation is specialized to this crate's small dimensions. Reference
//! \[14\] describes the broader interval standard; this module does not claim
//! IEEE 1788 conformance. The
//! [interval construction](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#outward-rounded-interval-expressions)
//! explains outward endpoints and determinant enclosures.

use crate::rounding::{compare_product_with_rounded, compare_quotient_with_rounded, two_sum_error};
use crate::{ArithmeticOperation, IntervalBound, IntervalOperand, LaError, Matrix};

/// Largest dimension supported by [`IntervalMatrix::det`] and
/// [`IntervalMatrix::det_sign`].
///
/// A subset-DP determinant needs `2^D` partial intervals. The implementation
/// reserves up to 128 entries inline, covering the geometry-oriented D ≤ 7 scope
/// without heap allocation.
pub const MAX_INTERVAL_MATRIX_DIM: usize = 7;

/// A closed finite binary64 interval `[lower, upper]`.
///
/// Construction keeps both endpoints finite and ordered. Arithmetic rounds
/// outward, so every successful result contains the exact-real result of the
/// corresponding operation on all represented inputs. Both IEEE-754 signed
/// zeros are accepted and canonicalized to `+0.0`; subnormal bounds are
/// retained.
///
/// This is a deliberately small proof-bearing surface, not a general-purpose
/// interval arithmetic package. Division requires a zero-free divisor.
///
/// # Examples
/// ```
/// use la_stack::prelude::*;
///
/// # fn main() -> Result<(), LaError> {
/// let difference = Interval::try_from_subtraction(1.0, 0.1)?;
/// let square = difference.try_square()?;
/// assert!(difference.lower() < difference.upper());
/// assert!(square.contains((1.0_f64 - 0.1).powi(2)));
/// # Ok(())
/// # }
/// ```
#[must_use]
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Interval {
    lower: f64,
    upper: f64,
}

/// Sign evidence from an outward-rounded interval determinant.
///
/// `Positive`, `Negative`, and `Zero` are proofs about every determinant
/// represented by the interval matrix. `Inconclusive` means the computed
/// enclosure overlaps zero and must not be interpreted as exact singularity.
#[must_use]
#[non_exhaustive]
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum IntervalDeterminantSign {
    /// The determinant interval is strictly greater than zero.
    Positive,
    /// The determinant interval is strictly less than zero.
    Negative,
    /// The determinant interval is exactly `[0, 0]`.
    Zero,
    /// The determinant interval contains zero and at least one nonzero value.
    Inconclusive,
}

/// Fixed-size square matrix of outward-rounded [`Interval`] entries.
///
/// Storage is the inline array `[[Interval; D]; D]`. Determinants use a
/// division-free Leibniz subset DP through D=7, so zero-containing pivot
/// intervals never require a special case and no heap allocation occurs.
/// [`solve`](Self::solve) uses interval Gaussian elimination with row pivoting
/// and explicitly reports inconclusive pivots. It has no determinant dimension
/// limit, uses O(D²) stack storage, and is intended for small systems.
///
/// # Examples
/// ```
/// use la_stack::prelude::*;
///
/// # fn main() -> Result<(), LaError> {
/// let matrix = IntervalMatrix::<3>::try_from_point_rows([
///     [0.0, 1.0, 0.0],
///     [1.0, 0.0, 0.0],
///     [0.0, 0.0, 1.0],
/// ])?;
/// assert_eq!(matrix.det_sign()?, IntervalDeterminantSign::Negative);
/// # Ok(())
/// # }
/// ```
#[must_use]
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct IntervalMatrix<const D: usize> {
    rows: [[Interval; D]; D],
}

/// Canonicalize either signed representation of real zero to `+0.0`.
#[inline]
const fn canonical_zero(value: f64) -> f64 {
    if value == 0.0 { 0.0 } else { value }
}

/// Turn a finite rounded sum into the tight adjacent-float enclosure implied by
/// its exact `FastTwoSum` residual.
#[inline]
const fn rounded_add_bounds(
    left: f64,
    right: f64,
    operation: ArithmeticOperation,
) -> Result<(f64, f64), LaError> {
    let rounded = left + right;
    if !rounded.is_finite() {
        return Err(LaError::interval_range_exhausted(operation));
    }

    let error = two_sum_error(left, right, rounded);
    // The ordered FastTwoSum residual is finite for finite operands and sum.
    let (lower, upper) = if error < 0.0 {
        (rounded.next_down(), rounded)
    } else if error > 0.0 {
        (rounded, rounded.next_up())
    } else {
        (rounded, rounded)
    };
    if !lower.is_finite() || !upper.is_finite() {
        return Err(LaError::interval_range_exhausted(operation));
    }

    Ok((canonical_zero(lower), canonical_zero(upper)))
}

/// Round only the endpoint used by a non-singleton interval sum.
#[inline]
const fn rounded_add_endpoint<const UPPER: bool>(
    left: f64,
    right: f64,
    operation: ArithmeticOperation,
) -> Result<f64, LaError> {
    let rounded = left + right;
    if !rounded.is_finite() {
        return Err(LaError::interval_range_exhausted(operation));
    }
    let error = two_sum_error(left, right, rounded);
    let endpoint = if UPPER && error > 0.0 {
        rounded.next_up()
    } else if !UPPER && error < 0.0 {
        rounded.next_down()
    } else {
        rounded
    };
    if !endpoint.is_finite() {
        return Err(LaError::interval_range_exhausted(operation));
    }
    Ok(canonical_zero(endpoint))
}

/// Turn a finite rounded product into the tight adjacent-float enclosure of the
/// exact binary64-input product.
// This kernel is called for each selected endpoint inside the determinant DP.
#[expect(
    clippy::inline_always,
    reason = "measured endpoint-call overhead in the interval determinant kernel"
)]
#[inline(always)]
const fn rounded_product_bounds(
    left: f64,
    right: f64,
    operation: ArithmeticOperation,
) -> Result<(f64, f64), LaError> {
    if left == 0.0 || right == 0.0 {
        return Ok((0.0, 0.0));
    }

    let rounded = left * right;
    if !rounded.is_finite() {
        return Err(LaError::interval_range_exhausted(operation));
    }

    let relation = compare_product_with_rounded(left, right, rounded);
    let (lower, upper) = if relation < 0 {
        (rounded.next_down(), rounded)
    } else if relation > 0 {
        (rounded, rounded.next_up())
    } else {
        (rounded, rounded)
    };
    if !lower.is_finite() || !upper.is_finite() {
        return Err(LaError::interval_range_exhausted(operation));
    }

    Ok((canonical_zero(lower), canonical_zero(upper)))
}

/// Enclose a quotient directly, including subnormal results and results rounded
/// to zero. Callers supply a nonzero denominator. An exact integer comparison
/// determines which adjacent float is needed, even when an FMA residual would
/// underflow or a reciprocal intermediate would overflow.
#[inline]
const fn rounded_quotient_bounds(
    numerator: f64,
    denominator: f64,
    operation: ArithmeticOperation,
) -> Result<(f64, f64), LaError> {
    let rounded = numerator / denominator;
    if !rounded.is_finite() {
        return Err(LaError::interval_range_exhausted(operation));
    }
    let relation = compare_quotient_with_rounded(numerator, denominator, rounded);
    let (lower, upper) = if relation < 0 {
        (rounded.next_down(), rounded)
    } else if relation > 0 {
        (rounded, rounded.next_up())
    } else {
        (rounded, rounded)
    };
    if !lower.is_finite() || !upper.is_finite() {
        return Err(LaError::interval_range_exhausted(operation));
    }
    Ok((canonical_zero(lower), canonical_zero(upper)))
}

impl Interval {
    /// Exact real zero.
    pub const ZERO: Self = Self {
        lower: 0.0,
        upper: 0.0,
    };

    /// Exact real one.
    pub const ONE: Self = Self {
        lower: 1.0,
        upper: 1.0,
    };

    /// Construct a closed interval from finite ordered bounds.
    ///
    /// Signed zero endpoints are canonicalized to `+0.0`.
    ///
    /// # Examples
    /// ```
    /// use core::assert_matches;
    /// use la_stack::prelude::*;
    ///
    /// # fn main() -> Result<(), LaError> {
    /// let range = Interval::try_new(-2.0, 3.0)?;
    /// assert!(range.contains(1.0));
    /// assert!(!range.contains(4.0));
    /// assert_matches!(
    ///     Interval::try_new(3.0, -2.0),
    ///     Err(LaError::InvertedInterval { lower: 3.0, upper: -2.0, .. })
    /// );
    /// # Ok(())
    /// # }
    /// ```
    ///
    /// # Errors
    /// Returns [`LaError::NonFinite`] when either endpoint is NaN or infinity.
    /// Returns [`LaError::InvertedInterval`] when `lower > upper`.
    #[inline]
    pub const fn try_new(lower: f64, upper: f64) -> Result<Self, LaError> {
        if !lower.is_finite() {
            return Err(LaError::non_finite_input_interval_bound(
                IntervalBound::Lower,
            ));
        }
        if !upper.is_finite() {
            return Err(LaError::non_finite_input_interval_bound(
                IntervalBound::Upper,
            ));
        }
        if lower > upper {
            return Err(LaError::inverted_interval(lower, upper));
        }
        Ok(Self::new_unchecked(lower, upper))
    }

    /// Construct a point interval from a finite binary64 value.
    ///
    /// This preserves the supplied value, including any earlier rounding.
    /// Use [`try_from_subtraction`](Self::try_from_subtraction) to enclose a
    /// subtraction before its rounding uncertainty is lost.
    ///
    /// # Examples
    /// ```
    /// use la_stack::prelude::*;
    ///
    /// # fn main() -> Result<(), LaError> {
    /// let half = Interval::point(0.5)?;
    /// assert_eq!((half.lower(), half.upper()), (0.5, 0.5));
    /// assert_eq!(half.try_add(&half)?, Interval::ONE);
    /// # Ok(())
    /// # }
    /// ```
    ///
    /// # Errors
    /// Returns [`LaError::NonFinite`] when `value` is NaN or infinity.
    #[inline]
    pub const fn point(value: f64) -> Result<Self, LaError> {
        match Self::try_new(value, value) {
            Ok(interval) => Ok(interval),
            Err(LaError::NonFinite { .. }) => Err(LaError::non_finite_input_scalar()),
            Err(error) => Err(error),
        }
    }

    /// Enclose the exact-real subtraction of two finite binary64 inputs.
    ///
    /// Unlike subtracting first and then calling [`point`](Self::point), this
    /// method preserves the rounding uncertainty introduced by the subtraction.
    ///
    /// # Examples
    /// ```
    /// use la_stack::prelude::*;
    ///
    /// # fn main() -> Result<(), LaError> {
    /// // The exact difference 1 - 2^-54 lies between adjacent binary64 values.
    /// let difference = Interval::try_from_subtraction(1.0, f64::EPSILON / 4.0)?;
    /// assert_eq!(difference.lower(), 1.0_f64.next_down());
    /// assert_eq!(difference.upper(), 1.0);
    ///
    /// // Subtracting first loses that uncertainty and produces a point at 1.
    /// let rounded = Interval::point(1.0 - f64::EPSILON / 4.0)?;
    /// assert_eq!(rounded, Interval::ONE);
    /// # Ok(())
    /// # }
    /// ```
    ///
    /// # Errors
    /// Returns [`LaError::NonFinite`] for a non-finite input, preserving whether
    /// it was the left or right operand. Returns
    /// [`LaError::IntervalRangeExhausted`] when the exact difference has no
    /// finite binary64 enclosure.
    #[inline]
    pub const fn try_from_subtraction(left: f64, right: f64) -> Result<Self, LaError> {
        if !left.is_finite() {
            return Err(LaError::non_finite_input_interval_operand(
                IntervalOperand::Left,
            ));
        }
        if !right.is_finite() {
            return Err(LaError::non_finite_input_interval_operand(
                IntervalOperand::Right,
            ));
        }
        if right == 0.0 {
            return Ok(Self::new_unchecked(left, left));
        }
        if left == 0.0 {
            return Ok(Self::new_unchecked(-right, -right));
        }
        match rounded_add_bounds(left, -right, ArithmeticOperation::IntervalSubtraction) {
            Ok((lower, upper)) => Ok(Self::new_unchecked(lower, upper)),
            Err(error) => Err(error),
        }
    }

    /// Return the finite lower bound.
    #[inline]
    #[must_use]
    pub const fn lower(self) -> f64 {
        self.lower
    }

    /// Return the finite upper bound.
    #[inline]
    #[must_use]
    pub const fn upper(self) -> f64 {
        self.upper
    }

    /// Return whether this interval contains the finite `value`.
    #[inline]
    #[must_use]
    pub const fn contains(self, value: f64) -> bool {
        value.is_finite() && self.lower <= value && value <= self.upper
    }

    /// Add two intervals with outward rounding.
    ///
    /// # Examples
    /// ```
    /// use la_stack::prelude::*;
    ///
    /// # fn main() -> Result<(), LaError> {
    /// let left = Interval::try_new(1.0, 2.0)?;
    /// let right = Interval::try_new(0.5, 1.0)?;
    /// assert_eq!(left.try_add(&right)?, Interval::try_new(1.5, 3.0)?);
    /// # Ok(())
    /// # }
    /// ```
    ///
    /// # Errors
    /// Returns [`LaError::IntervalRangeExhausted`] when the exact result range
    /// has no finite binary64 enclosure.
    #[inline]
    pub const fn try_add(&self, other: &Self) -> Result<Self, LaError> {
        self.try_add_for(other, ArithmeticOperation::IntervalAddition)
    }

    /// Multiply two intervals with outward rounding.
    ///
    /// For a square of the same represented value, prefer
    /// [`try_square`](Self::try_square), which can give a tighter enclosure.
    ///
    /// # Examples
    /// ```
    /// use la_stack::prelude::*;
    ///
    /// # fn main() -> Result<(), LaError> {
    /// let left = Interval::try_new(-2.0, 3.0)?;
    /// let right = Interval::try_new(-4.0, -1.0)?;
    /// assert_eq!(left.try_mul(&right)?, Interval::try_new(-12.0, 8.0)?);
    /// # Ok(())
    /// # }
    /// ```
    ///
    /// # Errors
    /// Returns [`LaError::IntervalRangeExhausted`] when the exact result range
    /// has no finite binary64 enclosure.
    #[inline]
    pub const fn try_mul(&self, other: &Self) -> Result<Self, LaError> {
        self.try_mul_for(other, ArithmeticOperation::IntervalMultiplication)
    }

    /// Enclose `1 / x` for every real `x` in this interval.
    ///
    /// # Examples
    /// ```
    /// use la_stack::prelude::*;
    /// # fn main() -> Result<(), LaError> {
    /// assert_eq!(Interval::try_new(-4.0, -2.0)?.try_recip()?,
    ///            Interval::try_new(-0.5, -0.25)?);
    /// # Ok(())
    /// # }
    /// ```
    ///
    /// # Errors
    /// Returns [`LaError::IntervalZeroDivisor`] if this interval contains zero,
    /// or [`LaError::IntervalRangeExhausted`] if no finite enclosure exists.
    #[inline]
    pub const fn try_recip(&self) -> Result<Self, LaError> {
        Self::ONE.try_div_for(self, ArithmeticOperation::IntervalReciprocal)
    }

    /// Enclose `x / y` for all independently represented `x` and `y`.
    ///
    /// Endpoint quotients are rounded outward directly, so a representable
    /// quotient does not require a representable reciprocal. Subnormal results
    /// retain their enclosure, including when ordinary division rounds to zero.
    /// Repeated operands are treated independently: `x / x` need not be `[1, 1]`.
    ///
    /// # Examples
    /// ```
    /// use la_stack::prelude::*;
    /// # fn main() -> Result<(), LaError> {
    /// let tiny = Interval::point(f64::from_bits(1))?;
    /// assert_eq!(tiny.try_div(&tiny)?, Interval::ONE);
    /// assert_eq!(Interval::try_new(-2.0, 6.0)?.try_div(&Interval::try_new(2.0, 4.0)?)?,
    ///            Interval::try_new(-1.0, 3.0)?);
    /// # Ok(())
    /// # }
    /// ```
    ///
    /// # Errors
    /// Returns [`LaError::IntervalZeroDivisor`] if `other` contains zero
    /// (even for a zero numerator), or [`LaError::IntervalRangeExhausted`] if
    /// no finite enclosure exists.
    #[inline]
    pub const fn try_div(&self, other: &Self) -> Result<Self, LaError> {
        self.try_div_for(other, ArithmeticOperation::IntervalDivision)
    }

    /// Negate an interval exactly by swapping and negating its endpoints.
    ///
    /// # Examples
    /// ```
    /// use la_stack::prelude::*;
    ///
    /// # fn main() -> Result<(), LaError> {
    /// let range = Interval::try_new(-2.0, 3.0)?;
    /// assert_eq!(range.negate(), Interval::try_new(-3.0, 2.0)?);
    /// # Ok(())
    /// # }
    /// ```
    #[inline]
    pub const fn negate(&self) -> Self {
        Self::new_unchecked(-self.upper, -self.lower)
    }

    /// Square an interval with outward rounding.
    ///
    /// An interval spanning zero has exact lower bound zero. The upper bound is
    /// the outward-rounded square of the endpoint with greatest magnitude.
    ///
    /// # Examples
    /// ```
    /// use la_stack::prelude::*;
    ///
    /// # fn main() -> Result<(), LaError> {
    /// let range = Interval::try_new(-2.0, 3.0)?;
    /// assert_eq!(range.try_square()?, Interval::try_new(0.0, 9.0)?);
    /// // Multiplication treats its two operands independently and is wider.
    /// assert_eq!(range.try_mul(&range)?, Interval::try_new(-6.0, 9.0)?);
    /// # Ok(())
    /// # }
    /// ```
    ///
    /// # Errors
    /// Returns [`LaError::IntervalRangeExhausted`] when the exact square range
    /// has no finite binary64 enclosure.
    #[expect(
        clippy::inline_always,
        reason = "measured call overhead in lifted interval assembly"
    )]
    #[inline(always)]
    pub const fn try_square(&self) -> Result<Self, LaError> {
        let operation = ArithmeticOperation::IntervalSquare;
        let left_square = match rounded_product_bounds(self.lower, self.lower, operation) {
            Ok(bounds) => bounds,
            Err(error) => return Err(error),
        };
        if self.is_point() {
            return Ok(Self::new_unchecked(left_square.0, left_square.1));
        }
        let right_square = match rounded_product_bounds(self.upper, self.upper, operation) {
            Ok(bounds) => bounds,
            Err(error) => return Err(error),
        };
        let lower = if self.lower <= 0.0 && self.upper >= 0.0 {
            0.0
        } else if left_square.0 < right_square.0 {
            left_square.0
        } else {
            right_square.0
        };
        let upper = if left_square.1 > right_square.1 {
            left_square.1
        } else {
            right_square.1
        };
        Ok(Self::new_unchecked(lower, upper))
    }

    /// Construct an interval after its finite ordered-bound invariant is known.
    #[inline]
    const fn new_unchecked(lower: f64, upper: f64) -> Self {
        Self {
            lower: canonical_zero(lower),
            upper: canonical_zero(upper),
        }
    }

    /// Add while attributing range failure to the owning public operation.
    #[expect(
        clippy::inline_always,
        reason = "measured result-copy overhead in the interval determinant kernel"
    )]
    #[inline(always)]
    const fn try_add_for(
        &self,
        other: &Self,
        operation: ArithmeticOperation,
    ) -> Result<Self, LaError> {
        if self.is_zero() {
            return Ok(*other);
        }
        if other.is_zero() {
            return Ok(*self);
        }
        if self.is_point() && other.is_point() {
            return match rounded_add_bounds(self.lower, other.lower, operation) {
                Ok((lower, upper)) => Ok(Self::new_unchecked(lower, upper)),
                Err(error) => Err(error),
            };
        }

        let lower = match rounded_add_endpoint::<false>(self.lower, other.lower, operation) {
            Ok(lower) => lower,
            Err(error) => return Err(error),
        };
        let upper = match rounded_add_endpoint::<true>(self.upper, other.upper, operation) {
            Ok(upper) => upper,
            Err(error) => return Err(error),
        };
        Ok(Self::new_unchecked(lower, upper))
    }

    /// Multiply while attributing range failure to the owning public operation.
    #[expect(
        clippy::inline_always,
        reason = "measured result-copy overhead in the interval determinant kernel"
    )]
    #[inline(always)]
    const fn try_mul_for(
        &self,
        other: &Self,
        operation: ArithmeticOperation,
    ) -> Result<Self, LaError> {
        if self.is_zero() || other.is_zero() {
            return Ok(Self::ZERO);
        }
        if self.is_one() {
            return Ok(*other);
        }
        if other.is_one() {
            return Ok(*self);
        }
        if self.is_point() && other.is_point() {
            return match rounded_product_bounds(self.lower, other.lower, operation) {
                Ok((lower, upper)) => Ok(Self::new_unchecked(lower, upper)),
                Err(error) => Err(error),
            };
        }

        self.try_mul_by_sign(other, operation)
    }

    /// Select only the endpoint products that can attain each range extremum.
    #[inline]
    const fn try_mul_by_sign(
        &self,
        other: &Self,
        operation: ArithmeticOperation,
    ) -> Result<Self, LaError> {
        let self_nonnegative = self.lower >= 0.0;
        let self_nonpositive = self.upper <= 0.0;
        let other_nonnegative = other.lower >= 0.0;
        let other_nonpositive = other.upper <= 0.0;

        if self_nonnegative {
            if other_nonnegative {
                return Self::try_product_extrema(
                    (self.lower, other.lower),
                    (self.upper, other.upper),
                    operation,
                );
            }
            if other_nonpositive {
                return Self::try_product_extrema(
                    (self.upper, other.lower),
                    (self.lower, other.upper),
                    operation,
                );
            }
            return Self::try_product_extrema(
                (self.upper, other.lower),
                (self.upper, other.upper),
                operation,
            );
        }
        if self_nonpositive {
            if other_nonnegative {
                return Self::try_product_extrema(
                    (self.lower, other.upper),
                    (self.upper, other.lower),
                    operation,
                );
            }
            if other_nonpositive {
                return Self::try_product_extrema(
                    (self.upper, other.upper),
                    (self.lower, other.lower),
                    operation,
                );
            }
            return Self::try_product_extrema(
                (self.lower, other.upper),
                (self.lower, other.lower),
                operation,
            );
        }
        if other_nonnegative {
            return Self::try_product_extrema(
                (self.lower, other.upper),
                (self.upper, other.upper),
                operation,
            );
        }
        if other_nonpositive {
            return Self::try_product_extrema(
                (self.upper, other.lower),
                (self.lower, other.lower),
                operation,
            );
        }

        let lower_left = match rounded_product_bounds(self.lower, other.upper, operation) {
            Ok(bounds) => bounds,
            Err(error) => return Err(error),
        };
        let lower_right = match rounded_product_bounds(self.upper, other.lower, operation) {
            Ok(bounds) => bounds,
            Err(error) => return Err(error),
        };
        let upper_left = match rounded_product_bounds(self.lower, other.lower, operation) {
            Ok(bounds) => bounds,
            Err(error) => return Err(error),
        };
        let upper_right = match rounded_product_bounds(self.upper, other.upper, operation) {
            Ok(bounds) => bounds,
            Err(error) => return Err(error),
        };
        let lower = if lower_left.0 < lower_right.0 {
            lower_left.0
        } else {
            lower_right.0
        };
        let upper = if upper_left.1 > upper_right.1 {
            upper_left.1
        } else {
            upper_right.1
        };
        Ok(Self::new_unchecked(lower, upper))
    }

    /// Enclose the selected exact lower and upper product extrema.
    #[inline]
    const fn try_product_extrema(
        lower_factors: (f64, f64),
        upper_factors: (f64, f64),
        operation: ArithmeticOperation,
    ) -> Result<Self, LaError> {
        let lower = match rounded_product_bounds(lower_factors.0, lower_factors.1, operation) {
            Ok((lower, _)) => lower,
            Err(error) => return Err(error),
        };
        let upper = match rounded_product_bounds(upper_factors.0, upper_factors.1, operation) {
            Ok((_, upper)) => upper,
            Err(error) => return Err(error),
        };
        Ok(Self::new_unchecked(lower, upper))
    }

    /// Divide while attributing range failure to the owning public operation.
    const fn try_div_for(
        &self,
        other: &Self,
        operation: ArithmeticOperation,
    ) -> Result<Self, LaError> {
        if other.contains(0.0) {
            return Err(LaError::IntervalZeroDivisor {
                lower: other.lower,
                upper: other.upper,
            });
        }
        // On a rectangle with zero-free denominator, extrema occur at corners.
        // Point operands avoid computing the same quotient four times.
        let numerators = [self.lower, self.upper];
        let denominators = [other.lower, other.upper];
        let numerator_count = if self.is_point() { 1 } else { 2 };
        let denominator_count = if other.is_point() { 1 } else { 2 };
        let mut lower = f64::MAX;
        let mut upper = -f64::MAX;
        let mut i = 0;
        while i < numerator_count {
            let mut j = 0;
            while j < denominator_count {
                let (lo, hi) =
                    match rounded_quotient_bounds(numerators[i], denominators[j], operation) {
                        Ok(bounds) => bounds,
                        Err(error) => return Err(error),
                    };
                if lo < lower {
                    lower = lo;
                }
                if hi > upper {
                    upper = hi;
                }
                j += 1;
            }
            i += 1;
        }
        Ok(Self::new_unchecked(lower, upper))
    }

    /// Enclose `self - factor * value` with the solve's error attribution.
    const fn try_eliminate(&self, factor: &Self, value: &Self) -> Result<Self, LaError> {
        let product = match factor.try_mul_for(value, ArithmeticOperation::IntervalSolve) {
            Ok(product) => product,
            Err(error) => return Err(error),
        };
        self.try_add_for(&product.negate(), ArithmeticOperation::IntervalSolve)
    }

    /// Lower bound on the magnitude of every represented value.
    const fn separation_from_zero(&self) -> f64 {
        if self.lower > 0.0 {
            self.lower
        } else if self.upper < 0.0 {
            -self.upper
        } else {
            0.0
        }
    }

    /// Return whether this interval is exactly real zero.
    #[inline]
    const fn is_zero(&self) -> bool {
        self.lower == 0.0 && self.upper == 0.0
    }

    /// Return whether this interval is exactly real one.
    #[inline]
    const fn is_one(&self) -> bool {
        self.lower.to_bits() == 1.0_f64.to_bits() && self.upper.to_bits() == 1.0_f64.to_bits()
    }

    /// Return whether this interval contains one binary64 point.
    #[inline]
    const fn is_point(&self) -> bool {
        self.lower.to_bits() == self.upper.to_bits()
    }
}

impl Default for Interval {
    #[inline]
    fn default() -> Self {
        Self::ZERO
    }
}

impl<const D: usize> IntervalMatrix<D> {
    /// Construct an interval matrix from already-validated interval rows.
    ///
    /// See [`det`](Self::det) for an example with non-point entries.
    #[inline]
    pub const fn from_rows(rows: [[Interval; D]; D]) -> Self {
        Self { rows }
    }

    /// Lift finite binary64 rows into point intervals.
    ///
    /// This preserves the stored binary64 values exactly; it does not recover
    /// uncertainty from arithmetic performed before this call.
    /// See [`IntervalMatrix`] for a determinant-sign example using this constructor.
    ///
    /// # Errors
    /// Returns [`LaError::NonFinite`] with matrix coordinates for the first NaN
    /// or infinity in row-major order.
    #[inline]
    pub const fn try_from_point_rows(rows: [[f64; D]; D]) -> Result<Self, LaError> {
        let mut intervals = [[Interval::ZERO; D]; D];
        let mut row = 0;
        while row < D {
            let mut column = 0;
            while column < D {
                let value = rows[row][column];
                if !value.is_finite() {
                    return Err(LaError::non_finite_input_matrix(row, column));
                }
                intervals[row][column] = Interval::new_unchecked(value, value);
                column += 1;
            }
            row += 1;
        }
        Ok(Self::from_rows(intervals))
    }

    /// Lift a finite [`Matrix`] into point intervals.
    ///
    /// Earlier rounded expression construction is not enclosed; use interval
    /// operations while constructing derived coefficients when that uncertainty
    /// belongs in the proof.
    ///
    /// # Examples
    /// ```
    /// use la_stack::prelude::*;
    ///
    /// # fn main() -> Result<(), LaError> {
    /// let matrix = Matrix::<2>::try_from_rows([[2.0, 0.0], [0.0, 3.0]])?;
    /// let intervals = IntervalMatrix::from_matrix(&matrix);
    /// assert_eq!(intervals.det()?, Interval::point(6.0)?);
    /// # Ok(())
    /// # }
    /// ```
    #[inline]
    pub const fn from_matrix(matrix: &Matrix<D>) -> Self {
        let matrix_rows = matrix.as_rows();
        let mut intervals = [[Interval::ZERO; D]; D];
        let mut row = 0;
        while row < D {
            let mut column = 0;
            while column < D {
                let value = matrix_rows[row][column];
                intervals[row][column] = Interval::new_unchecked(value, value);
                column += 1;
            }
            row += 1;
        }
        Self::from_rows(intervals)
    }

    /// All-zero interval matrix.
    #[inline]
    pub const fn zero() -> Self {
        Self::from_rows([[Interval::ZERO; D]; D])
    }

    /// Identity interval matrix.
    #[inline]
    pub const fn identity() -> Self {
        let mut matrix = Self::zero();
        let mut index = 0;
        while index < D {
            matrix.rows[index][index] = Interval::ONE;
            index += 1;
        }
        matrix
    }

    /// Borrow the row-major interval storage.
    #[inline]
    #[must_use]
    pub const fn as_rows(&self) -> &[[Interval; D]; D] {
        &self.rows
    }

    /// Consume this matrix and return its row-major interval storage.
    #[inline]
    #[must_use]
    pub const fn into_rows(self) -> [[Interval; D]; D] {
        self.rows
    }

    /// Get an interval entry with bounds checking.
    #[inline]
    #[must_use]
    pub const fn get(&self, row: usize, column: usize) -> Option<Interval> {
        if row < D && column < D {
            Some(self.rows[row][column])
        } else {
            None
        }
    }

    /// Get an interval entry while preserving index context on failure.
    ///
    /// See [`set`](Self::set) for an example of mutation and checked access.
    ///
    /// # Errors
    /// Returns [`LaError::IndexOutOfBounds`] when either index is not `< D`.
    #[inline]
    pub const fn try_get(&self, row: usize, column: usize) -> Result<Interval, LaError> {
        if row < D && column < D {
            Ok(self.rows[row][column])
        } else {
            Err(LaError::index_out_of_bounds(row, column, D))
        }
    }

    /// Set an interval entry with bounds checking.
    ///
    /// Validation is unnecessary for the value because [`Interval`] already
    /// carries the finite ordered-bound proof. An invalid index leaves the
    /// matrix unchanged.
    ///
    /// # Examples
    /// ```
    /// use core::assert_matches;
    /// use la_stack::prelude::*;
    ///
    /// # fn main() -> Result<(), LaError> {
    /// let mut matrix = IntervalMatrix::<2>::identity();
    /// let range = Interval::try_new(2.0, 3.0)?;
    /// matrix.set(0, 0, range)?;
    /// assert_eq!(matrix.try_get(0, 0)?, range);
    /// let before = matrix;
    /// assert_matches!(
    ///     matrix.set(2, 0, Interval::ZERO),
    ///     Err(LaError::IndexOutOfBounds { row: 2, col: 0, dim: 2, .. })
    /// );
    /// assert_eq!(matrix, before);
    /// # Ok(())
    /// # }
    /// ```
    ///
    /// # Errors
    /// Returns [`LaError::IndexOutOfBounds`] when either index is not `< D`.
    #[inline]
    pub const fn set(&mut self, row: usize, column: usize, value: Interval) -> Result<(), LaError> {
        if row >= D || column >= D {
            return Err(LaError::index_out_of_bounds(row, column, D));
        }
        self.rows[row][column] = value;
        Ok(())
    }

    /// Enclose the solution of every real system `A x = b` admitted by these
    /// coefficient intervals and `rhs`.
    ///
    /// On success every admitted `A` is nonsingular and the unique exact-real
    /// solution satisfies `x[i] ∈ result[i]`. Entries vary independently; any
    /// correlations introduced during expression assembly are safely widened.
    /// The enclosure need not be tight. Wide intervals and ill-conditioning
    /// can prevent certification even when every admitted matrix is invertible.
    ///
    /// Uses outward-rounded Gaussian elimination and back substitution, with
    /// row pivots maximizing their minimum distance from zero (first row wins
    /// ties). No tolerance admits a zero-containing pivot. See `REFERENCES.md`
    /// \[20\] and the
    /// [enclosure proof](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#interval-linear-solves).
    ///
    /// Both shapes are fixed by `D`. All const dimensions are supported, with
    /// O(D³) work and O(D²) inline scratch space intended for small systems;
    /// D=0 returns the empty solution. Inputs are unchanged on all outcomes.
    /// This operation is available without the `exact` feature.
    ///
    /// # Examples
    /// ```
    /// use la_stack::prelude::*;
    /// # fn main() -> Result<(), LaError> {
    /// let matrix = IntervalMatrix::try_from_point_rows([[0.0, 2.0], [3.0, 1.0]])?;
    /// let rhs = [Interval::point(4.0)?, Interval::point(5.0)?];
    /// let x = matrix.solve(&rhs)?;
    /// assert!(x[0].contains(1.0) && x[1].contains(2.0));
    /// # Ok(())
    /// # }
    /// ```
    ///
    /// # Errors
    /// Returns [`LaError::IntervalSolveInconclusive`] at the first column whose
    /// remaining pivot intervals all contain zero, including the final column.
    /// This never establishes singularity, even for point inputs or `[0, 0]`
    /// pivots; this method never returns [`LaError::Singular`].
    /// Returns [`LaError::IntervalRangeExhausted`] with
    /// [`ArithmeticOperation::IntervalSolve`] if an intermediate or result has
    /// no finite enclosure. Both outcomes require the caller's fallback when a
    /// certified answer is needed. Invalid raw endpoints are rejected by the
    /// interval constructors before this method can be called.
    pub const fn solve(&self, rhs: &[Interval; D]) -> Result<[Interval; D], LaError> {
        let mut rows = self.rows;
        let mut work_rhs = *rhs;
        let mut column = 0;
        while column < D {
            let mut pivot_row = column;
            let mut separation = 0.0;
            let mut row = column;
            while row < D {
                let candidate = rows[row][column].separation_from_zero();
                if candidate > separation {
                    separation = candidate;
                    pivot_row = row;
                }
                row += 1;
            }
            if separation == 0.0 {
                return Err(LaError::IntervalSolveInconclusive { pivot_col: column });
            }
            rows.swap(column, pivot_row);
            work_rhs.swap(column, pivot_row);
            row = column + 1;
            while row < D {
                if !rows[row][column].is_zero() {
                    let factor = match rows[row][column]
                        .try_div_for(&rows[column][column], ArithmeticOperation::IntervalSolve)
                    {
                        Ok(factor) => factor,
                        Err(error) => return Err(error),
                    };
                    let mut j = column + 1;
                    while j < D {
                        rows[row][j] = match rows[row][j].try_eliminate(&factor, &rows[column][j]) {
                            Ok(value) => value,
                            Err(error) => return Err(error),
                        };
                        j += 1;
                    }
                    work_rhs[row] = match work_rhs[row].try_eliminate(&factor, &work_rhs[column]) {
                        Ok(value) => value,
                        Err(error) => return Err(error),
                    };
                    // The exact row operation cancels this entry identically
                    // for every admitted system, despite interval dependency.
                    rows[row][column] = Interval::ZERO;
                }
                row += 1;
            }
            column += 1;
        }
        let mut solution = [Interval::ZERO; D];
        let mut row = D;
        while row > 0 {
            row -= 1;
            let mut residual = work_rhs[row];
            let mut column = row + 1;
            while column < D {
                residual = match residual.try_eliminate(&rows[row][column], &solution[column]) {
                    Ok(value) => value,
                    Err(error) => return Err(error),
                };
                column += 1;
            }
            solution[row] =
                match residual.try_div_for(&rows[row][row], ArithmeticOperation::IntervalSolve) {
                    Ok(value) => value,
                    Err(error) => return Err(error),
                };
        }
        Ok(solution)
    }

    /// Enclose the determinant with division-free subset dynamic programming.
    ///
    /// For each column subset, the DP stores the determinant interval of the
    /// leading rows and those columns. This evaluates the Leibniz expansion in
    /// `D × 2^(D-1)` products and additions without choosing or dividing by a
    /// pivot. The returned interval therefore encloses every exact-real
    /// determinant represented by the input intervals, subject only to an
    /// explicit range failure.
    ///
    /// The D=0 determinant follows the empty-product convention and is `[1, 1]`.
    ///
    /// Use [`det_sign`](Self::det_sign) when only sign evidence is needed.
    ///
    /// # Examples
    /// ```
    /// use la_stack::prelude::*;
    ///
    /// # fn main() -> Result<(), LaError> {
    /// // Every represented diagonal matrix has a determinant in [8, 15].
    /// let matrix = IntervalMatrix::<2>::from_rows([
    ///     [Interval::try_new(2.0, 3.0)?, Interval::ZERO],
    ///     [Interval::ZERO, Interval::try_new(4.0, 5.0)?],
    /// ]);
    /// assert_eq!(matrix.det()?, Interval::try_new(8.0, 15.0)?);
    /// # Ok(())
    /// # }
    /// ```
    ///
    /// # Errors
    /// Returns [`LaError::UnsupportedDimension`] for D>7. Returns
    /// [`LaError::IntervalRangeExhausted`] with interval-determinant provenance
    /// when an exact intermediate has no finite binary64 enclosure; callers can
    /// then proceed to an exact or higher-range fallback.
    #[inline]
    pub const fn det(&self) -> Result<Interval, LaError> {
        // Stable Rust cannot express [Interval; 1 << D]. Dispatch only the
        // workspace extent; every arm uses the identical arithmetic tree.
        match D {
            0 | 1 => self.det_with_workspace::<2>(),
            2 => self.det_with_workspace::<4>(),
            3 => self.det_with_workspace::<8>(),
            4 => self.det_with_workspace::<16>(),
            5 => self.det_with_workspace::<32>(),
            6 => self.det_with_workspace::<64>(),
            7 => self.det_with_workspace::<128>(),
            _ => Err(LaError::unsupported_dimension(D, MAX_INTERVAL_MATRIX_DIM)),
        }
    }

    /// Evaluate with the state count established by the dimension dispatch.
    ///
    /// Only `det` calls this helper, after proving D is supported and STATES
    /// covers all `2^D` subsets. Its D=0 and D=1 arms reserve two slots. These
    /// preconditions keep shifts and indexing in range while preserving the
    /// same expansion order for every workspace size.
    #[inline]
    const fn det_with_workspace<const STATES: usize>(&self) -> Result<Interval, LaError> {
        let state_count = 1_usize << D;
        let mut partials = [Interval::ZERO; STATES];
        partials[0] = Interval::ONE;
        let operation = ArithmeticOperation::IntervalDeterminant;

        let mut subset = 1;
        while subset < state_count {
            let row = subset.count_ones() as usize - 1;
            let mut sum = Interval::ZERO;
            // Removing the lowest set bit visits active columns in the same
            // ascending order as the full scan, preserving every interval
            // operation and the first range-exhaustion error.
            let mut remaining = subset;
            while remaining != 0 {
                let column = remaining.trailing_zeros() as usize;
                let column_bit = 1_usize << column;
                remaining &= remaining - 1;
                let previous = subset ^ column_bit;
                let mut term =
                    match partials[previous].try_mul_for(&self.rows[row][column], operation) {
                        Ok(term) => term,
                        Err(error) => return Err(error),
                    };
                if !remaining.count_ones().is_multiple_of(2) {
                    term = term.negate();
                }
                sum = match sum.try_add_for(&term, operation) {
                    Ok(next_sum) => next_sum,
                    Err(error) => return Err(error),
                };
            }
            partials[subset] = sum;
            subset += 1;
        }

        Ok(partials[state_count - 1])
    }

    /// Return proof-bearing determinant sign evidence.
    ///
    /// An interval strictly on one side of zero proves that sign. Only the
    /// singleton interval `[0, 0]` proves `Zero`; every other overlap with zero
    /// is [`IntervalDeterminantSign::Inconclusive`].
    ///
    /// # Examples
    /// ```
    /// use la_stack::prelude::*;
    ///
    /// # fn main() -> Result<(), LaError> {
    /// let mut matrix = IntervalMatrix::<2>::identity();
    /// assert_eq!(matrix.det_sign()?, IntervalDeterminantSign::Positive);
    ///
    /// matrix.set(0, 0, Interval::try_new(-1.0, 1.0)?)?;
    /// // This range includes nonsingular matrices of both signs; a caller
    /// // needs tighter or exact input before it can decide singularity.
    /// assert_eq!(matrix.det_sign()?, IntervalDeterminantSign::Inconclusive);
    ///
    /// matrix.set(0, 0, Interval::ZERO)?;
    /// assert_eq!(matrix.det_sign()?, IntervalDeterminantSign::Zero);
    /// # Ok(())
    /// # }
    /// ```
    ///
    /// # Errors
    /// Propagates the dimension and arithmetic range failures from
    /// [`det`](Self::det).
    #[inline]
    pub const fn det_sign(&self) -> Result<IntervalDeterminantSign, LaError> {
        let determinant = match self.det() {
            Ok(determinant) => determinant,
            Err(error) => return Err(error),
        };
        if determinant.lower > 0.0 {
            Ok(IntervalDeterminantSign::Positive)
        } else if determinant.upper < 0.0 {
            Ok(IntervalDeterminantSign::Negative)
        } else if determinant.lower == 0.0 && determinant.upper == 0.0 {
            Ok(IntervalDeterminantSign::Zero)
        } else {
            Ok(IntervalDeterminantSign::Inconclusive)
        }
    }
}

impl<const D: usize> Default for IntervalMatrix<D> {
    #[inline]
    fn default() -> Self {
        Self::zero()
    }
}

#[cfg(test)]
mod tests {
    use core::assert_matches;

    use pastey::paste;

    use super::*;
    use crate::{IntervalBound, IntervalOperand, NonFiniteLocation, NonFiniteOrigin};

    #[test]
    fn point_and_bounds_enforce_interval_invariants() {
        assert_eq!(Interval::point(-0.0).unwrap().lower().to_bits(), 0);
        assert_eq!(Interval::try_new(-0.0, 0.0).unwrap(), Interval::ZERO);
        assert_matches!(
            Interval::point(f64::NAN),
            Err(LaError::NonFinite {
                location: NonFiniteLocation::Scalar,
                origin: NonFiniteOrigin::Input,
                ..
            })
        );
        assert_matches!(
            Interval::try_new(2.0, 1.0),
            Err(LaError::InvertedInterval {
                lower: 2.0,
                upper: 1.0,
                ..
            })
        );
    }

    #[test]
    fn constructors_preserve_non_finite_input_locations() {
        for value in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
            assert_eq!(
                Interval::try_new(value, 0.0),
                Err(LaError::NonFinite {
                    location: NonFiniteLocation::IntervalBound {
                        bound: IntervalBound::Lower,
                    },
                    origin: NonFiniteOrigin::Input,
                })
            );
            assert_eq!(
                Interval::try_new(0.0, value),
                Err(LaError::NonFinite {
                    location: NonFiniteLocation::IntervalBound {
                        bound: IntervalBound::Upper,
                    },
                    origin: NonFiniteOrigin::Input,
                })
            );
            assert_eq!(
                Interval::try_from_subtraction(value, 0.0),
                Err(LaError::NonFinite {
                    location: NonFiniteLocation::IntervalOperand {
                        operand: IntervalOperand::Left,
                    },
                    origin: NonFiniteOrigin::Input,
                })
            );
            assert_eq!(
                Interval::try_from_subtraction(0.0, value),
                Err(LaError::NonFinite {
                    location: NonFiniteLocation::IntervalOperand {
                        operand: IntervalOperand::Right,
                    },
                    origin: NonFiniteOrigin::Input,
                })
            );
        }

        assert_eq!(
            Interval::try_new(f64::NAN, f64::INFINITY),
            Err(LaError::NonFinite {
                location: NonFiniteLocation::IntervalBound {
                    bound: IntervalBound::Lower,
                },
                origin: NonFiniteOrigin::Input,
            })
        );
        assert_eq!(
            Interval::try_from_subtraction(f64::NAN, f64::INFINITY),
            Err(LaError::NonFinite {
                location: NonFiniteLocation::IntervalOperand {
                    operand: IntervalOperand::Left,
                },
                origin: NonFiniteOrigin::Input,
            })
        );

        let rows = [[0.0, f64::NAN], [f64::INFINITY, 0.0]];
        assert_eq!(
            IntervalMatrix::<2>::try_from_point_rows(rows),
            Err(LaError::NonFinite {
                location: NonFiniteLocation::MatrixCell { row: 0, col: 1 },
                origin: NonFiniteOrigin::Input,
            })
        );
    }

    #[test]
    fn exact_operations_remain_point_intervals() -> Result<(), LaError> {
        let one = Interval::point(1.0)?;
        let two = Interval::point(2.0)?;
        assert_eq!(one.try_add(&two)?, Interval::point(3.0)?);
        assert_eq!(two.try_mul(&two)?, Interval::point(4.0)?);
        assert_eq!(Interval::try_from_subtraction(3.0, 2.0)?, one);
        assert_eq!(
            Interval::try_new(-2.0, -1.0)?.negate(),
            Interval::try_new(1.0, 2.0)?
        );
        Ok(())
    }

    #[test]
    fn inexact_operations_expand_only_in_the_required_direction() -> Result<(), LaError> {
        let subtraction = Interval::try_from_subtraction(1.0, 0.1)?;
        let rounded_subtraction = 1.0_f64 - 0.1;
        assert_eq!(
            subtraction,
            Interval::try_new(rounded_subtraction.next_down(), rounded_subtraction)?
        );

        let product = Interval::point(0.1)?.try_mul(&Interval::point(0.2)?)?;
        let rounded_product = 0.1_f64 * 0.2;
        assert_eq!(
            product,
            Interval::try_new(rounded_product.next_down(), rounded_product)?
        );

        let below_one = 1.0 - f64::EPSILON;
        let above_one = 1.0 + f64::EPSILON;
        let binade_boundary = Interval::point(below_one)?.try_mul(&Interval::point(above_one)?)?;
        assert_eq!(
            binade_boundary,
            Interval::try_new(1.0_f64.next_down(), 1.0)?
        );
        Ok(())
    }

    #[test]
    fn cancellation_preserves_an_exact_ulp_difference() -> Result<(), LaError> {
        let next = 1.0_f64.next_up();
        let difference = Interval::try_from_subtraction(next, 1.0)?;
        assert_eq!(difference, Interval::point(f64::EPSILON)?);
        Ok(())
    }

    #[test]
    fn subtraction_by_or_from_zero_preserves_exact_endpoint_bits() -> Result<(), LaError> {
        for value in [
            -f64::MAX,
            -1.0,
            -f64::MIN_POSITIVE,
            -f64::from_bits(1),
            -0.0,
            0.0,
            f64::from_bits(1),
            f64::MIN_POSITIVE,
            1.0,
            f64::MAX,
        ] {
            for zero in [0.0, -0.0] {
                for (left, right, expected) in [(value, zero, value), (zero, value, -value)] {
                    let result = Interval::try_from_subtraction(left, right)?;
                    let expected_bits = if expected == 0.0 {
                        0
                    } else {
                        expected.to_bits()
                    };
                    assert_eq!(result.lower().to_bits(), expected_bits);
                    assert_eq!(result.upper().to_bits(), expected_bits);
                }
            }
        }
        Ok(())
    }

    #[test]
    fn singleton_squares_keep_tight_bounds_including_underflow() -> Result<(), LaError> {
        for (value, lower, upper) in [
            (0.0, 0.0_f64, 0.0),
            (-0.0, 0.0, 0.0),
            (2.0, 4.0, 4.0),
            (-2.0, 4.0, 4.0),
            (f64::from_bits(1), 0.0, f64::from_bits(1)),
            (-f64::from_bits(1), 0.0, f64::from_bits(1)),
        ] {
            let squared = Interval::point(value)?.try_square()?;
            assert_eq!(squared.lower().to_bits(), lower.to_bits());
            assert_eq!(squared.upper().to_bits(), upper.to_bits());
        }
        // Squaring the exact stored binary64 0.1 lies strictly between these
        // adjacent values; the expected bits come from its rational square.
        for value in [0.1, -0.1] {
            let squared = Interval::point(value)?.try_square()?;
            assert_eq!(squared.lower().to_bits(), 0x3f84_7ae1_47ae_147b);
            assert_eq!(squared.upper().to_bits(), 0x3f84_7ae1_47ae_147c);
        }
        Ok(())
    }

    #[test]
    fn same_sign_squares_select_the_nearest_and_farthest_endpoints() -> Result<(), LaError> {
        let expected = Interval::try_new(4.0, 9.0)?;
        for (lower, upper) in [(-3.0, -2.0), (2.0, 3.0)] {
            assert_eq!(Interval::try_new(lower, upper)?.try_square()?, expected);
        }
        Ok(())
    }

    #[test]
    fn wide_addition_rounds_only_the_required_endpoints() -> Result<(), LaError> {
        let left = Interval::try_new(-1.0, 1.0)?;
        // A quarter ULP rounds both sums inward; three quarters rounds them
        // outward already. Both exact dyadic sums have the same tight bounds.
        for radius in [f64::EPSILON / 4.0, 3.0 * f64::EPSILON / 4.0] {
            let right = Interval::try_new(-radius, radius)?;
            assert_eq!(
                left.try_add(&right)?,
                Interval::try_new((-1.0_f64).next_down(), 1.0_f64.next_up())?
            );
        }
        assert_eq!(
            left.try_add(&Interval::try_new(-2.0, 2.0)?)?,
            Interval::try_new(-3.0, 3.0)?
        );
        Ok(())
    }

    #[test]
    fn underflowed_product_still_encloses_the_positive_exact_result() -> Result<(), LaError> {
        let least_subnormal = f64::from_bits(1);
        let product = Interval::point(least_subnormal)?.try_mul(&Interval::point(0.5)?)?;
        assert_eq!(product, Interval::try_new(0.0, least_subnormal)?);
        Ok(())
    }

    #[test]
    fn range_failure_preserves_interval_operation() -> Result<(), LaError> {
        let error = Interval::point(f64::MAX)?
            .try_mul(&Interval::point(2.0)?)
            .unwrap_err();
        assert_eq!(
            error,
            LaError::IntervalRangeExhausted {
                operation: ArithmeticOperation::IntervalMultiplication,
            }
        );
        Ok(())
    }

    #[test]
    fn rounded_maximum_detects_exact_sum_beyond_finite_range() -> Result<(), LaError> {
        let maximum = Interval::point(f64::MAX)?;
        let tiny = Interval::point(f64::MIN_POSITIVE)?;
        assert_eq!(
            maximum.try_add(&maximum),
            Err(LaError::IntervalRangeExhausted {
                operation: ArithmeticOperation::IntervalAddition,
            })
        );
        assert_eq!(
            maximum.try_add(&tiny),
            Err(LaError::IntervalRangeExhausted {
                operation: ArithmeticOperation::IntervalAddition,
            })
        );
        let nonnegative = Interval::try_new(0.0, f64::MAX)?;
        assert_eq!(
            nonnegative.try_add(&nonnegative),
            Err(LaError::IntervalRangeExhausted {
                operation: ArithmeticOperation::IntervalAddition,
            })
        );
        let positive = Interval::try_new(f64::from_bits(1), 1.0)?;
        for (left, right) in [(maximum, positive), (maximum.negate(), positive.negate())] {
            assert_eq!(
                left.try_add(&right),
                Err(LaError::IntervalRangeExhausted {
                    operation: ArithmeticOperation::IntervalAddition,
                }),
            );
        }

        let finite_difference = maximum.try_add(&tiny.negate())?;
        assert_eq!(finite_difference.upper().to_bits(), f64::MAX.to_bits());
        assert!(finite_difference.lower() < finite_difference.upper());
        Ok(())
    }

    #[test]
    fn subtraction_and_square_preserve_distinct_range_operations() -> Result<(), LaError> {
        assert_eq!(
            Interval::try_from_subtraction(f64::MAX, -f64::MIN_POSITIVE),
            Err(LaError::IntervalRangeExhausted {
                operation: ArithmeticOperation::IntervalSubtraction,
            })
        );
        assert_eq!(
            Interval::point(f64::MAX)?.try_square(),
            Err(LaError::IntervalRangeExhausted {
                operation: ArithmeticOperation::IntervalSquare,
            })
        );
        assert_eq!(
            Interval::try_new(-1.0, f64::MAX)?.try_square(),
            Err(LaError::IntervalRangeExhausted {
                operation: ArithmeticOperation::IntervalSquare,
            })
        );
        Ok(())
    }

    #[test]
    fn determinant_overflow_reports_interval_determinant_range_failure() -> Result<(), LaError> {
        let matrix = IntervalMatrix::<2>::try_from_point_rows([[f64::MAX, 0.0], [0.0, 2.0]])?;
        assert_eq!(
            matrix.det(),
            Err(LaError::IntervalRangeExhausted {
                operation: ArithmeticOperation::IntervalDeterminant,
            })
        );

        let accumulating =
            IntervalMatrix::<2>::try_from_point_rows([[f64::MAX, f64::MAX], [-1.0, 1.0]])?;
        assert_eq!(
            accumulating.det(),
            Err(LaError::IntervalRangeExhausted {
                operation: ArithmeticOperation::IntervalDeterminant,
            })
        );
        assert_eq!(
            accumulating.det_sign(),
            Err(LaError::IntervalRangeExhausted {
                operation: ArithmeticOperation::IntervalDeterminant,
            })
        );
        Ok(())
    }

    #[test]
    fn determinant_reports_intermediate_exhaustion_before_exact_cancellation() -> Result<(), LaError>
    {
        let matrix = IntervalMatrix::<2>::try_from_point_rows([[f64::MAX, f64::MAX], [2.0, 2.0]])?;
        assert_eq!(
            matrix.det(),
            Err(LaError::IntervalRangeExhausted {
                operation: ArithmeticOperation::IntervalDeterminant,
            })
        );
        Ok(())
    }

    #[test]
    fn square_spanning_zero_has_exact_zero_lower_bound() -> Result<(), LaError> {
        let square = Interval::try_new(-2.0, 3.0)?.try_square()?;
        assert_eq!(square, Interval::try_new(0.0, 9.0)?);
        Ok(())
    }

    #[test]
    fn multiplication_selects_correct_extrema_in_every_sign_quadrant() -> Result<(), LaError> {
        for (left, right, expected) in [
            ((2.0, 3.0), (4.0, 5.0), (8.0, 15.0)),
            ((2.0, 3.0), (-5.0, -4.0), (-15.0, -8.0)),
            ((2.0, 3.0), (-5.0, 4.0), (-15.0, 12.0)),
            ((-3.0, -2.0), (4.0, 5.0), (-15.0, -8.0)),
            ((-3.0, -2.0), (-5.0, -4.0), (8.0, 15.0)),
            ((-3.0, -2.0), (-5.0, 4.0), (-12.0, 15.0)),
            ((-3.0, 2.0), (4.0, 5.0), (-15.0, 10.0)),
            ((-3.0, 2.0), (-5.0, -4.0), (-10.0, 15.0)),
            ((-3.0, 2.0), (-5.0, 4.0), (-12.0, 15.0)),
        ] {
            let product = Interval::try_new(left.0, left.1)?
                .try_mul(&Interval::try_new(right.0, right.1)?)?;
            assert_eq!(product, Interval::try_new(expected.0, expected.1)?);
        }

        assert_eq!(
            Interval::ZERO.try_mul(&Interval::try_new(-f64::MAX, f64::MAX)?)?,
            Interval::ZERO
        );
        Ok(())
    }

    #[test]
    fn multiplication_rejects_unrepresentable_selected_extrema() -> Result<(), LaError> {
        let half_maximum = f64::MAX / 2.0;
        for (left, right) in [
            ((half_maximum, f64::MAX), (-2.0, -1.0)),
            ((half_maximum, f64::MAX), (1.0, 2.0)),
            ((-f64::MAX, 1.0), (-1.0, 2.0)),
            ((-1.0, f64::MAX), (-2.0, 1.0)),
            ((-f64::MAX, 1.0), (-2.0, 1.0)),
            ((-1.0, f64::MAX), (-1.0, 2.0)),
        ] {
            let result =
                Interval::try_new(left.0, left.1)?.try_mul(&Interval::try_new(right.0, right.1)?);
            assert_eq!(
                result,
                Err(LaError::IntervalRangeExhausted {
                    operation: ArithmeticOperation::IntervalMultiplication,
                }),
                "left={left:?}, right={right:?}"
            );
        }
        Ok(())
    }

    macro_rules! gen_interval_identity_tests {
        ($d:literal) => {
            paste! {
                #[test]
                fn [<interval_identity_sign_is_positive_ $d d>]() {
                    let matrix = IntervalMatrix::<$d>::identity();
                    assert_eq!(matrix.det(), Ok(Interval::ONE));
                    assert_eq!(
                        matrix.det_sign(),
                        Ok(IntervalDeterminantSign::Positive)
                    );
                }
            }
        };
    }

    gen_interval_identity_tests!(0);
    gen_interval_identity_tests!(1);
    gen_interval_identity_tests!(2);
    gen_interval_identity_tests!(3);
    gen_interval_identity_tests!(4);
    gen_interval_identity_tests!(5);
    gen_interval_identity_tests!(6);
    gen_interval_identity_tests!(7);

    #[test]
    fn smallest_determinants_preserve_const_evaluation_and_interval_signs() -> Result<(), LaError> {
        const EMPTY: Result<Interval, LaError> = IntervalMatrix::<0>::zero().det();
        const WIDE: Interval = match Interval::try_new(-f64::MAX, f64::MAX) {
            Ok(interval) => interval,
            Err(_) => panic!("finite ordered literal bounds"),
        };
        const SINGLE: Result<Interval, LaError> = IntervalMatrix::from_rows([[WIDE]]).det();
        assert_eq!(EMPTY, Ok(Interval::ONE));
        assert_eq!(SINGLE, Ok(WIDE));
        for (lower, upper, sign) in [
            (
                -f64::MAX,
                -f64::from_bits(1),
                IntervalDeterminantSign::Negative,
            ),
            (
                f64::from_bits(1),
                f64::MAX,
                IntervalDeterminantSign::Positive,
            ),
            (0.0, 0.0, IntervalDeterminantSign::Zero),
            (-f64::MAX, f64::MAX, IntervalDeterminantSign::Inconclusive),
        ] {
            let interval = Interval::try_new(lower, upper)?;
            let matrix = IntervalMatrix::from_rows([[interval]]);
            assert_eq!(matrix.det()?, interval);
            assert_eq!(matrix.det_sign()?, sign);
        }
        Ok(())
    }

    #[test]
    fn determinant_sign_handles_row_swap_and_exact_singularity() -> Result<(), LaError> {
        let swapped = IntervalMatrix::<3>::try_from_point_rows([
            [0.0, 1.0, 0.0],
            [1.0, 0.0, 0.0],
            [0.0, 0.0, 1.0],
        ])?;
        assert_eq!(swapped.det_sign()?, IntervalDeterminantSign::Negative);

        let singular = IntervalMatrix::<3>::try_from_point_rows([
            [1.0, 2.0, 3.0],
            [1.0, 2.0, 3.0],
            [0.0, 0.0, 1.0],
        ])?;
        assert_eq!(singular.det_sign()?, IntervalDeterminantSign::Zero);
        Ok(())
    }

    #[test]
    fn wide_determinant_interval_is_inconclusive() -> Result<(), LaError> {
        let matrix = IntervalMatrix::<2>::from_rows([
            [Interval::ONE, Interval::ZERO],
            [Interval::ZERO, Interval::try_new(-1.0, 1.0)?],
        ]);
        assert_eq!(matrix.det_sign()?, IntervalDeterminantSign::Inconclusive);
        Ok(())
    }

    #[test]
    fn determinant_rejects_dimensions_above_supported_stack_dp() {
        assert_matches!(
            IntervalMatrix::<8>::identity().det(),
            Err(LaError::UnsupportedDimension {
                requested: 8,
                max: MAX_INTERVAL_MATRIX_DIM,
                ..
            })
        );
    }

    #[test]
    fn matrix_accessors_preserve_validated_storage() -> Result<(), LaError> {
        let source = Matrix::<2>::identity();
        let mut intervals = IntervalMatrix::from_matrix(&source);
        let value = Interval::try_new(2.0, 3.0)?;
        intervals.set(0, 1, value)?;
        assert_eq!(intervals.get(0, 1), Some(value));
        assert_eq!(intervals.get(2, 0), None);
        assert_eq!(intervals.try_get(0, 1)?, value);
        assert_matches!(
            intervals.try_get(2, 0),
            Err(LaError::IndexOutOfBounds {
                row: 2,
                col: 0,
                dim: 2,
                ..
            })
        );
        assert_eq!(intervals.as_rows()[0][1], value);
        assert_eq!(intervals.into_rows()[0][1], value);
        Ok(())
    }

    #[test]
    fn rejected_matrix_set_is_failure_atomic() -> Result<(), LaError> {
        let mut matrix = IntervalMatrix::<2>::identity();
        let before = matrix;
        let value = Interval::try_new(2.0, 3.0)?;

        assert_eq!(
            matrix.set(2, 0, value),
            Err(LaError::IndexOutOfBounds {
                row: 2,
                col: 0,
                dim: 2,
            })
        );
        assert_eq!(matrix, before);
        assert_eq!(
            matrix.set(0, 2, value),
            Err(LaError::IndexOutOfBounds {
                row: 0,
                col: 2,
                dim: 2,
            })
        );
        assert_eq!(matrix, before);
        Ok(())
    }
}