la-stack 0.4.6

#![forbid(unsafe_code)]

//! Allocation-free scaled products for floating-point factor diagonals.
//!
//! Binary64 decomposition and rounding follow `REFERENCES.md` \[9-10\].
//! Mantissa normalization and a separate exponent preserve intermediate range;
//! they do not remove rounding in earlier mantissa multiplications. See the
//! [scaled determinant product description](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#scaled-determinant-products)
//! for the replay policy and deferred final-factor rounding.

const SIGN_MASK: u64 = 1_u64 << 63;
const FRACTION_BITS: u32 = 52;
const FRACTION_MASK: u64 = (1_u64 << FRACTION_BITS) - 1;
const EXPONENT_MASK: u64 = 0x7ff;
const EXPONENT_BIAS: i128 = 1023;
const MIN_NORMAL_EXPONENT: i128 = -1022;
const MIN_SUBNORMAL_EXPONENT: i128 = -1074;

/// One direct product step with its range proof attached.
#[derive(Clone, Copy, Debug, PartialEq)]
pub(crate) struct RangeCheckedProduct {
    product: f64,
    range_preserved: bool,
}

impl RangeCheckedProduct {
    /// Return the directly accumulated product.
    #[inline]
    pub(crate) const fn product(self) -> f64 {
        self.product
    }

    /// Return whether direct accumulation preserved the non-zero product's range.
    #[inline]
    pub(crate) const fn range_preserved(self) -> bool {
        self.range_preserved
    }
}

/// Multiply one non-zero direct product step and attach its range proof.
///
/// Successful LU and LDLT construction proves every diagonal factor is finite
/// and non-zero. Starting from `±1.0`, a normal result therefore proves that
/// direct accumulation has not overflowed or lost range through gradual
/// underflow. Callers combine every step's proof and replay the complete product
/// with scaling if any step fails, keeping the success path branch-free between
/// factors.
#[inline]
pub(crate) const fn range_checked_product(accumulator: f64, factor: f64) -> RangeCheckedProduct {
    let product = accumulator * factor;
    let product_exponent = (product.to_bits() >> FRACTION_BITS) & EXPONENT_MASK;
    // Subtracting one maps the valid normal fields 1..=0x7fe to
    // 0..=0x7fd. Zero wraps high and 0x7ff maps to the exclusive upper
    // bound, so the common normal path needs one unsigned comparison.
    RangeCheckedProduct {
        product,
        range_preserved: product_exponent.wrapping_sub(1) < EXPONENT_MASK - 1,
    }
}

/// One non-zero finite factor normalized as `mantissa × 2^exponent`.
#[derive(Clone, Copy)]
struct NormalizedFactor {
    mantissa: f64,
    exponent: i128,
}

/// A finite product kept as `(-1)^negative × mantissa × 2^exponent`.
///
/// Non-zero finite factors are normalized to `1 ≤ mantissa < 2` before
/// multiplication. Consequently, no intermediate mantissa multiplication can
/// underflow or overflow; only [`Self::finish`] decides whether the final
/// rounded result is finite. The most recent factor stays deferred so a final
/// subnormal product can be formed directly in the binary64 destination range
/// without first rounding it as a normal mantissa.
pub(crate) struct ScaledProduct {
    mantissa: f64,
    exponent: i128,
    pending_factor: Option<NormalizedFactor>,
    negative: bool,
    zero: bool,
    non_finite: bool,
}

impl ScaledProduct {
    /// Start an empty product with the requested initial sign.
    #[inline]
    pub(crate) const fn new(negative: bool) -> Self {
        Self {
            mantissa: 1.0,
            exponent: 0,
            pending_factor: None,
            negative,
            zero: false,
            non_finite: false,
        }
    }

    /// Multiply by one factor while retaining a normalized mantissa.
    ///
    /// A non-finite factor is recorded so [`Self::finish`] returns `None`,
    /// even if another factor is zero.
    #[inline]
    pub(crate) const fn multiply(&mut self, factor: f64) {
        let bits = factor.to_bits();
        self.negative ^= bits & SIGN_MASK != 0;

        let magnitude = bits & !SIGN_MASK;
        let biased_exponent = (magnitude >> FRACTION_BITS) & EXPONENT_MASK;
        let fraction = magnitude & FRACTION_MASK;

        if biased_exponent == EXPONENT_MASK {
            self.non_finite = true;
            return;
        }
        if biased_exponent == 0 && fraction == 0 {
            self.zero = true;
            return;
        }
        if self.zero {
            return;
        }

        let (factor_mantissa, factor_exponent) = if biased_exponent == 0 {
            // A subnormal value is `fraction × 2^-1074`. Move its highest set
            // bit to the binary64 hidden-bit position to obtain a mantissa in
            // [1, 2), and compensate in the exponent.
            let highest_bit = fraction.ilog2();
            let shift = FRACTION_BITS - highest_bit;
            let significand = fraction << shift;
            (
                f64::from_bits((1023_u64 << FRACTION_BITS) | (significand & FRACTION_MASK)),
                (highest_bit as i128) + MIN_SUBNORMAL_EXPONENT,
            )
        } else {
            (
                f64::from_bits((1023_u64 << FRACTION_BITS) | fraction),
                (biased_exponent as i128) - EXPONENT_BIAS,
            )
        };

        if let Some(pending) = self.pending_factor {
            self.absorb_factor(pending);
        }
        self.pending_factor = Some(NormalizedFactor {
            mantissa: factor_mantissa,
            exponent: factor_exponent,
        });
    }

