cadmpeg-ir 0.6.2

// SPDX-License-Identifier: Apache-2.0
//! Finite sums with an exact-product fallback for floating-point range loss.

use crate::scalar::{FiniteReal, NonZeroReal};

/// A finite dot product with an exact-product fallback for range loss or cancellation.
///
/// A non-finite input, and a sum outside the finite range, leave the dot
/// product non-finite; it then carries the plain left-to-right sum of the
/// products.
pub(crate) fn finite_dot<const N: usize>(
    coefficients: [f64; N],
    components: [f64; N],
) -> Result<FiniteReal, f64> {
    let products = std::array::from_fn(|index| coefficients[index] * components[index]);
    let plain = || {
        products
            .into_iter()
            .reduce(|sum, product| sum + product)
            .unwrap_or(0.0)
    };
    if coefficients
        .iter()
        .chain(&components)
        .any(|value| !value.is_finite())
    {
        return Err(plain());
    }
    if let Some(value) = fast_dot(coefficients, components, products) {
        return Ok(value);
    }
    let mut sum = ExactSignedSum::default();
    for (coefficient, component) in coefficients.into_iter().zip(components) {
        sum.add_product(coefficient, component);
    }
    sum.finite_sum().ok_or_else(plain)
}

// The smallest exponent a finite `f64` significand carries. Every exponent
// this module handles is biased by it, which keeps the biased value unsigned
// and makes the product shift below a `usize` with no conversion.
const MIN_SIGNIFICAND_EXPONENT: i32 = -1074;

// Four finite factors need at most 212 significand bits. The accumulator
// includes their entire exponent range, 64 low guard bits for rounded scaled
// terms, and enough high carry bits for any addressable collection of terms.
const EXACT_PRODUCT_EXPONENT: i32 = 4 * MIN_SIGNIFICAND_EXPONENT - 64;
const EXACT_SUM_WORDS: usize = 138;
const EXACT_SUM_WORDS_I32: i32 = 138;
const SUBNORMAL_UNIT_BIT: usize = 3286;

#[derive(Clone, Copy)]
pub(crate) struct ExactSignedSum {
    positive: [u64; EXACT_SUM_WORDS],
    negative: [u64; EXACT_SUM_WORDS],
    window: Window,
    invalid: bool,
}

impl Default for ExactSignedSum {
    fn default() -> Self {
        Self {
            positive: [0; EXACT_SUM_WORDS],
            negative: [0; EXACT_SUM_WORDS],
            window: Window::EMPTY,
            invalid: false,
        }
    }
}

/// The accumulator words `low..high`. Every word outside them is zero on both
/// sides, so a scan of the accumulator reads only these.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
struct Window {
    low: usize,
    high: usize,
}

impl Window {
    const EMPTY: Self = Self {
        low: EXACT_SUM_WORDS,
        high: 0,
    };

    /// Widen the window over the words `low..high`.
    fn include(&mut self, low: usize, high: usize) {
        self.low = self.low.min(low);
        self.high = self.high.max(high);
    }

    fn words(self) -> std::ops::Range<usize> {
        self.low..self.high.max(self.low)
    }
}

/// The largest biased significand exponent a finite `f64` states. The field is
/// eleven bits wide and a finite value never reaches the all-ones field, and
/// biasing by [`MIN_SIGNIFICAND_EXPONENT`] subtracts one more.
const MAX_BIASED_SIGNIFICAND_EXPONENT: i32 = 0x7fe - 1;

/// The widest significand a finite `f64` states, the hidden bit included.
///
/// [`finite_significand`] proves it. Its subnormal branch answers a nonzero
/// `fraction`, which is below `1 << 52`, and its normal branch answers
/// `(1 << 52) | fraction`. The highest set bit is therefore at most 52 and the
/// width at most 53.
const MAX_SIGNIFICAND_BITS: i32 = 53;

/// The least power of two a [`ScaledExponent`] holds. `ExactSignedSum::finish`
/// states `EXACT_PRODUCT_EXPONENT + highest_bit + 1` with `highest_bit` at
/// least zero, which is below everything `scaled_finite` can state.
const MIN_SCALED_EXPONENT: i32 = EXACT_PRODUCT_EXPONENT + 1;

