g_math 0.6.7

//! Exact wide accumulators shared by the exact predicates, the certified
//! intervals and the fused quadratic form.
//!
//! One trait, [`Wide`], over the profile's wide integer types: signed
//! multiplication on magnitudes with a loud bit-length budget, checked
//! addition, sign test. One alias module, [`acc`], selecting the accumulator
//! per profile from the storage width `W`: `Orient` (orient2d `2W+2` bits,
//! orient3d `3W+3`, the fused quadratic form `3W-2+2 log2 n`) and `Circle`
//! (incircle `4W+5`, insphere `5W+6`). One fused kernel,
//! [`quadratic_form_exact`], and the three narrowings from `3F` fractional
//! bits down to `F`: floor and ceil for an enclosure, nearest with ties
//! toward positive infinity for a scalar (the binary house rule).
//!
//! Nothing here is public. The predicates, `Interval` and `fused` are the
//! surfaces; this module is the arithmetic they share so that the width
//! budget lives in one place.

use super::interval::exact_product;
use super::{FixedMatrix, FixedVector};
use crate::fixed_point::core_types::errors::OverflowDetected;
use crate::fixed_point::universal::fasc::stack_evaluator::{BinaryStorage, ComputeStorage};

#[cfg(table_format = "q16_16")]
use crate::fixed_point::frac_config;
#[cfg(table_format = "q32_32")]
use crate::fixed_point::i256::mul_i128_to_i256;
#[cfg(any(table_format = "q16_16", table_format = "q32_32"))]
use crate::fixed_point::I256;
#[cfg(any(table_format = "q32_32", table_format = "q64_64"))]
use crate::fixed_point::I512;
#[cfg(any(table_format = "q64_64", table_format = "q128_128"))]
use crate::fixed_point::I1024;
#[cfg(any(table_format = "q128_128", table_format = "q256_256"))]
use crate::fixed_point::I2048;

// ============================================================================
// Exact accumulator arithmetic
// ============================================================================

/// Exact signed integer arithmetic on a fixed-width accumulator.
///
/// Products are formed on magnitudes and the sign reapplied, so the wide
/// types' `Mul` implementations are only ever entered with non-negative
/// operands (a schoolbook or a truncated unsigned product, both exact there).
/// Each product first asserts that the operands' bit lengths fit the type,
/// which is the width budget made loud. Sums are checked: a running sum that
/// leaves the accumulator is a typed `TierOverflow`, never a wrap.
pub(crate) trait Wide:
    Copy + PartialEq + std::ops::Add<Output = Self> + std::ops::Sub<Output = Self> + std::ops::Mul<Output = Self> + std::ops::Neg<Output = Self>
{
    const BITS: u32;
    fn zero() -> Self;
    fn is_negative(self) -> bool;
    /// Significant bits of a NON-NEGATIVE value.
    fn bit_length(self) -> u32;
    /// Checked signed addition on the accumulator (`None` on overflow).
    fn checked_add_wide(self, rhs: Self) -> Option<Self>;

    #[inline]
    fn mul_exact(self, rhs: Self) -> Self {
        let neg = self.is_negative() != rhs.is_negative();
        let a = if self.is_negative() { -self } else { self };
        let b = if rhs.is_negative() { -rhs } else { rhs };
        assert!(
            a.bit_length() + b.bit_length() <= Self::BITS - 1,
            "exact predicate: accumulator width budget violated"
        );
        let p = a * b;
        if neg { -p } else { p }
    }

    #[inline]
    fn add_exact(self, rhs: Self) -> Result<Self, OverflowDetected> {
        self.checked_add_wide(rhs).ok_or(OverflowDetected::TierOverflow)
    }
}

#[inline]
fn words_bit_length(words: &[u64]) -> u32 {
    for i in (0..words.len()).rev() {
        if words[i] != 0 {
            return i as u32 * 64 + (64 - words[i].leading_zeros());
        }
    }
    0
}

impl Wide for i128 {
    const BITS: u32 = 128;
    #[inline] fn zero() -> Self { 0 }
    #[inline] fn is_negative(self) -> bool { self < 0 }
    #[inline] fn bit_length(self) -> u32 { 128 - self.leading_zeros() }
    #[inline] fn checked_add_wide(self, rhs: Self) -> Option<Self> { i128::checked_add(self, rhs) }
}