    /// Fold one normalized factor into the running product.
    #[inline]
    const fn absorb_factor(&mut self, factor: NormalizedFactor) {
        self.mantissa *= factor.mantissa;
        self.exponent += factor.exponent;

        // Both operands were in [1, 2), and even the two largest binary64
        // mantissas multiply to a value that rounds below 4.0. Therefore one
        // factor-of-two normalization is sufficient to preserve
        // `mantissa < 2`.
        if self.mantissa >= 2.0 {
            self.mantissa *= 0.5;
            self.exponent += 1;
        }
    }

    /// Finalize the accumulated mantissa and pending factor as binary64.
    ///
    /// Returns `None` if any factor was non-finite or the accumulated result
    /// rounds outside the finite binary64 range. Magnitudes below that range
    /// round to a signed zero or subnormal value with round-to-nearest,
    /// ties-to-even semantics.
    /// Earlier mantissa products have already rounded, so this does not
    /// guarantee correct rounding of the exact product of all original factors.
    #[inline]
    #[expect(
        clippy::cast_possible_truncation,
        clippy::cast_sign_loss,
        reason = "the preceding bounds prove the normal and deferred-final biased exponents fit u64"
    )]
    pub(crate) const fn finish(mut self) -> Option<f64> {
        if self.non_finite {
            return None;
        }

        let sign = if self.negative { SIGN_MASK } else { 0 };
        if self.zero {
            return Some(f64::from_bits(sign));
        }
        let Some(pending) = self.pending_factor else {
            return Some(f64::from_bits(sign | (1023_u64 << FRACTION_BITS)));
        };

        let final_exponent = self.exponent + pending.exponent;

        // The product of two mantissas in [1, 2) is strictly below 4. Values
        // below this exponent are therefore strictly below half the least
        // subnormal and round to signed zero.
        if final_exponent < MIN_SUBNORMAL_EXPONENT - 2 {
            return Some(f64::from_bits(sign));
        }

        if final_exponent < MIN_NORMAL_EXPONENT {
            // Scale both normal operands so their single multiplication lands
            // directly in the final subnormal range. Multiplying normalized
            // mantissas first would round once to 53 bits here and a second time
            // to the much coarser subnormal grid below.
            let left = f64::from_bits(
                (1_u64 << FRACTION_BITS) | (self.mantissa.to_bits() & FRACTION_MASK),
            );
            let right_biased_exponent =
                (final_exponent - MIN_NORMAL_EXPONENT + EXPONENT_BIAS) as u64;
            let right = f64::from_bits(
                (right_biased_exponent << FRACTION_BITS)
                    | (pending.mantissa.to_bits() & FRACTION_MASK),
            );
            let magnitude = left * right;
            return Some(f64::from_bits(sign | magnitude.to_bits()));
        }

        self.absorb_factor(pending);
        if self.exponent > EXPONENT_BIAS {
            return None;
        }

        let mantissa_bits = self.mantissa.to_bits();
        let fraction = mantissa_bits & FRACTION_MASK;
        let biased_exponent = (self.exponent + EXPONENT_BIAS) as u64;
        Some(f64::from_bits(
            sign | (biased_exponent << FRACTION_BITS) | fraction,
        ))
    }
}

#[cfg(test)]
mod tests {
    use super::{SIGN_MASK, ScaledProduct, range_checked_product};

    const TWO_NEG_800: f64 = f64::from_bits(223_u64 << 52);
    const TWO_POS_800: f64 = f64::from_bits(1823_u64 << 52);
    const SUBNORMAL_ROUNDING_LEFT: f64 = f64::from_bits(0x3cb2_e219_27ac_435a);
    const SUBNORMAL_ROUNDING_RIGHT: f64 = f64::from_bits(0x0014_55e5_f80b_50eb);

    /// Return the exact bits produced by a two-factor scaled product.
    fn scaled_product_bits(left: f64, right: f64) -> Option<u64> {
        let mut product = ScaledProduct::new(false);
        product.multiply(left);
        product.multiply(right);
        product.finish().map(f64::to_bits)
    }

    #[test]
    fn mantissa_product_is_renormalized_once() {
        assert_eq!(scaled_product_bits(1.5, 1.5), Some(2.25_f64.to_bits()));
    }

    #[test]
    fn signed_zero_tracks_later_factor_signs() {
        let mut product = ScaledProduct::new(false);
        product.multiply(-0.0);
        product.multiply(-2.0);

        assert_eq!(product.finish().map(f64::to_bits), Some(0.0_f64.to_bits()));
    }