/// The greatest power of two a [`ScaledExponent`] holds.
/// `ExactSignedSum::finish` states `EXACT_PRODUCT_EXPONENT + highest_bit + 2`
/// in its rounding-carry arm, with `highest_bit` at most one below the
/// accumulator's bit width.
const MAX_SCALED_EXPONENT: i32 = EXACT_PRODUCT_EXPONENT + EXACT_SUM_WORDS_I32 * 64 + 1;

/// The power of two that scales a [`ScaledValue`]'s mantissa back to the value.
///
/// The field is private and the type has no constructor of its own, so the two
/// construction expressions in [`ScaledValue::of_nonzero`] and
/// [`ExactSignedSum::finish`] are the only values that exist. Neither reads an
/// exponent from outside, and the four `const` assertions below state that
/// both reach only `MIN_SCALED_EXPONENT..=MAX_SCALED_EXPONENT`.
#[derive(Clone, Copy, PartialEq, Eq, PartialOrd, Ord)]
pub(crate) struct ScaledExponent(i32);

/// `ScaledValue::of_nonzero` states the exponent of the value's own
/// significand, `MIN_SIGNIFICAND_EXPONENT + biased + bits`, with `biased` from
/// zero to `MAX_BIASED_SIGNIFICAND_EXPONENT` and `bits` from one to
/// `MAX_SIGNIFICAND_BITS`: the raise of a subnormal value is taken back off.
const _: () = assert!(MIN_SIGNIFICAND_EXPONENT + 1 >= MIN_SCALED_EXPONENT);
const _: () = assert!(
    MIN_SIGNIFICAND_EXPONENT + MAX_BIASED_SIGNIFICAND_EXPONENT + MAX_SIGNIFICAND_BITS
        <= MAX_SCALED_EXPONENT
);

/// `ExactSignedSum::finish` states `EXACT_PRODUCT_EXPONENT + highest_bit + 1`,
/// and `+ 2` in its rounding-carry arm, with `highest_bit` indexing an
/// accumulator `EXACT_SUM_WORDS` words of 64 bits wide.
const _: () = assert!(EXACT_PRODUCT_EXPONENT + 1 >= MIN_SCALED_EXPONENT);
const _: () =
    assert!(EXACT_PRODUCT_EXPONENT + (EXACT_SUM_WORDS_I32 * 64 - 1) + 2 <= MAX_SCALED_EXPONENT);

impl ScaledExponent {
    /// The exponent shift that rescales a value in this frame into `other`'s.
    ///
    /// Plain `-`: both sides lie in `MIN_SCALED_EXPONENT..=MAX_SCALED_EXPONENT`,
    /// so the difference is at most `MAX_SCALED_EXPONENT - MIN_SCALED_EXPONENT`
    /// in magnitude.
    fn difference(self, other: Self) -> i32 {
        self.0 - other.0
    }
}

#[derive(Clone, Copy)]
pub(crate) struct ScaledValue {
    sign: f64,
    mantissa: f64,
    pub(crate) exponent: ScaledExponent,
}

impl ScaledValue {
    /// Twice a positive interval width whose half-width is finite. The
    /// exponent is at most one above a finite binary64 exponent.
    pub(crate) fn doubled(self) -> Self {
        Self {
            exponent: ScaledExponent(self.exponent.0 + 1),
            ..self
        }
    }

    pub(crate) fn exponent(self) -> i32 {
        self.exponent.0
    }

    pub(crate) fn rescale(self, exponent: i32) -> Option<FiniteReal> {
        super::scale_power_of_two(
            self.sign * self.mantissa,
            self.exponent.0.checked_sub(exponent)?,
        )
    }

    /// The value, rescaled into the frame whose exponent is `scale_exponent`.
    pub(crate) fn scaled_by(self, scale_exponent: ScaledExponent) -> f64 {
        super::scale_power_of_two(
            self.sign * self.mantissa,
            self.exponent.difference(scale_exponent),
        )
        .map_or(self.sign * f64::INFINITY, FiniteReal::get)
    }
    /// The value, or the signed infinity it overflows to. The mantissa is
    /// finite and nonzero, so only the scaling can leave the finite range.
    pub(crate) fn finite(self) -> Result<FiniteReal, f64> {
        super::scale_power_of_two(self.sign * self.mantissa, self.exponent.0)
            .ok_or(self.sign * f64::INFINITY)
    }