#[cfg(any(table_format = "q16_16", table_format = "q32_32"))]
impl Wide for I256 {
    const BITS: u32 = 256;
    #[inline] fn zero() -> Self { I256::zero() }
    #[inline] fn is_negative(self) -> bool { I256::is_negative(self) }
    #[inline] fn bit_length(self) -> u32 { words_bit_length(&self.words) }
    #[inline] fn checked_add_wide(self, rhs: Self) -> Option<Self> { I256::checked_add(self, rhs) }
}

#[cfg(any(table_format = "q32_32", table_format = "q64_64"))]
impl Wide for I512 {
    const BITS: u32 = 512;
    #[inline] fn zero() -> Self { I512::zero() }
    #[inline] fn is_negative(self) -> bool { I512::is_negative(self) }
    #[inline] fn bit_length(self) -> u32 { words_bit_length(&self.words) }
    #[inline] fn checked_add_wide(self, rhs: Self) -> Option<Self> { I512::checked_add(self, rhs) }
}

#[cfg(any(table_format = "q64_64", table_format = "q128_128"))]
impl Wide for I1024 {
    const BITS: u32 = 1024;
    #[inline] fn zero() -> Self { I1024::zero() }
    #[inline] fn is_negative(self) -> bool { (self.words[15] as i64) < 0 }
    #[inline] fn bit_length(self) -> u32 { words_bit_length(&self.words) }
    #[inline] fn checked_add_wide(self, rhs: Self) -> Option<Self> { I1024::checked_add(self, rhs) }
}

#[cfg(any(table_format = "q128_128", table_format = "q256_256"))]
impl Wide for I2048 {
    const BITS: u32 = 2048;
    #[inline] fn zero() -> Self { I2048::zero() }
    #[inline] fn is_negative(self) -> bool { (self.words[31] as i64) < 0 }
    #[inline] fn bit_length(self) -> u32 { words_bit_length(&self.words) }
    #[inline] fn checked_add_wide(self, rhs: Self) -> Option<Self> { I2048::checked_add(self, rhs) }
}

// Per-profile accumulator selection, from the storage width W:
//   orient2d 2W+2, orient3d 3W+3, quadratic form 3W-2+2log2(n) -> Orient;
//   incircle 4W+5, insphere 5W+6 -> Circle.
// `orient` widens a storage value exactly.
#[cfg(table_format = "q16_16")]
pub(crate) mod acc {
    pub type Orient = i128; // W = 32: 66 / 99 bits
    pub type Circle = super::I256; // 133 / 166 bits
    #[inline] pub fn orient(v: super::BinaryStorage) -> Orient { v as i128 }
    #[inline] pub fn circle(v: super::BinaryStorage) -> Circle { super::I256::from_i128(v as i128) }
}
#[cfg(table_format = "q32_32")]
pub(crate) mod acc {
    pub type Orient = super::I256; // W = 64: 130 / 195 bits
    pub type Circle = super::I512; // 261 / 326 bits
    #[inline] pub fn orient(v: super::BinaryStorage) -> Orient { super::I256::from_i128(v as i128) }
    #[inline] pub fn circle(v: super::BinaryStorage) -> Circle { super::I512::from_i128(v as i128) }
}
#[cfg(table_format = "q64_64")]
pub(crate) mod acc {
    pub type Orient = super::I512; // W = 128: 258 / 387 bits
    pub type Circle = super::I1024; // 517 / 646 bits
    #[inline] pub fn orient(v: super::BinaryStorage) -> Orient { super::I512::from_i128(v) }
    #[inline] pub fn circle(v: super::BinaryStorage) -> Circle { super::I1024::from_i128(v) }
}
#[cfg(table_format = "q128_128")]
pub(crate) mod acc {
    pub type Orient = super::I1024; // W = 256: 514 / 771 bits
    pub type Circle = super::I2048; // 1029 / 1286 bits
    #[inline] pub fn orient(v: super::BinaryStorage) -> Orient { super::I1024::from_i256(v) }
    #[inline] pub fn circle(v: super::BinaryStorage) -> Circle { super::I2048::from_i256(v) }
}
#[cfg(table_format = "q256_256")]
pub(crate) mod acc {
    pub type Orient = super::I2048; // W = 512: 1026 / 1539 bits; circle predicates exceed 2048
    #[inline] pub fn orient(v: super::BinaryStorage) -> Orient { super::I2048::from_i512(v) }
}