    #[test]
    fn empty_product_preserves_initial_sign() {
        assert_eq!(
            ScaledProduct::new(false).finish().map(f64::to_bits),
            Some(1.0_f64.to_bits())
        );
        assert_eq!(
            ScaledProduct::new(true).finish().map(f64::to_bits),
            Some((-1.0_f64).to_bits())
        );
    }

    #[test]
    fn non_finite_factors_make_the_product_unrepresentable() {
        for factor in [f64::INFINITY, f64::NEG_INFINITY, f64::NAN] {
            let mut product = ScaledProduct::new(false);
            product.multiply(factor);

            assert_eq!(product.finish(), None);
        }
    }

    #[test]
    fn non_finite_factors_remain_unrepresentable_before_or_after_zero() {
        for non_finite in [f64::INFINITY, f64::NEG_INFINITY, f64::NAN] {
            for zero in [0.0, -0.0] {
                for middle_factors in [[non_finite, zero], [zero, non_finite]] {
                    let mut product = ScaledProduct::new(false);
                    product.multiply(1.5);
                    for factor in middle_factors {
                        product.multiply(factor);
                    }
                    product.multiply(-2.0);

                    assert_eq!(product.finish(), None);
                }
            }
        }
    }

    #[test]
    fn balanced_extreme_factors_do_not_depend_on_storage_order() {
        let mut forward = ScaledProduct::new(false);
        for factor in [TWO_NEG_800, TWO_NEG_800, TWO_POS_800, TWO_POS_800] {
            forward.multiply(factor);
        }

        let mut reverse = ScaledProduct::new(false);
        for factor in [TWO_POS_800, TWO_POS_800, TWO_NEG_800, TWO_NEG_800] {
            reverse.multiply(factor);
        }

        assert_eq!(forward.finish(), Some(1.0));
        assert_eq!(reverse.finish(), Some(1.0));
    }

    #[test]
    fn final_range_decision_distinguishes_underflow_and_overflow() {
        let mut underflow = ScaledProduct::new(true);
        underflow.multiply(TWO_NEG_800);
        underflow.multiply(TWO_NEG_800);
        assert_eq!(underflow.finish().map(f64::to_bits), Some(1_u64 << 63));

        let mut overflow = ScaledProduct::new(false);
        overflow.multiply(TWO_POS_800);
        overflow.multiply(TWO_POS_800);
        assert_eq!(overflow.finish(), None);
    }

    #[test]
    fn final_subnormal_product_is_rounded_once_in_const_evaluation() {
        const POSITIVE: Option<f64> = {
            let mut product = ScaledProduct::new(false);
            product.multiply(SUBNORMAL_ROUNDING_LEFT);
            product.multiply(SUBNORMAL_ROUNDING_RIGHT);
            product.finish()
        };
        const NEGATIVE: Option<f64> = {
            let mut product = ScaledProduct::new(true);
            product.multiply(SUBNORMAL_ROUNDING_LEFT);
            product.multiply(SUBNORMAL_ROUNDING_RIGHT);
            product.finish()
        };

        assert_eq!(
            (SUBNORMAL_ROUNDING_LEFT * SUBNORMAL_ROUNDING_RIGHT).to_bits(),
            1
        );
        assert_eq!(POSITIVE.map(f64::to_bits), Some(1));
        assert_eq!(NEGATIVE.map(f64::to_bits), Some(SIGN_MASK | 1));
    }

    #[test]
    fn final_subnormal_product_is_rounded_once_after_earlier_range_loss() {
        let factors = [
            TWO_NEG_800,
            TWO_NEG_800,
            TWO_POS_800,
            TWO_POS_800,
            SUBNORMAL_ROUNDING_LEFT,
            SUBNORMAL_ROUNDING_RIGHT,
        ];
        let mut product = ScaledProduct::new(false);
        for factor in factors {
            product.multiply(factor);
        }

        assert!(!range_checked_product(TWO_NEG_800, TWO_NEG_800).range_preserved());
        assert_eq!(product.finish().map(f64::to_bits), Some(1));
    }

    #[test]
    fn final_subnormal_product_preserves_ties_to_even_at_range_boundaries() {
        let least_subnormal = f64::from_bits(1);
        let three_subnormals = f64::from_bits(3);
        let largest_below_one = f64::from_bits(0x3fef_ffff_ffff_ffff);

        assert_eq!(scaled_product_bits(least_subnormal, 0.5), Some(0));
        assert_eq!(scaled_product_bits(three_subnormals, 0.5), Some(2));
        assert_eq!(
            scaled_product_bits(f64::MIN_POSITIVE, largest_below_one),
            Some(f64::MIN_POSITIVE.to_bits())
        );
    }

    #[test]
    fn direct_product_range_proof_distinguishes_ordinary_values_from_range_loss() {
        let ordinary = range_checked_product(1.5, 2.0);
        assert_eq!(ordinary.product().to_bits(), 3.0_f64.to_bits());
        assert!(ordinary.range_preserved());

        for step in [
            range_checked_product(TWO_NEG_800, TWO_NEG_800),
            range_checked_product(f64::MIN_POSITIVE, 0.5),
            range_checked_product(TWO_POS_800, TWO_POS_800),
        ] {
            assert!(!step.range_preserved());
        }
    }
}