    /// `self / denominator`, or the signed infinity the quotient overflows
    /// to. Both mantissas are in `[0.5, 1)`, so their ratio is finite and
    /// nonzero, and the exponent difference lies in the `ScaledExponent` span:
    /// only the final scaling can leave the finite range.
    pub(crate) fn quotient(self, denominator: Self) -> Result<FiniteReal, f64> {
        let ratio = self.sign * denominator.sign * (self.mantissa / denominator.mantissa);
        super::scale_power_of_two(ratio, self.exponent.difference(denominator.exponent))
            .ok_or(ratio.signum() * f64::INFINITY)
    }

    /// Divide by up to five nonzero extended-range factors before rounding to
    /// binary64, doubling the quotient when `doubled` is set. Each normalized
    /// mantissa is in `[0.5, 1)`, so the small product of mantissa quotients
    /// stays in range independently of exponent. A quotient that overflows
    /// carries its signed infinity.
    ///
    /// Plain `+` and `-`: every exponent lies in
    /// `MIN_SCALED_EXPONENT..=MAX_SCALED_EXPONENT`, and six of them and a
    /// doubling stay far inside `i32`.
    pub(crate) fn quotient_by_factors<const N: usize>(
        self,
        denominators: [Self; N],
        doubled: bool,
    ) -> Result<FiniteReal, f64> {
        const { assert!(N <= 5) };
        let mut mantissa = self.sign * self.mantissa;
        let mut exponent = self.exponent.0 + i32::from(doubled);
        for denominator in denominators {
            mantissa /= denominator.sign * denominator.mantissa;
            exponent -= denominator.exponent.0;
        }
        super::scale_power_of_two(mantissa, exponent).ok_or(mantissa.signum() * f64::INFINITY)
    }

    /// `self / denominator` times `2^exponent_shift`, without rounding either
    /// operand into the binary64 range first: the ratio of the two mantissas
    /// lies in `(0.5, 2)`, and only the final scaling rounds. The caller
    /// states the bound that keeps the product finite.
    ///
    /// Plain `+`: the exponent difference lies in the `ScaledExponent` span,
    /// and a shift within the binary64 exponent range keeps the sum far inside
    /// `i32`.
    pub(crate) fn quotient_shifted_product(self, denominator: Self, exponent_shift: i32) -> f64 {
        super::power_of_two_product(
            self.sign * denominator.sign * (self.mantissa / denominator.mantissa),
            self.exponent.difference(denominator.exponent) + exponent_shift,
        )
    }

    /// The scaled form of a nonzero finite value: its sign, the mantissa of
    /// its magnitude in `[0.5, 1)`, and the power of two between them.
    ///
    /// A subnormal magnitude is first raised by `2^RAISE`, which is exact and
    /// makes it normal, and the raise is taken back off the exponent. A normal
    /// magnitude's significand is `(1 << 52) | fraction`, whose highest bit is
    /// set, so its 53 bits form the mantissa without a search. Every nonzero
    /// finite value has this form, so nothing is checked.
    pub(crate) fn of_nonzero(value: NonZeroReal) -> Self {
        // The least subnormal magnitude is `2^-1074`; raised, it is `2^-1010`,
        // above the least normal magnitude `2^-1022`.
        const RAISE: i32 = 64;
        let value = value.get();
        let (normal, raised) = if value.abs() < f64::MIN_POSITIVE {
            (value * 2.0_f64.powi(RAISE), RAISE)
        } else {
            (value, 0)
        };
        let bits = normal.to_bits();
        // The mask keeps eleven bits, so the field is an `i32` by
        // construction; a normal field is at least one.
        let field_bits = ((bits >> 52) & 0x7ff).to_le_bytes();
        // endian-exception: reconstructed-scalar
        let field = i32::from(u16::from_le_bytes([field_bits[0], field_bits[1]]));
        let significand = (1_u64 << 52) | (bits & ((1_u64 << 52) - 1));
        Self {
            sign: if bits >> 63 != 0 { -1.0 } else { 1.0 },
            // The hidden bit and fraction encode the exact normalized mantissa.
            mantissa: f64::from_bits((1022_u64 << 52) | (significand & ((1_u64 << 52) - 1))),
            exponent: ScaledExponent(
                field - 1 + MIN_SIGNIFICAND_EXPONENT + MAX_SIGNIFICAND_BITS - raised,
            ),
        }
    }