// ============================================================================
// The fused quadratic form: one exact value at 3F fractional bits
// ============================================================================

/// A storage value sign-extended to the compute width (no shift; exact).
#[inline]
pub(crate) fn widen_storage(v: BinaryStorage) -> ComputeStorage {
    #[cfg(table_format = "q16_16")]
    {
        v as i64
    }
    #[cfg(table_format = "q32_32")]
    {
        v as i128
    }
    #[cfg(table_format = "q64_64")]
    {
        crate::fixed_point::I256::from_i128(v)
    }
    #[cfg(table_format = "q128_128")]
    {
        crate::fixed_point::I512::from_i256(v)
    }
    #[cfg(table_format = "q256_256")]
    {
        crate::fixed_point::I1024::from_i512(v)
    }
}

/// Sign of a compute-width value (the compute types do not all expose
/// `is_negative`; the top word's sign bit is the definition).
#[inline]
#[cfg(any(table_format = "q64_64", table_format = "q128_128", table_format = "q256_256"))]
fn compute_is_negative(v: ComputeStorage) -> bool {
    (v.words[v.words.len() - 1] as i64) < 0
}

/// Exact product of two compute-width values on the accumulator (compute
/// width doubled). The wide `mul_to_*` family is UNSIGNED by convention, so
/// it is applied to magnitudes and the sign reapplied; the magnitudes'
/// bit lengths are asserted against the accumulator first (the width budget
/// made loud, as in [`Wide::mul_exact`]).
#[inline]
pub(crate) fn widen_product(a: ComputeStorage, b: ComputeStorage) -> acc::Orient {
    #[cfg(table_format = "q16_16")]
    {
        (a as i128) * (b as i128)
    }
    #[cfg(table_format = "q32_32")]
    {
        // sign-correct by itself
        mul_i128_to_i256(a, b)
    }
    #[cfg(any(table_format = "q64_64", table_format = "q128_128", table_format = "q256_256"))]
    {
        let (a_neg, b_neg) = (compute_is_negative(a), compute_is_negative(b));
        let neg = a_neg != b_neg;
        let abs_a = if a_neg { -a } else { a };
        let abs_b = if b_neg { -b } else { b };
        assert!(
            words_bit_length(&abs_a.words) + words_bit_length(&abs_b.words) <= acc::Orient::BITS - 1,
            "exact predicate: accumulator width budget violated"
        );
        #[cfg(table_format = "q64_64")]
        let p = abs_a.mul_to_i512(abs_b);
        #[cfg(table_format = "q128_128")]
        let p = abs_a.mul_to_i1024(abs_b);
        #[cfg(table_format = "q256_256")]
        let p = abs_a.mul_to_i2048(abs_b);
        if neg { -p } else { p }
    }
}

/// Exact `v^T M v` at `3F` fractional bits on the orientation accumulator.
///
/// Every term is an exact triple product: `v_i v_j` is formed exactly at the
/// compute tier (`2F`) and multiplied on the accumulator by the exact entry
/// `m_ii` (or, off the diagonal, by `m_ij + m_ji` at the compute width, which
/// pairs the two symmetric terms into one wide multiply; no symmetry of `M`
/// is assumed). The sum is checked. The value returned has not been rounded
/// at all; the caller narrows it once with [`narrow_triple_floor`],
/// [`narrow_triple_ceil`] or [`narrow_triple_nearest`].
///
/// Budget: `|v_i v_j| < 2^(2W-2)`, `|m_ij + m_ji| < 2^W`, so each product is
/// below `3W-2` bits and the sum of `n(n+1)/2` of them below
/// `3W-2+2 log2 n`; the accumulator has `4W` (or more) bits and every
/// multiply asserts its own budget.
///
/// Panics if `m` is not square or its size differs from `v`.
pub(crate) fn quadratic_form_exact(v: &FixedVector, m: &FixedMatrix) -> Result<acc::Orient, OverflowDetected> {
    let n = v.len();
    assert!(m.is_square() && m.rows() == n, "quadratic_form: dimension mismatch");
    let mut sum = acc::Orient::zero();
    for i in 0..n {
        let vi = v[i].raw();
        let vv = exact_product(vi, vi);
        sum = sum.add_exact(widen_product(widen_storage(m.get(i, i).raw()), vv))?;
        for j in (i + 1)..n {
            // W+1 bits at most: exact at the compute width (2W).
            let pair = widen_storage(m.get(i, j).raw()) + widen_storage(m.get(j, i).raw());
            let vv = exact_product(vi, v[j].raw());
            sum = sum.add_exact(widen_product(pair, vv))?;
        }
    }
    Ok(sum)
}