    /// The scaled form of the exact product of nonzero finite factors,
    /// rounded once as [`ExactSignedSum::finish`] rounds the same product.
    ///
    /// The accumulator holds the one product on the side its sign names, and
    /// every factor's significand is nonzero, so that side is the product's
    /// magnitude and has a highest set bit. Nothing is checked.
    pub(crate) fn product_of_nonzero<const N: usize>(factors: [NonZeroReal; N]) -> Option<Self> {
        let mut sum = ExactSignedSum::default();
        sum.add_factors(factors.map(NonZeroReal::get));
        if sum.invalid {
            return None;
        }
        let negative = factors
            .iter()
            .fold(false, |negative, factor| negative ^ (factor.get() < 0.0));
        let magnitude = if negative {
            &sum.negative
        } else {
            &sum.positive
        };
        rounded_magnitude(
            negative,
            magnitude,
            sum.window.low,
            highest_bit(magnitude, sum.window)?,
        )
    }
}

/// The index of the highest set bit of a nonzero magnitude whose nonzero
/// words lie in `window`.
///
/// `word` indexes `EXACT_SUM_WORDS` words of 64 bits, so a bit index of the
/// accumulator is a `u16`. `checked_ilog2` answers `None` for a zero word,
/// which `rposition` has already excluded; the `?` states that rather than
/// asserting it.
fn highest_bit(magnitude: &[u64; EXACT_SUM_WORDS], window: Window) -> Option<u16> {
    let words = window.words();
    let word = words.start + magnitude[words].iter().rposition(|value| *value != 0)?;
    u16::try_from(
        word * 64 + cadmpeg_core::decode::index_from_u32(magnitude[word].checked_ilog2()?),
    )
    .ok()
}

/// The sign, the integer significand, and the significand's exponent biased by
/// [`MIN_SIGNIFICAND_EXPONENT`], of a finite `f64`.
///
/// The biased exponent is a `u16`: the field is eleven bits wide, and a finite
/// value never reaches the all-ones field, so the biased value is at most 2045.
fn finite_significand(value: f64) -> Option<(bool, u64, u16)> {
    if !value.is_finite() {
        return None;
    }
    let bits = value.to_bits();
    let negative = bits >> 63 != 0;
    // The mask keeps eleven bits, so the field is a `u16` by construction.
    let exponent_field = u16::try_from((bits >> 52) & 0x7ff).ok()?;
    let fraction = bits & ((1_u64 << 52) - 1);
    if exponent_field == 0 {
        (fraction != 0).then_some((negative, fraction, 0))
    } else {
        // A normal value's exponent is `field - 1075`, which is
        // `field - 1` once biased by `MIN_SIGNIFICAND_EXPONENT`.
        Some((negative, (1_u64 << 52) | fraction, exponent_field - 1))
    }
}

/// Add `value` at word `index`, carrying upwards, and widen `window` over every
/// word written.
fn add_word(words: &mut [u64; EXACT_SUM_WORDS], window: &mut Window, index: usize, value: u64) {
    let low = index;
    let (sum, mut carry) = words[index].overflowing_add(value);
    words[index] = sum;
    let mut index = index + 1;
    while carry {
        let (sum, next) = words[index].overflowing_add(1);
        words[index] = sum;
        carry = next;
        index += 1;
    }
    window.include(low, index);
}

fn add_shifted(words: &mut [u64; EXACT_SUM_WORDS], window: &mut Window, value: u128, shift: usize) {
    let bytes = value.to_le_bytes();
    // endian-exception: reconstructed-scalar
    let low = u64::from_le_bytes([
        bytes[0], bytes[1], bytes[2], bytes[3], bytes[4], bytes[5], bytes[6], bytes[7],
    ]);
    // endian-exception: reconstructed-scalar
    let high = u64::from_le_bytes([
        bytes[8], bytes[9], bytes[10], bytes[11], bytes[12], bytes[13], bytes[14], bytes[15],
    ]);
    for (limb_index, limb) in [low, high].into_iter().enumerate() {
        if limb == 0 {
            continue;
        }
        let bit = shift + limb_index * 64;
        let word = bit / 64;
        let remainder = bit % 64;
        if remainder == 0 {
            add_word(words, window, word, limb);
        } else {
            add_word(words, window, word, limb << remainder);
            add_word(words, window, word + 1, limb >> (64 - remainder));
        }
    }
}

impl ExactSignedSum {
    pub(crate) fn add_product(&mut self, left: f64, right: f64) {
        self.add_factors([left, right]);
    }

    pub(crate) fn add_factors<const N: usize>(&mut self, factors: [f64; N]) {
        // All callers supply at most four factors, the tensor-product point case.
        assert!(N <= 4);
        let mut product = [1_u64, 0, 0, 0];
        let mut negative = false;
        let mut exponent = 0_i32;
        for factor in factors {
            let Some((sign, significand, biased)) = finite_significand(factor) else {
                return;
            };
            negative ^= sign;
            exponent += MIN_SIGNIFICAND_EXPONENT + i32::from(biased);
            let mut carry = 0_u128;
            for word in &mut product {
                let value = u128::from(*word) * u128::from(significand) + carry;
                let bytes = value.to_le_bytes();
                // endian-exception: reconstructed-scalar
                *word = u64::from_le_bytes([
                    bytes[0], bytes[1], bytes[2], bytes[3], bytes[4], bytes[5], bytes[6], bytes[7],
                ]);
                carry = value >> 64;
            }
        }
        let Ok(shift) = usize::try_from(exponent - EXACT_PRODUCT_EXPONENT) else {
            self.invalid = true;
            return;
        };
        let target = if negative {
            &mut self.negative
        } else {
            &mut self.positive
        };
        add_shifted(
            target,
            &mut self.window,
            u128::from(product[0]) | (u128::from(product[1]) << 64),
            shift,
        );
        add_shifted(
            target,
            &mut self.window,
            u128::from(product[2]) | (u128::from(product[3]) << 64),
            shift + 128,
        );
    }

    /// Add a rounded extended-range sum times a finite coefficient. Rational
    /// quotient derivatives use this to subtract weight derivatives before division.
    pub(crate) fn add_scaled_product(
        &mut self,
        value: Option<ScaledValue>,
        factor: FiniteReal,
    ) -> Option<()> {
        let Some(value) = value else {
            return Some(());
        };
        let Some((negative, significand, exponent)) =
            finite_significand(value.sign * value.mantissa)
        else {
            return Some(());
        };
        let Some((factor_negative, factor_significand, factor_exponent)) =
            finite_significand(factor.get())
        else {
            return Some(());
        };
        let product = u128::from(significand) * u128::from(factor_significand);
        let shift = value.exponent.0
            + 2 * MIN_SIGNIFICAND_EXPONENT
            + i32::from(exponent)
            + i32::from(factor_exponent)
            - EXACT_PRODUCT_EXPONENT;
        let shift = usize::try_from(shift).ok()?;
        if shift + 128 >= EXACT_SUM_WORDS * 64 {
            return None;
        }
        let target = if negative ^ factor_negative {
            &mut self.negative
        } else {
            &mut self.positive
        };
        add_shifted(target, &mut self.window, product, shift);
        Some(())
    }

    pub(crate) fn finish(mut self) -> Option<ScaledValue> {
        if self.invalid {
            return None;
        }
        let window = self.window;
        let (negative, magnitude) = self.signed_difference()?;
        rounded_magnitude(
            negative,
            magnitude,
            window.low,
            highest_bit(magnitude, window)?,
        )
    }