/// Split a `3F`-scaled accumulator into its floor at the storage scale, the
/// "any discarded bit set" flag and the round bit (bit `2F-1`).
///
/// The arithmetic right shift is floor in two's complement. The floor is
/// fits-checked against the storage tier; a value beyond it is a
/// `TierOverflow`.
#[inline]
fn split_triple(a: acc::Orient) -> Result<(BinaryStorage, bool, bool), OverflowDetected> {
    #[cfg(table_format = "q16_16")]
    {
        let shift = 2 * frac_config::FRAC_BITS;
        let mask = (1i128 << shift) - 1;
        let inexact = (a & mask) != 0;
        let round_bit = ((a >> (shift - 1)) & 1) == 1;
        let shifted = a >> shift;
        if shifted > i32::MAX as i128 || shifted < i32::MIN as i128 {
            return Err(OverflowDetected::TierOverflow);
        }
        Ok((shifted as i32, inexact, round_bit))
    }
    #[cfg(table_format = "q32_32")]
    {
        // I256 >> 64: the discarded bits are exactly words[0]
        let inexact = a.words[0] != 0;
        let round_bit = (a.words[0] >> 63) == 1;
        let shifted = a >> 64u32;
        if !shifted.fits_in_i128() {
            return Err(OverflowDetected::TierOverflow);
        }
        let s = shifted.as_i128();
        if s > i64::MAX as i128 || s < i64::MIN as i128 {
            return Err(OverflowDetected::TierOverflow);
        }
        Ok((s as i64, inexact, round_bit))
    }
    #[cfg(table_format = "q64_64")]
    {
        // I512 >> 128: the discarded bits are words[0..2]
        let inexact = (a.words[0] | a.words[1]) != 0;
        let round_bit = (a.words[1] >> 63) == 1;
        let shifted = a >> 128usize;
        if !shifted.fits_in_i128() {
            return Err(OverflowDetected::TierOverflow);
        }
        Ok((shifted.as_i128(), inexact, round_bit))
    }
    #[cfg(table_format = "q128_128")]
    {
        // I1024 >> 256: the discarded bits are words[0..4]
        let inexact = (a.words[0] | a.words[1] | a.words[2] | a.words[3]) != 0;
        let round_bit = (a.words[3] >> 63) == 1;
        let shifted = a >> 256usize;
        if !shifted.fits_in_i256() {
            return Err(OverflowDetected::TierOverflow);
        }
        Ok((shifted.as_i256(), inexact, round_bit))
    }
    #[cfg(table_format = "q256_256")]
    {
        // I2048 >> 512: the discarded bits are words[0..8]
        let mut low = 0u64;
        for w in &a.words[0..8] {
            low |= *w;
        }
        let inexact = low != 0;
        let round_bit = (a.words[7] >> 63) == 1;
        let shifted = a >> 512usize;
        if !shifted.fits_in_i512() {
            return Err(OverflowDetected::TierOverflow);
        }
        Ok((shifted.as_i512(), inexact, round_bit))
    }
}

#[inline]
fn storage_one() -> BinaryStorage {
    #[cfg(any(table_format = "q16_16", table_format = "q32_32", table_format = "q64_64"))]
    {
        1
    }
    #[cfg(table_format = "q128_128")]
    {
        crate::fixed_point::I256::from_i128(1)
    }
    #[cfg(table_format = "q256_256")]
    {
        crate::fixed_point::I512::from_i128(1)
    }
}

/// Narrow a `3F`-scaled exact value to the storage scale, rounding toward
/// negative infinity. Fits-checked; never wraps.
#[inline]
pub(crate) fn narrow_triple_floor(a: acc::Orient) -> Result<BinaryStorage, OverflowDetected> {
    let (floor, _, _) = split_triple(a)?;
    Ok(floor)
}