    /// Round a sum in the subnormal range directly onto the binary64 grid.
    /// Rounding through a 53-bit scaled mantissa first can discard a tail that
    /// breaks a tie at half of the smallest subnormal.
    pub(crate) fn finite_sum(mut self) -> Option<FiniteReal> {
        if self.invalid {
            return None;
        }
        let window = self.window;
        let Some((negative, magnitude)) = self.signed_difference() else {
            return Some(FiniteReal::ZERO);
        };
        let highest = highest_bit(magnitude, window)?;
        let highest_bit = usize::from(highest);
        if highest_bit >= SUBNORMAL_UNIT_BIT + 52 {
            return rounded_magnitude(negative, magnitude, window.low, highest)?
                .finite()
                .ok();
        }
        let mut units = 0_u64;
        if highest_bit >= SUBNORMAL_UNIT_BIT {
            units = bit_field(
                magnitude,
                SUBNORMAL_UNIT_BIT,
                highest_bit + 1 - SUBNORMAL_UNIT_BIT,
            );
        }
        let guard = bit_is_set(magnitude, SUBNORMAL_UNIT_BIT - 1);
        let sticky = any_bit_below(magnitude, window.low, SUBNORMAL_UNIT_BIT - 1);
        if guard && (sticky || units & 1 != 0) {
            units += 1;
        }
        if units == 0 {
            return Some(FiniteReal::ZERO);
        }
        FiniteReal::subnormal_units(negative, units)
    }

    /// Returns the sign and the magnitude of `positive - negative`, or `None`
    /// when the two accumulators are equal. The magnitude is written over the
    /// larger side, and its nonzero words lie in the window.
    ///
    /// The comparison that names the larger accumulator is the only route to
    /// the subtraction below, so the minuend is never the smaller side and the
    /// loop cannot leave a borrow in the highest word. Splitting the
    /// comparison from the subtraction would state that relation as an
    /// assertion instead of holding it.
    fn signed_difference(&mut self) -> Option<(bool, &[u64; EXACT_SUM_WORDS])> {
        let words = self.window.words();
        let is_negative = self.positive[words.clone()]
            .iter()
            .zip(&self.negative[words.clone()])
            .rev()
            .find_map(|(left, right)| match left.cmp(right) {
                std::cmp::Ordering::Greater => Some(false),
                std::cmp::Ordering::Less => Some(true),
                std::cmp::Ordering::Equal => None,
            })?;
        let (larger, smaller) = if is_negative {
            (&mut self.negative, &self.positive)
        } else {
            (&mut self.positive, &self.negative)
        };
        let mut borrow = false;
        for index in words {
            let (difference, first_borrow) = larger[index].overflowing_sub(smaller[index]);
            let (difference, second_borrow) = difference.overflowing_sub(u64::from(borrow));
            larger[index] = difference;
            borrow = first_borrow || second_borrow;
        }
        Some((is_negative, larger))
    }
}

fn bit_is_set(words: &[u64; EXACT_SUM_WORDS], bit: usize) -> bool {
    words[bit / 64] & (1_u64 << (bit % 64)) != 0
}

/// The `count` bits of `words` from bit `low` upwards, `count` at most 64,
/// with bit `low` as the least significant.
fn bit_field(words: &[u64; EXACT_SUM_WORDS], low: usize, count: usize) -> u64 {
    let (word, offset) = (low / 64, low % 64);
    let mut field = words[word] >> offset;
    if offset != 0 && word + 1 < EXACT_SUM_WORDS {
        field |= words[word + 1] << (64 - offset);
    }
    if count < 64 {
        field & ((1_u64 << count) - 1)
    } else {
        field
    }
}

/// Whether any bit of `words` below `bit` is set, when every word below word
/// `low` is zero.
fn any_bit_below(words: &[u64; EXACT_SUM_WORDS], low: usize, bit: usize) -> bool {
    let (word, offset) = (bit / 64, bit % 64);
    words[low.min(word)..word].iter().any(|value| *value != 0)
        || words[word] & ((1_u64 << offset) - 1) != 0
}

/// The scaled form of an accumulator magnitude whose highest set bit is
/// `highest_bit`, rounded to 53 bits, half to even, with `negative` as its
/// sign. Every word below word `low` is zero.
fn rounded_magnitude(
    negative: bool,
    magnitude: &[u64; EXACT_SUM_WORDS],
    low: usize,
    highest_bit: u16,
) -> Option<ScaledValue> {
    let keep = (highest_bit + 1).min(53);
    let mut significand = bit_field(
        magnitude,
        usize::from(highest_bit + 1 - keep),
        usize::from(keep),
    );
    if keep == 53 {
        let guard_bit = highest_bit
            .checked_sub(keep)
            .is_some_and(|bit| bit_is_set(magnitude, usize::from(bit)));
        let sticky = highest_bit
            .checked_sub(keep)
            .is_some_and(|bit| any_bit_below(magnitude, low, usize::from(bit)));
        if guard_bit && (sticky || significand & 1 != 0) {
            significand += 1;
            if significand == 1_u64 << 53 {
                return Some(ScaledValue {
                    sign: if negative { -1.0 } else { 1.0 },
                    mantissa: 0.5,
                    exponent: ScaledExponent(EXACT_PRODUCT_EXPONENT + i32::from(highest_bit) + 2),
                });
            }
        }
    }
    Some(ScaledValue {
        sign: if negative { -1.0 } else { 1.0 },
        mantissa: cadmpeg_core::convert::f64_from_u64(significand)?
            * super::normal_power_of_two(-i32::from(keep))?,
        exponent: ScaledExponent(EXACT_PRODUCT_EXPONENT + i32::from(highest_bit) + 1),
    })
}

/// The scaled form of a finite nonzero value; a zero or non-finite value has
/// none.
pub(crate) fn scaled_finite(value: f64) -> Option<ScaledValue> {
    NonZeroReal::new(value).map(ScaledValue::of_nonzero)
}

/// `value / denominator`, multiplied by each factor through one quotient that
/// keeps the extended exponent range. A knot span below the normal range makes
/// the plain quotient overflow although every product here is finite, and one
/// quotient states the same rounding for every factor.
///
/// A zero value, and a zero denominator that a repeated knot states, make every
/// product zero. `None` states a product that no finite `f64` holds.
pub(crate) fn scaled_ratio_products<const N: usize>(
    value: FiniteReal,
    denominator: FiniteReal,
    factors: [FiniteReal; N],
) -> Option<[FiniteReal; N]> {
    // A finite value has a scaled form exactly when it is not zero.
    let (Some(value), Some(denominator)) =
        (scaled_finite(value.get()), scaled_finite(denominator.get()))
    else {
        return Some([FiniteReal::ZERO; N]);
    };
    // Both significands are in `[0.5, 1)`, so the ratio is in `(0.5, 2)` and
    // every product below stays normal until `scale_power_of_two` states it.
    let ratio = value.sign * denominator.sign * (value.mantissa / denominator.mantissa);
    // Plain `+`: the difference lies inside the `ScaledExponent` span and a
    // factor's exponent inside the `f64` range.
    let exponent = value.exponent.difference(denominator.exponent);
    let mut products = [FiniteReal::ZERO; N];
    for (product, factor) in products.iter_mut().zip(factors) {
        let Some(factor) = scaled_finite(factor.get()) else {
            continue;
        };
        *product = super::scale_power_of_two(
            ratio * factor.sign * factor.mantissa,
            exponent + factor.exponent.0,
        )?;
    }
    Some(products)
}

pub(crate) fn fast_dot<const N: usize>(
    coefficients: [f64; N],
    components: [f64; N],
    products: [f64; N],
) -> Option<FiniteReal> {
    if products.iter().any(|product| !product.is_finite()) {
        return None;
    }
    if coefficients.into_iter().zip(components).zip(products).any(
        |((coefficient, component), product)| {
            // A product below the normal range may have rounded before it
            // reached this check, even when it is nonzero. Recompute every
            // such column from the exact significands.
            coefficient != 0.0 && component != 0.0 && product.abs() < f64::MIN_POSITIVE
        },
    ) {
        return None;
    }
    let scale = products
        .iter()
        .fold(0.0_f64, |scale, product| scale.max(product.abs()));
    if scale == 0.0 {
        return Some(FiniteReal::ZERO);
    }
    let sum = FiniteReal::new(products.into_iter().sum::<f64>())?;
    if sum.get().abs() / scale < f64::EPSILON {
        return None;
    }
    Some(sum)
}

/// The outcome of [`product_sum`].
pub(crate) enum ProductSum {
    /// A term is absent, or a factor is not finite. The sum has no value.
    Undefined,
    /// The terms cancel to exactly zero.
    Zero,
    /// The sum, in the extended exponent range.
    Value(ScaledValue),
}

impl ProductSum {
    /// The two outcomes `ExactSignedSum::finish` and `scaled_finite` state,
    /// whose `None` is an exact zero.
    fn from_scaled(value: Option<ScaledValue>) -> Self {
        match value {
            Some(value) => Self::Value(value),
            None => Self::Zero,
        }
    }
}