/// Narrow a `3F`-scaled exact value to the storage scale, rounding toward
/// positive infinity. The bump above the floor is checked: a floor at the
/// storage maximum has no ceiling in storage, and that is a `TierOverflow`.
#[inline]
pub(crate) fn narrow_triple_ceil(a: acc::Orient) -> Result<BinaryStorage, OverflowDetected> {
    let (floor, inexact, _) = split_triple(a)?;
    if inexact {
        floor.checked_add(storage_one()).ok_or(OverflowDetected::TierOverflow)
    } else {
        Ok(floor)
    }
}

/// Narrow a `3F`-scaled exact value to the storage scale, rounding to nearest
/// with ties toward positive infinity: floor plus the round bit, the binary
/// house rule, on every profile. The bump is checked.
#[inline]
pub(crate) fn narrow_triple_nearest(a: acc::Orient) -> Result<BinaryStorage, OverflowDetected> {
    let (floor, _, round_bit) = split_triple(a)?;
    if round_bit {
        floor.checked_add(storage_one()).ok_or(OverflowDetected::TierOverflow)
    } else {
        Ok(floor)
    }
}

/// `a / 2^shift` rounded to nearest with ties toward positive infinity (floor
/// plus the round bit), narrowed to storage (checked). For an exact product
/// whose scale carries a power-of-two normalisation, so the shift is not one
/// of the fixed `2F` or `3F`. (Realtime does the same on a plain `i128` in
/// `fused`.)
#[cfg(not(table_format = "q16_16"))]
pub(crate) fn narrow_shifted_nearest(a: acc::Orient, shift: u32) -> Result<BinaryStorage, OverflowDetected> {
    if shift == 0 {
        return narrow_orient_to_storage(a).ok_or(OverflowDetected::TierOverflow);
    }
    let round_bit = (a.words[((shift - 1) / 64) as usize] >> ((shift - 1) % 64)) & 1 == 1;
    #[cfg(table_format = "q32_32")]
    let floor = a >> shift;
    #[cfg(any(table_format = "q64_64", table_format = "q128_128", table_format = "q256_256"))]
    let floor = a >> shift as usize;
    let floor = narrow_orient_to_storage(floor).ok_or(OverflowDetected::TierOverflow)?;
    if round_bit {
        floor.checked_add(storage_one()).ok_or(OverflowDetected::TierOverflow)
    } else {
        Ok(floor)
    }
}

/// `num / den` rounded to nearest, ties toward positive infinity, narrowed
/// to storage (checked): `num` an exact wide value, `den` a nonzero compute
/// raw. With `num` at `3F` fractional bits and `den` at `2F` the quotient is
/// at the storage scale. Lets a caller divide an exact product by a compute
/// value without forming a quotient that may exceed the compute tier.
pub(crate) fn divide_to_storage_nearest(num: acc::Orient, den: ComputeStorage) -> Result<BinaryStorage, OverflowDetected> {
    narrow_orient_to_storage(divide_nearest(num, den)?).ok_or(OverflowDetected::TierOverflow)
}

/// `num / den` rounded to nearest, ties toward positive infinity, narrowed to
/// the compute tier (checked): with `num` at `4F` fractional bits and `den`
/// at `2F` the quotient is a compute raw. The compute-tier twin of
/// [`divide_to_storage_nearest`], for algorithms that keep their state at the
/// compute tier between steps.
pub(crate) fn divide_to_compute_nearest(num: acc::Orient, den: ComputeStorage) -> Result<ComputeStorage, OverflowDetected> {
    narrow_orient_to_compute(divide_nearest(num, den)?).ok_or(OverflowDetected::TierOverflow)
}

/// Exact `sum a_i b_i` of compute raws at `4F` fractional bits (checked).
pub(crate) fn exact_dot_compute(a: &[ComputeStorage], b: &[ComputeStorage]) -> Result<acc::Orient, OverflowDetected> {
    assert_eq!(a.len(), b.len(), "exact_dot_compute: length mismatch");
    let mut sum = <acc::Orient as Wide>::zero();
    for (x, y) in a.iter().zip(b) {
        sum = sum.add_exact(widen_product(*x, *y))?;
    }
    Ok(sum)
}

/// `init - sum a_i b_i` of compute raws: the sum exact at `4F`, `init`
/// widened exactly, one rounding to the compute tier (nearest, ties toward
/// +infinity, checked). The substitution and elimination kernel of the
/// compute-tier factorizations.
pub(crate) fn exact_sub_dot_compute(
    init: ComputeStorage, a: &[ComputeStorage], b: &[ComputeStorage],
) -> Result<ComputeStorage, OverflowDetected> {
    let one = crate::fixed_point::universal::fasc::stack_evaluator::compute::make_compute_int(1);
    let dot = exact_dot_compute(a, b)?;
    narrow_product_to_compute(widen_product(init, one).add_exact(-dot)?)
}

/// `num / den` in the accumulator type, nearest with ties toward +infinity.
fn divide_nearest(num: acc::Orient, den: ComputeStorage) -> Result<acc::Orient, OverflowDetected> {
    #[cfg(table_format = "q16_16")]
    let (q, r, d) = {
        let d = den as i128;
        if d == 0 { return Err(OverflowDetected::DivisionByZero); }
        (num / d, num % d, d)
    };
    #[cfg(table_format = "q32_32")]
    let (q, r, d) = {
        let d = I256::from_i128(den);
        if d == I256::zero() { return Err(OverflowDetected::DivisionByZero); }
        (num / d, num % d, d)
    };
    #[cfg(table_format = "q64_64")]
    let (q, r, d) = {
        let d = I512::from_i256(den);
        if d == I512::zero() { return Err(OverflowDetected::DivisionByZero); }
        (num / d, num % d, d)
    };
    #[cfg(table_format = "q128_128")]
    let (q, r, d) = {
        let d = I1024::from_i512(den);
        if d == I1024::zero() { return Err(OverflowDetected::DivisionByZero); }
        (num / d, num % d, d)
    };
    #[cfg(table_format = "q256_256")]
    let (q, r, d) = {
        let d = I2048::from_i1024(den);
        if d == I2048::zero() { return Err(OverflowDetected::DivisionByZero); }
        let (q, r) = crate::fixed_point::domains::binary_fixed::i2048::i2048_divmod(num, d);
        (q, r, d)
    };
    // truncated quotient; move one unit away from zero past the half (at the
    // half only for a positive quotient: ties toward +infinity)
    let positive = num.is_negative() == d.is_negative();
    let r_abs = if r.is_negative() { -r } else { r };
    let d_abs = if d.is_negative() { -d } else { d };
    let twice = r_abs + r_abs;
    let bump = if positive { twice >= d_abs } else { twice > d_abs };
    let q = if bump { if positive { q + unit_orient() } else { q - unit_orient() } } else { q };
    Ok(q)
}

#[inline]
fn unit_orient() -> acc::Orient {
    #[cfg(table_format = "q16_16")]
    { 1i128 }
    #[cfg(table_format = "q32_32")]
    { I256::from_i128(1) }
    #[cfg(table_format = "q64_64")]
    { I512::from_i128(1) }
    #[cfg(table_format = "q128_128")]
    { I1024::from_i128(1) }
    #[cfg(table_format = "q256_256")]
    { I2048::from_i128(1) }
}

/// An `Orient` value as a compute raw, `None` if it does not fit.
fn narrow_orient_to_compute(q: acc::Orient) -> Option<ComputeStorage> {
    #[cfg(table_format = "q16_16")]
    { i64::try_from(q).ok() }
    #[cfg(table_format = "q32_32")]
    { if q.fits_in_i128() { Some(q.as_i128()) } else { None } }
    #[cfg(table_format = "q64_64")]
    { if q.fits_in_i256() { Some(q.as_i256()) } else { None } }
    #[cfg(table_format = "q128_128")]
    { if q.fits_in_i512() { Some(q.as_i512()) } else { None } }
    #[cfg(table_format = "q256_256")]
    { if q.fits_in_i1024() { Some(q.as_i1024()) } else { None } }
}

/// An `Orient` value as storage, `None` if it does not fit.
fn narrow_orient_to_storage(q: acc::Orient) -> Option<BinaryStorage> {
    #[cfg(table_format = "q16_16")]
    { i32::try_from(q).ok() }
    #[cfg(table_format = "q32_32")]
    { if q.fits_in_i128() { i64::try_from(q.as_i128()).ok() } else { None } }
    #[cfg(table_format = "q64_64")]
    { if q.fits_in_i128() { Some(q.as_i128()) } else { None } }
    #[cfg(table_format = "q128_128")]
    { if q.fits_in_i256() { Some(q.as_i256()) } else { None } }
    #[cfg(table_format = "q256_256")]
    { if q.fits_in_i512() { Some(crate::fixed_point::I512::from_words([q.words[0], q.words[1], q.words[2], q.words[3], q.words[4], q.words[5], q.words[6], q.words[7]])) } else { None } }
}