/// Sum products without losing a finite result to an intermediate exponent.
/// Ordinary inputs use f64 arithmetic; range loss and cancellation replay the
/// same terms through the exact accumulator.
pub(crate) fn product_sum<const N: usize>(
    terms: impl Iterator<Item = Option<[f64; N]>> + Clone,
) -> ProductSum {
    let mut sum = 0.0_f64;
    let mut largest = 0.0_f64;
    let mut exact = false;
    for factors in terms.clone() {
        let Some(factors) = factors else {
            return ProductSum::Undefined;
        };
        if factors.iter().any(|value| !value.is_finite()) {
            return ProductSum::Undefined;
        }
        if factors.contains(&0.0) {
            continue;
        }
        let mut product = 1.0;
        for factor in factors {
            product *= factor;
            exact |= !product.is_normal();
        }
        sum += product;
        largest = largest.max(product.abs());
        exact |= !sum.is_finite();
    }
    exact |= largest != 0.0 && sum.abs() / largest < f64::EPSILON;
    if !exact {
        return ProductSum::from_scaled(scaled_finite(sum));
    }
    let mut sum = ExactSignedSum::default();
    for factors in terms {
        let Some(factors) = factors else {
            return ProductSum::Undefined;
        };
        sum.add_factors(factors);
    }
    ProductSum::from_scaled(sum.finish())
}

/// The f64 pass of [`product_sum`] over several sums whose terms share their
/// leading factors. Lane `k` of a term is the product of the shared factors
/// and then `lanes[k]`, formed in that order, so every product, partial sum
/// and range test of a lane is the one `product_sum` forms for the terms
/// `[shared.., lanes[k]]`.
#[derive(Clone, Copy)]
pub(crate) struct LaneProductSums<const LANES: usize> {
    sum: [f64; LANES],
    largest: [f64; LANES],
    exact: [bool; LANES],
}

impl<const LANES: usize> Default for LaneProductSums<LANES> {
    fn default() -> Self {
        Self {
            sum: [0.0; LANES],
            largest: [0.0; LANES],
            exact: [false; LANES],
        }
    }
}

impl<const LANES: usize> LaneProductSums<LANES> {
    /// Add one term. `None` states a factor that is not finite, for which
    /// `product_sum` leaves the sum undefined.
    pub(crate) fn add<const SHARED: usize>(
        &mut self,
        shared: [f64; SHARED],
        lanes: [f64; LANES],
    ) -> Option<()> {
        if shared.iter().chain(&lanes).any(|value| !value.is_finite()) {
            return None;
        }
        if shared.contains(&0.0) {
            return Some(());
        }
        let mut product = 1.0_f64;
        let mut exact = false;
        for factor in shared {
            product *= factor;
            exact |= !product.is_normal();
        }
        for (lane, factor) in lanes.into_iter().enumerate() {
            if factor == 0.0 {
                continue;
            }
            let lane_product = product * factor;
            self.exact[lane] |= exact || !lane_product.is_normal();
            self.sum[lane] += lane_product;
            self.largest[lane] = self.largest[lane].max(lane_product.abs());
            self.exact[lane] |= !self.sum[lane].is_finite();
        }
        Some(())
    }

    /// Each lane's sum where the f64 pass states it, or its replay.
    pub(crate) fn finish(&self) -> [LaneSum; LANES] {
        std::array::from_fn(|lane| {
            let (sum, largest) = (self.sum[lane], self.largest[lane]);
            let exact = self.exact[lane] || (largest != 0.0 && sum.abs() / largest < f64::EPSILON);
            if exact {
                LaneSum::Replay
            } else {
                LaneSum::Settled(scaled_finite(sum))
            }
        })
    }
}

/// A lane of [`LaneProductSums`] after its f64 pass.
#[derive(Clone, Copy)]
pub(crate) enum LaneSum {
    /// The f64 pass states the sum. `None` is an exact zero.
    Settled(Option<ScaledValue>),
    /// A product left the normal range or the sum cancelled: `product_sum`
    /// replays the lane's terms through [`ExactSignedSum`].
    Replay,
}

#[cfg(test)]
mod tests;