/// Narrow an exact product of two compute raws (`4F` fractional bits, from
/// [`widen_product`]) to the compute scale (`2F`), rounding to nearest with
/// ties toward positive infinity. The result is fits-checked against the
/// compute width; a value beyond it is a `TierOverflow`, never a wrap.
#[inline]
pub(crate) fn narrow_product_to_compute(a: acc::Orient) -> Result<ComputeStorage, OverflowDetected> {
    #[cfg(table_format = "q16_16")]
    let (floor, round_bit, unit) = {
        let shift = 2 * frac_config::FRAC_BITS;
        let round_bit = ((a >> (shift - 1)) & 1) == 1;
        let shifted = a >> shift;
        if shifted > i64::MAX as i128 || shifted < i64::MIN as i128 {
            return Err(OverflowDetected::TierOverflow);
        }
        (shifted as i64, round_bit, 1i64)
    };
    #[cfg(table_format = "q32_32")]
    let (floor, round_bit, unit) = {
        // I256 >> 64: the discarded bits are exactly words[0]
        let round_bit = (a.words[0] >> 63) == 1;
        let shifted = a >> 64u32;
        if !shifted.fits_in_i128() {
            return Err(OverflowDetected::TierOverflow);
        }
        (shifted.as_i128(), round_bit, 1i128)
    };
    #[cfg(table_format = "q64_64")]
    let (floor, round_bit, unit) = {
        // I512 >> 128: the discarded bits are words[0..2]
        let round_bit = (a.words[1] >> 63) == 1;
        let shifted = a >> 128usize;
        if !shifted.fits_in_i256() {
            return Err(OverflowDetected::TierOverflow);
        }
        (shifted.as_i256(), round_bit, crate::fixed_point::I256::from_i128(1))
    };
    #[cfg(table_format = "q128_128")]
    let (floor, round_bit, unit) = {
        // I1024 >> 256: the discarded bits are words[0..4]
        let round_bit = (a.words[3] >> 63) == 1;
        let shifted = a >> 256usize;
        if !shifted.fits_in_i512() {
            return Err(OverflowDetected::TierOverflow);
        }
        (shifted.as_i512(), round_bit, crate::fixed_point::I512::from_i128(1))
    };
    #[cfg(table_format = "q256_256")]
    let (floor, round_bit, unit) = {
        // I2048 >> 512: the discarded bits are words[0..8]
        let round_bit = (a.words[7] >> 63) == 1;
        let shifted = a >> 512usize;
        if !shifted.fits_in_i1024() {
            return Err(OverflowDetected::TierOverflow);
        }
        (shifted.as_i1024(), round_bit, crate::fixed_point::I1024::from_i128(1))
    };
    if round_bit {
        floor.checked_add(unit).ok_or(OverflowDetected::TierOverflow)
    } else {
        Ok(floor)
    }
}


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

    fn one_third_raw() -> BinaryStorage {
        (FixedPoint::one() / FixedPoint::from_int(3)).raw()
    }

    /// A storage value `x` widened to the accumulator, times `2^(2F)`, narrows
    /// back to `x` under all three rules with nothing discarded.
    #[test]
    fn exact_multiples_narrow_identically_under_every_rule() {
        let x = one_third_raw();
        let v = FixedVector::from_slice(&[FixedPoint::one()]);
        let mut m = FixedMatrix::new(1, 1);
        m.set(0, 0, FixedPoint::from_raw(x));
        let a = quadratic_form_exact(&v, &m).unwrap();
        assert_eq!(narrow_triple_floor(a).unwrap(), x);
        assert_eq!(narrow_triple_ceil(a).unwrap(), x);
        assert_eq!(narrow_triple_nearest(a).unwrap(), x);
    }

    /// `v = [1/2]`, `m = [[2 ulp]]`: the exact value is exactly half an ulp,
    /// a tie. Floor 0, ceil 1, nearest rounds the tie toward +infinity: 1.
    /// The negative twin: floor -1, ceil 0, nearest 0.
    #[test]
    fn constructed_ties_round_toward_positive_infinity() {
        let half = FixedPoint::one() / FixedPoint::from_int(2);
        let two_ulp = FixedPoint::from_raw(storage_one() + storage_one());
        let v = FixedVector::from_slice(&[half]);
        let mut m = FixedMatrix::new(1, 1);
        m.set(0, 0, two_ulp);
        let a = quadratic_form_exact(&v, &m).unwrap();
        let zero = FixedPoint::ZERO.raw();
        assert_eq!(narrow_triple_floor(a).unwrap(), zero);
        assert_eq!(narrow_triple_ceil(a).unwrap(), storage_one());
        assert_eq!(narrow_triple_nearest(a).unwrap(), storage_one());
        m.set(0, 0, -two_ulp);
        let a = quadratic_form_exact(&v, &m).unwrap();
        assert_eq!(narrow_triple_floor(a).unwrap(), -storage_one());
        assert_eq!(narrow_triple_ceil(a).unwrap(), zero);
        assert_eq!(narrow_triple_nearest(a).unwrap(), zero);
    }

    /// Dividing by the compute-tier one (2^(2F)) is narrowing by 2^(2F): the
    /// quotient rounding must equal narrow_triple_nearest, ties included
    /// (both signs).
    #[test]
    fn divide_by_compute_one_matches_the_nearest_narrowing() {
        use crate::fixed_point::universal::fasc::stack_evaluator::compute::make_compute_int;
        let one_c = make_compute_int(1);
        let half = FixedPoint::one() / FixedPoint::from_int(2);
        let mut state: u64 = 0x5DEE_CE66_D1CE_4E5B;
        let mut cases = 0;
        // multipliers below 128: storage integers up to 24 fraction bits
        for k in -30i32..=30 {
            for tweak in [0i32, 1, 2, 3] {
                state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
                let ulps = (state >> 40) as i32 % 16 + k * 3 + tweak;
                let v = FixedVector::from_slice(&[half]);
                let mut m = FixedMatrix::new(1, 1);
                m.set(0, 0, FixedPoint::from_raw(storage_one()) * FixedPoint::from_int(ulps));
                let a = quadratic_form_exact(&v, &m).unwrap();
                assert_eq!(divide_to_storage_nearest(a, one_c), narrow_triple_nearest(a), "{ulps}");
                cases += 1;
            }
        }
        assert!(cases > 200);
    }

    /// The pairing `(m_ij + m_ji) v_i v_j` must not assume symmetry: an
    /// antisymmetric off-diagonal contributes nothing, a one-sided one
    /// contributes its full product.
    #[test]
    fn off_diagonal_terms_are_paired_without_assuming_symmetry() {
        let v = FixedVector::from_slice(&[FixedPoint::from_int(3), FixedPoint::from_int(5)]);
        let mut m = FixedMatrix::new(2, 2);
        m.set(0, 1, FixedPoint::from_int(7));
        m.set(1, 0, FixedPoint::from_int(-7));
        let a = quadratic_form_exact(&v, &m).unwrap();
        assert_eq!(narrow_triple_nearest(a).unwrap(), FixedPoint::ZERO.raw());
        m.set(1, 0, FixedPoint::ZERO);
        let a = quadratic_form_exact(&v, &m).unwrap();
        assert_eq!(narrow_triple_nearest(a).unwrap(), FixedPoint::from_int(105).raw());
    }

    /// A result beyond the storage tier is a typed overflow, never a wrap.
    #[test]
    fn a_result_beyond_storage_is_a_typed_overflow() {
        #[cfg(any(table_format = "q16_16", table_format = "q32_32", table_format = "q64_64"))]
        let big = FixedPoint::from_raw(BinaryStorage::MAX);
        #[cfg(table_format = "q128_128")]
        let big = FixedPoint::from_raw(crate::fixed_point::I256::max_value());
        #[cfg(table_format = "q256_256")]
        let big = FixedPoint::from_raw(crate::fixed_point::I512::max_value());
        let v = FixedVector::from_slice(&[big]);
        let mut m = FixedMatrix::new(1, 1);
        m.set(0, 0, big);
        let a = quadratic_form_exact(&v, &m).unwrap();
        assert_eq!(narrow_triple_floor(a), Err(OverflowDetected::TierOverflow));
        assert_eq!(narrow_triple_ceil(a), Err(OverflowDetected::TierOverflow));
        assert_eq!(narrow_triple_nearest(a), Err(OverflowDetected::TierOverflow));
    }
}