g_math 0.6.2

//! Linear algebra helpers with compute-tier precision.
//!
//! Core routines:
//! - `compute_tier_dot`: accumulates dot products at tier N+1 (double width)
//! - `compute_tier_sub_dot_raw`: fused init-minus-dot at compute tier
//! - `Rotation`, `householder_vector`, `reflect`: orthogonal transforms whose
//!   coefficients stay at the compute tier, each output narrowed once
//!
//! These are the matrix-operation analog of BinaryCompute chain persistence.

use super::FixedPoint;
use crate::fixed_point::universal::fasc::stack_evaluator::BinaryStorage;
#[cfg(table_format = "q16_16")]
use crate::fixed_point::frac_config;
#[cfg(table_format = "q64_64")]
use crate::fixed_point::I256;

#[cfg(table_format = "q128_128")]
use crate::fixed_point::{I256, I512};

#[cfg(table_format = "q256_256")]
use crate::fixed_point::{I512, I1024};

// Q32.32 and Q16.16 use native integer types only — no I256/I512/I1024 imports needed.
// BinaryStorage = i64 / i32, ComputeStorage = i128 / i64 (already native).

// Re-export ComputeStorage for fused operations
pub(crate) use crate::fixed_point::universal::fasc::stack_evaluator::ComputeStorage;
use crate::fixed_point::universal::fasc::stack_evaluator::compute::downscale_to_storage;

// Re-export fused sincos for imperative-module consumers (lie_group, etc.)
pub(crate) use crate::fixed_point::universal::fasc::stack_evaluator::compute::sincos_at_compute_tier;

// ============================================================================
// Rounding helper: round-to-nearest when downscaling from compute to storage
// ============================================================================

/// Downscale a compute-tier accumulator to storage tier with round-to-nearest.
/// This matches FASC's `downscale_to_storage` behavior, NOT truncation.
///
/// Panics if the result exceeds the storage tier's range, matching the
/// infallible imperative transcendentals. The previous shift-and-cast
/// fallback silently wrapped, which on narrow profiles turned an oversized
/// fused result into garbage (e.g. a negative squared distance).
#[inline]
pub(crate) fn round_to_storage(acc: ComputeStorage) -> BinaryStorage {
    downscale_to_storage(acc)
        .expect("round_to_storage: result exceeds storage tier range")
}

/// Upscale a storage-tier value to compute tier (shift left by FRAC_BITS).
#[inline]
pub(crate) fn upscale_to_compute(val: BinaryStorage) -> ComputeStorage {
    #[cfg(table_format = "q64_64")]
    { I256::from_i128(val) << 64usize }
    #[cfg(table_format = "q32_32")]
    { (val as i128) << 32 }
    #[cfg(table_format = "q16_16")]
    { (val as i64) << frac_config::FRAC_BITS }
    #[cfg(table_format = "q128_128")]
    { I512::from_i256(val) << 128usize }
    #[cfg(table_format = "q256_256")]
    { I1024::from_i512(val) << 256usize }
}

// ============================================================================
// Compute-tier dot product
// ============================================================================

/// Dot product accumulated at tier N+1 (compute tier).
///
/// Each product a_i * b_i is computed at double width without truncation.
/// The entire sum is accumulated at double width. Only one rounding step
/// occurs at the very end when truncating back to storage tier.
///
/// For an n-element dot product, this gives 1 ULP of rounding error instead
/// of the n ULP that storage-tier accumulation would produce.
///
/// Panics if the slices have different lengths.
pub fn compute_tier_dot(a: &[FixedPoint], b: &[FixedPoint]) -> FixedPoint {
    assert_eq!(a.len(), b.len(), "compute_tier_dot: length mismatch");

    #[cfg(table_format = "q64_64")]
    {
        // i128 × i128 → I256 (Q128.128), accumulate in I256, shift >> 64
        let mut acc = I256::zero();
        for i in 0..a.len() {
            let a_wide = I256::from_i128(a[i].raw());
            let b_wide = I256::from_i128(b[i].raw());
            // I256 * I256 → I256 (no overflow: inputs are sign-extended i128,
            // so the mathematical product fits in 256 bits)
            acc = acc + (a_wide * b_wide);
        }
        // Shift right by FRAC_BITS (64) to convert from Q128.128 to Q64.64
        FixedPoint::from_raw((acc >> 64u32).as_i128())
    }

    #[cfg(table_format = "q32_32")]
    {
        // i64 × i64 → i128 (Q64.64), accumulate in i128, shift >> 32
        let mut acc: i128 = 0;
        for i in 0..a.len() {
            acc += (a[i].raw() as i128) * (b[i].raw() as i128);
        }
        FixedPoint::from_raw((acc >> 32) as i64)
    }

    #[cfg(table_format = "q16_16")]
    {
        // i32 × i32 → i64, accumulate in i64, shift >> FRAC_BITS
        let mut acc: i64 = 0;
        for i in 0..a.len() {
            acc += (a[i].raw() as i64) * (b[i].raw() as i64);
        }
        FixedPoint::from_raw((acc >> frac_config::FRAC_BITS) as i32)
    }

    #[cfg(table_format = "q128_128")]
    {
        // I256 × I256 → I512 (Q256.256), accumulate in I512, shift >> 128
        let mut acc = I512::zero();
        for i in 0..a.len() {
            let a_raw = a[i].raw();
            let b_raw = b[i].raw();
            // Signed widening multiply: I256 × I256 → I512
            let a_neg = a_raw.is_negative();
            let b_neg = b_raw.is_negative();
            let result_neg = a_neg != b_neg;
            let abs_a = if a_neg { -a_raw } else { a_raw };
            let abs_b = if b_neg { -b_raw } else { b_raw };
            let product = abs_a.mul_to_i512(abs_b);
            let signed_product = if result_neg { -product } else { product };
            acc = acc + signed_product;
        }
        FixedPoint::from_raw((acc >> 128usize).as_i256())
    }

    #[cfg(table_format = "q256_256")]
    {
        // I512 × I512 → I1024 (Q512.512), accumulate in I1024, shift >> 256
        let mut acc = I1024::zero();
        for i in 0..a.len() {
            let a_raw = a[i].raw();
            let b_raw = b[i].raw();
            // Signed widening multiply: I512 × I512 → I1024
            let a_neg = a_raw.is_negative();
            let b_neg = b_raw.is_negative();
            let result_neg = a_neg != b_neg;
            let abs_a = if a_neg { -a_raw } else { a_raw };
            let abs_b = if b_neg { -b_raw } else { b_raw };
            let product = abs_a.mul_to_i1024(abs_b);
            let signed_product = if result_neg { -product } else { product };
            acc = acc + signed_product;
        }
        FixedPoint::from_raw((acc >> 256usize).as_i512())
    }
}

/// Compute-tier multiply-accumulate: acc += a_i * b_i for matrix operations.
///
/// Same as `compute_tier_dot` but takes raw BinaryStorage slices for
/// internal use where the FixedPoint wrapper would add unnecessary overhead.
#[inline]
pub(crate) fn compute_tier_dot_raw(a: &[BinaryStorage], b: &[BinaryStorage]) -> BinaryStorage {
    assert_eq!(a.len(), b.len(), "compute_tier_dot_raw: length mismatch");

    #[cfg(table_format = "q64_64")]
    {
        let mut acc = I256::zero();
        for i in 0..a.len() {
            acc = acc + (I256::from_i128(a[i]) * I256::from_i128(b[i]));
        }
        round_to_storage(acc)
    }

    #[cfg(table_format = "q32_32")]
    {
        // i64 × i64 → i128, accumulate in i128
        let mut acc: i128 = 0;
        for i in 0..a.len() {
            acc += (a[i] as i128) * (b[i] as i128);
        }
        round_to_storage(acc)
    }

    #[cfg(table_format = "q16_16")]
    {
        // i32 × i32 → i64, accumulate in i64
        let mut acc: i64 = 0;
        for i in 0..a.len() {
            acc += (a[i] as i64) * (b[i] as i64);
        }
        round_to_storage(acc)
    }

    #[cfg(table_format = "q128_128")]
    {
        let mut acc = I512::zero();
        for i in 0..a.len() {
            let a_neg = a[i].is_negative();
            let b_neg = b[i].is_negative();
            let result_neg = a_neg != b_neg;
            let abs_a = if a_neg { -a[i] } else { a[i] };
            let abs_b = if b_neg { -b[i] } else { b[i] };
            let product = abs_a.mul_to_i512(abs_b);
            acc = acc + if result_neg { -product } else { product };
        }
        round_to_storage(acc)
    }

    #[cfg(table_format = "q256_256")]
    {
        let mut acc = I1024::zero();
        for i in 0..a.len() {
            let a_neg = a[i].is_negative();
            let b_neg = b[i].is_negative();
            let result_neg = a_neg != b_neg;
            let abs_a = if a_neg { -a[i] } else { a[i] };
            let abs_b = if b_neg { -b[i] } else { b[i] };
            let product = abs_a.mul_to_i1024(abs_b);
            acc = acc + if result_neg { -product } else { product };
        }
        round_to_storage(acc)
    }
}

// ============================================================================
// Fused init-minus-dot at compute tier
// ============================================================================

/// Compute `init - dot(a, b)` entirely at compute tier (tier N+1).
///
/// Widens `init` to compute-tier format, subtracts each product a_i * b_i
/// accumulated at compute tier, then rounds back to storage tier once.
///
/// This is the core primitive for Gaussian elimination, forward/back
/// substitution, and Cholesky inner sums. It avoids the n ULP accumulation
/// error that would result from element-wise storage-tier operations.
///
/// Panics if `a.len() != b.len()`.
pub(crate) fn compute_tier_sub_dot_raw(
    init: BinaryStorage,
    a: &[BinaryStorage],
    b: &[BinaryStorage],
) -> BinaryStorage {
    assert_eq!(a.len(), b.len(), "compute_tier_sub_dot_raw: length mismatch");

    let acc = compute_tier_sub_dot_compute(init, a, b);
    round_to_storage(acc)
}

/// Same as `compute_tier_sub_dot_raw` but returns the result at compute tier
/// (ComputeStorage) WITHOUT downscaling. Used for fused operations where the
/// compute-tier intermediate feeds directly into sqrt or divide at compute tier.
pub(crate) fn compute_tier_sub_dot_compute(
    init: BinaryStorage,
    a: &[BinaryStorage],
    b: &[BinaryStorage],
) -> ComputeStorage {
    assert_eq!(a.len(), b.len(), "compute_tier_sub_dot_compute: length mismatch");

    #[cfg(table_format = "q64_64")]
    {
        let mut acc = I256::from_i128(init) << 64usize;
        for i in 0..a.len() {
            acc = acc - (I256::from_i128(a[i]) * I256::from_i128(b[i]));
        }
        acc
    }

    #[cfg(table_format = "q32_32")]
    {
        // i64 upscaled to i128, then subtract i64×i64→i128 products
        let mut acc: i128 = (init as i128) << 32;
        for i in 0..a.len() {
            acc -= (a[i] as i128) * (b[i] as i128);
        }
        acc
    }

    #[cfg(table_format = "q16_16")]
    {
        // i32 upscaled to i64, then subtract i32×i32→i64 products
        let mut acc: i64 = (init as i64) << frac_config::FRAC_BITS;
        for i in 0..a.len() {
            acc -= (a[i] as i64) * (b[i] as i64);
        }
        acc
    }

    #[cfg(table_format = "q128_128")]
    {
        let mut acc = I512::from_i256(init) << 128usize;
        for i in 0..a.len() {
            let a_neg = a[i].is_negative();
            let b_neg = b[i].is_negative();
            let result_neg = a_neg != b_neg;
            let abs_a = if a_neg { -a[i] } else { a[i] };
            let abs_b = if b_neg { -b[i] } else { b[i] };
            let product = abs_a.mul_to_i512(abs_b);
            acc = acc - if result_neg { -product } else { product };
        }
        acc
    }

    #[cfg(table_format = "q256_256")]
    {
        let mut acc = I1024::from_i512(init) << 256usize;
        for i in 0..a.len() {
            let a_neg = a[i].is_negative();
            let b_neg = b[i].is_negative();
            let result_neg = a_neg != b_neg;
            let abs_a = if a_neg { -a[i] } else { a[i] };
            let abs_b = if b_neg { -b[i] } else { b[i] };
            let product = abs_a.mul_to_i1024(abs_b);
            acc = acc - if result_neg { -product } else { product };
        }
        acc
    }
}

// ============================================================================
// Compute-tier orthogonal transforms for the iterative decompositions
// ============================================================================
//
// Jacobi, Golub-Kahan and Francis converge only if the transforms they apply
// inject less rounding noise than their convergence tests resolve. A rotation
// coefficient or Householder factor rounded to storage precision injects about
// |x| ulp into every entry it touches (the rotation is no longer orthogonal to
// storage precision). Here every coefficient stays at the compute tier and
// every transformed entry is narrowed once from an exact accumulator, so a
// transform adds at most about half an ulp per entry. Every step is checked:
// leaving the storage range is a `TierOverflow`, never a wrap.

use super::interval::exact_product;
use super::wide_acc::{acc, narrow_product_to_compute, narrow_triple_nearest, widen_product, widen_storage, Wide};
use crate::fixed_point::core_types::errors::OverflowDetected;
use crate::fixed_point::universal::fasc::stack_evaluator::compute::{
    compute_checked_add, compute_checked_divide, compute_divide, compute_is_negative,
    compute_is_zero, compute_multiply, compute_negate, make_compute_int, sqrt_at_compute_tier,
};

/// Consecutive non-improving iterations after which an iterative decomposition
/// is taken to sit at its precision floor.
pub(crate) const STAGNATION_SWEEPS: usize = 5;

/// Absolute noise floor of the iterative decompositions: four quanta.
///
/// An off-diagonal entry at or below it cannot be told apart from the rounding
/// the iterations inject, and dropping it moves every eigen- or singular value
/// by at most that much (Weyl). An exact zero eigen- or singular value is
/// computed as a block of such noise, which a purely relative test floored at
/// one quantum never deflates.
#[inline]
pub(crate) fn noise_floor() -> FixedPoint {
    FixedPoint::from_raw(noise_floor_raw())
}

/// Deflation bound for an entry beside diagonal magnitude `magnitude`: the
/// tight relative bound, floored at the absolute noise floor.
#[inline]
pub(crate) fn deflation_threshold(magnitude: FixedPoint) -> FixedPoint {
    convergence_threshold_tight(magnitude).max(noise_floor())
}

/// The looser sqrt(quantum) relative bound, floored at the noise floor,
/// accepted only once an iteration has stopped improving.
///
/// On realtime the shared [`convergence_threshold`] shifts by 8 whatever
/// `GMATH_FRAC_BITS` is. That is sqrt(quantum) at Q16.16 and looser above it,
/// but below it the bound is tighter than the precision floor: at Q22.10 it is
/// 4 ulps beside entries near one, under the rounding floor of a Francis step
/// there, and the iteration stalls until its budget runs out. So the shift is
/// `min(8, F/2)`.
#[inline]
pub(crate) fn stagnation_threshold(magnitude: FixedPoint) -> FixedPoint {
    #[cfg(table_format = "q16_16")]
    {
        let shift = (frac_config::FRAC_BITS / 2).min(8);
        FixedPoint::from_raw(magnitude.abs().raw() >> shift).max(noise_floor())
    }
    #[cfg(not(table_format = "q16_16"))]
    {
        convergence_threshold(magnitude).max(noise_floor())
    }
}

/// Magnitude of a compute-tier value.
#[inline]
pub(crate) fn compute_abs(v: ComputeStorage) -> ComputeStorage {
    if compute_is_negative(&v) { compute_negate(v) } else { v }
}

#[inline]
fn narrowed(value: acc::Orient) -> Result<FixedPoint, OverflowDetected> {
    Ok(FixedPoint::from_raw(narrow_triple_nearest(value)?))
}

/// `c x` for a compute-tier coefficient `c` and a storage value `x`, narrowed
/// once from the exact product.
#[inline]
pub(crate) fn scale_by(c: ComputeStorage, x: FixedPoint) -> Result<FixedPoint, OverflowDetected> {
    narrowed(widen_product(c, widen_storage(x.raw())))
}

/// A plane rotation `[cs sn; -sn cs]` whose coefficients stay at the compute
/// tier (`2 * FRAC_BITS` fractional bits).
#[derive(Clone, Copy)]
pub(crate) struct Rotation {
    cs: ComputeStorage,
    sn: ComputeStorage,
}

impl Rotation {
    /// A rotation from coefficients already at the compute tier.
    #[inline]
    pub(crate) fn from_parts(cs: ComputeStorage, sn: ComputeStorage) -> Self {
        Rotation { cs, sn }
    }

    /// The rotation taking `(a, b)` to `(r, 0)`: `cs = a / |r|`, `sn = b / |r|`
    /// (`cs = 1`, `sn = 0` when `b` is zero). `a` and `b` are compute raws at a
    /// common scale. Ratio form with one square root; neither input is squared,
    /// and no trig function is involved.
    pub(crate) fn zeroing_compute(a: ComputeStorage, b: ComputeStorage) -> Result<Self, OverflowDetected> {
        let one = make_compute_int(1);
        let zero = make_compute_int(0);
        if compute_is_zero(&b) {
            return Ok(Rotation { cs: one, sn: zero });
        }
        if compute_is_zero(&a) {
            let sn = if compute_is_negative(&b) { compute_negate(one) } else { one };
            return Ok(Rotation { cs: zero, sn });
        }
        if compute_abs(b) > compute_abs(a) {
            let tau = compute_divide(a, b)?;
            let inv = compute_divide(one, sqrt_at_compute_tier(compute_add(one, compute_multiply(tau, tau))))?;
            let sn = if compute_is_negative(&b) { compute_negate(inv) } else { inv };
            Ok(Rotation { cs: compute_multiply(sn, tau), sn })
        } else {
            let tau = compute_divide(b, a)?;
            let inv = compute_divide(one, sqrt_at_compute_tier(compute_add(one, compute_multiply(tau, tau))))?;
            let cs = if compute_is_negative(&a) { compute_negate(inv) } else { inv };
            Ok(Rotation { cs, sn: compute_multiply(cs, tau) })
        }
    }

    /// The rotation taking the storage values `(a, b)` to `(r, 0)`.
    #[inline]
    pub(crate) fn zeroing(a: FixedPoint, b: FixedPoint) -> Result<Self, OverflowDetected> {
        Self::zeroing_compute(upscale_to_compute(a.raw()), upscale_to_compute(b.raw()))
    }

    /// `(cs x + sn y, -sn x + cs y)`, each narrowed once from its exact value.
    pub(crate) fn apply(&self, x: FixedPoint, y: FixedPoint) -> Result<(FixedPoint, FixedPoint), OverflowDetected> {
        let (xw, yw) = (widen_storage(x.raw()), widen_storage(y.raw()));
        let first = widen_product(self.cs, xw).add_exact(widen_product(self.sn, yw))?;
        let second = widen_product(compute_negate(self.sn), xw).add_exact(widen_product(self.cs, yw))?;
        Ok((narrowed(first)?, narrowed(second)?))
    }

    /// `(cs x + sn y, -sn x + cs y)` on compute raws, each rounded once at the
    /// compute tier from its exact value.
    pub(crate) fn apply_compute(
        &self, x: ComputeStorage, y: ComputeStorage,
    ) -> Result<(ComputeStorage, ComputeStorage), OverflowDetected> {
        let first = widen_product(self.cs, x).add_exact(widen_product(self.sn, y))?;
        let second = widen_product(compute_negate(self.sn), x).add_exact(widen_product(self.cs, y))?;
        Ok((narrow_product_to_compute(first)?, narrow_product_to_compute(second)?))
    }

    /// `cs x + sn y` alone, on compute raws.
    pub(crate) fn combine_compute(&self, x: ComputeStorage, y: ComputeStorage) -> Result<ComputeStorage, OverflowDetected> {
        narrow_product_to_compute(widen_product(self.cs, x).add_exact(widen_product(self.sn, y))?)
    }

    /// `cs x` on a compute raw.
    #[inline]
    pub(crate) fn cos_times_compute(&self, x: ComputeStorage) -> Result<ComputeStorage, OverflowDetected> {
        compute_product(self.cs, x)
    }

    /// `sn x` on a compute raw.
    #[inline]
    pub(crate) fn sin_times_compute(&self, x: ComputeStorage) -> Result<ComputeStorage, OverflowDetected> {
        compute_product(self.sn, x)
    }

    /// `-sn x` on a compute raw.
    #[inline]
    pub(crate) fn neg_sin_times_compute(&self, x: ComputeStorage) -> Result<ComputeStorage, OverflowDetected> {
        compute_product(compute_negate(self.sn), x)
    }
}

/// `a b` of two compute raws, rounded once at the compute tier from the exact
/// product. A product beyond the compute tier is a `TierOverflow`.
#[inline]
pub(crate) fn compute_product(a: ComputeStorage, b: ComputeStorage) -> Result<ComputeStorage, OverflowDetected> {
    narrow_product_to_compute(widen_product(a, b))
}

/// Exact `sum a_i b_i` of storage values as a compute raw (`2 * FRAC_BITS`
/// fractional bits). A sum beyond the compute tier is a `TierOverflow`.
pub(crate) fn exact_dot(a: &[BinaryStorage], b: &[BinaryStorage]) -> Result<ComputeStorage, OverflowDetected> {
    assert_eq!(a.len(), b.len(), "exact_dot: length mismatch");
    let mut acc = make_compute_int(0);
    for i in 0..a.len() {
        acc = compute_checked_add(acc, exact_product(a[i], b[i]))?;
    }
    Ok(acc)
}

/// Householder direction `v = x - alpha e_1` for `x`, with
/// `alpha = -sign(x_0) ||x||` (no cancellation in `v_0`), and the exact `v.v`.
/// `None` when `x` is zero. `||x||` is a compute-tier square root of the exact
/// sum of squares, so a column whose squared norm exceeds the storage range is
/// still reflected; only a norm (or `v_0`) beyond storage is a `TierOverflow`.
pub(crate) fn householder_vector(
    x: &[BinaryStorage],
) -> Result<Option<(Vec<BinaryStorage>, ComputeStorage)>, OverflowDetected> {
    let xx = exact_dot(x, x)?;
    if compute_is_zero(&xx) {
        return Ok(None);
    }
    let norm = FixedPoint::from_raw(downscale_to_storage(sqrt_at_compute_tier(xx))?);
    let alpha = if FixedPoint::from_raw(x[0]).is_negative() { norm } else { -norm };
    let mut v = x.to_vec();
    v[0] = x[0].checked_sub(alpha.raw()).ok_or(OverflowDetected::TierOverflow)?;
    let vv = exact_dot(&v, &v)?;
    if compute_is_zero(&vv) {
        return Ok(None);
    }
    Ok(Some((v, vv)))
}

/// Reflect `w` in the hyperplane orthogonal to `v`: `w - 2 v (v.w) / (v.v)`.
///
/// `v` is an integer direction at any scale (the scale cancels) and `v_dot_v`
/// its exact `v.v`. The factor `2 (v.w) / (v.v)` is formed once at the compute
/// tier from exact sums, and each output entry is narrowed once from its exact
/// product with `v_k`.
pub(crate) fn reflect(
    w: &mut [BinaryStorage], v: &[BinaryStorage], v_dot_v: ComputeStorage,
) -> Result<(), OverflowDetected> {
    let vw = exact_dot(v, w)?;
    if compute_is_zero(&vw) {
        return Ok(());
    }
    let factor = compute_checked_divide(compute_checked_add(vw, vw)?, v_dot_v)?;
    for (wk, vk) in w.iter_mut().zip(v) {
        let update = narrow_triple_nearest(widen_product(factor, widen_storage(*vk)))?;
        *wk = wk.checked_sub(update).ok_or(OverflowDetected::TierOverflow)?;
    }
    Ok(())
}

#[cfg(table_format = "q64_64")]
fn noise_floor_raw() -> BinaryStorage { 4i128 }
#[cfg(table_format = "q32_32")]
fn noise_floor_raw() -> BinaryStorage { 4i64 }
#[cfg(table_format = "q16_16")]
fn noise_floor_raw() -> BinaryStorage { 4i32 }
#[cfg(table_format = "q128_128")]
fn noise_floor_raw() -> BinaryStorage { I256::from_i128(4) }
#[cfg(table_format = "q256_256")]
fn noise_floor_raw() -> BinaryStorage { I512::from_i128(4) }

// ============================================================================
// Convergence threshold for fixed-point iterative algorithms
// ============================================================================

/// Compute the convergence threshold for iterative algorithms.
///
/// In floating-point, convergence is tested against machine epsilon.
/// In fixed-point, we use `magnitude >> (FRAC_BITS / 2)`, which gives
/// sqrt(quantum) relative precision: the tightest achievable by iterative
/// multiply-based algorithms at storage tier.
///
/// Floored at 1 quantum (the smallest nonzero representable value).
///
/// Profile-dependent precision:
/// - Q64.64:  ~2^-32 relative (~9.5 decimal digits)
/// - Q128.128: ~2^-64 relative (~19 decimal digits)
/// - Q256.256: ~2^-128 relative (~38 decimal digits)
pub(crate) fn convergence_threshold(magnitude: FixedPoint) -> FixedPoint {
    let quantum = FixedPoint::from_raw(quantum_raw());
    let shifted = magnitude.abs().raw() >> half_frac_bits();
    let result = FixedPoint::from_raw(shifted);
    if result.is_zero() { quantum } else { result }
}

/// Tighter convergence threshold for compute-tier iterative algorithms.
///
/// When all rotation/accumulation steps happen at tier N+1 (via
/// `Rotation`, `reflect`, `compute_tier_dot_raw`), each step introduces
/// only 1 ULP of error, NOT √quantum. So we can converge to
/// `magnitude >> (2 * FRAC_BITS / 3)` instead of `>> (FRAC_BITS / 2)`.
///
/// Profile-dependent precision:
/// - Q64.64:  ~2^-42 relative (~12.6 decimal digits)
/// - Q128.128: ~2^-85 relative (~25.6 decimal digits)
/// - Q256.256: ~2^-170 relative (~51.2 decimal digits)
pub(crate) fn convergence_threshold_tight(magnitude: FixedPoint) -> FixedPoint {
    let quantum = FixedPoint::from_raw(quantum_raw());
    let shifted = magnitude.abs().raw() >> two_thirds_frac_bits();
    let result = FixedPoint::from_raw(shifted);
    if result.is_zero() { quantum } else { result }
}

#[cfg(table_format = "q64_64")]
fn two_thirds_frac_bits() -> u32 { 42 }
#[cfg(table_format = "q32_32")]
fn two_thirds_frac_bits() -> u32 { 21 }
#[cfg(table_format = "q16_16")]
fn two_thirds_frac_bits() -> u32 { 10 }
#[cfg(table_format = "q128_128")]
fn two_thirds_frac_bits() -> u32 { 85 }
#[cfg(table_format = "q256_256")]
fn two_thirds_frac_bits() -> usize { 170 }

#[cfg(table_format = "q64_64")]
fn half_frac_bits() -> u32 { 32 }
#[cfg(table_format = "q32_32")]
fn half_frac_bits() -> u32 { 16 }
#[cfg(table_format = "q16_16")]
fn half_frac_bits() -> u32 { 8 }
#[cfg(table_format = "q128_128")]
fn half_frac_bits() -> u32 { 64 }
#[cfg(table_format = "q256_256")]
fn half_frac_bits() -> usize { 128 }

#[cfg(table_format = "q64_64")]
fn quantum_raw() -> BinaryStorage { 1i128 }
#[cfg(table_format = "q32_32")]
fn quantum_raw() -> BinaryStorage { 1i64 }
#[cfg(table_format = "q16_16")]
fn quantum_raw() -> BinaryStorage { 1i32 }
#[cfg(table_format = "q128_128")]
fn quantum_raw() -> BinaryStorage { I256::from_i128(1) }
#[cfg(table_format = "q256_256")]
fn quantum_raw() -> BinaryStorage { I512::from_i128(1) }

// ============================================================================
// Compute-tier trit-weighted operations (zero-multiply dot product)
// Infrastructure ready for FASC ternary integration.
// ============================================================================

#[allow(unused_imports)]
use crate::fixed_point::domains::balanced_ternary::trit_packing::Trit;
// upscale_to_compute is defined locally above — no import needed

/// Trit-weighted dot product accumulated at compute tier (tier N+1).
///
/// For each trit in `packed_trits` (5 per byte, base-3 encoding):
///   - `+1` → `acc += widen(values[i])`
///   - `-1` → `acc -= widen(values[i])`
///   - ` 0` → skip (no operation)
///
/// The accumulator runs at ComputeStorage width. A single downscale
/// with rounding occurs at the very end. **Zero multiplications** in the
/// inner loop: only add/sub/skip.
///
/// After downscaling, the result is multiplied by `scale` (per-block
/// dequantization factor). The final multiply also uses compute-tier
/// intermediate to preserve precision.
///
/// # Arguments
/// - `packed_trits`: 5 trits per byte, base-3 encoded ({-1,0,+1} → {0,1,2})
/// - `num_elements`: exact number of trits (may be less than 5 × packed.len())
/// - `values`: activation vector in BinaryStorage format
/// - `scale`: per-block scale factor in BinaryStorage format
///
/// # Panics
/// Panics if `values.len() < num_elements`.
#[allow(dead_code)]
pub fn compute_tier_trit_dot_raw(
    packed_trits: &[u8],
    num_elements: usize,
    values: &[BinaryStorage],
    scale: BinaryStorage,
) -> BinaryStorage {
    assert!(values.len() >= num_elements, "compute_tier_trit_dot_raw: values shorter than num_elements");

    // Accumulate at compute tier (tier N+1) for full precision
    let mut acc = compute_zero();
    let mut trit_idx = 0;

    for &byte in packed_trits {
        if trit_idx >= num_elements {
            break;
        }

        // Unpack 5 trits from this byte (most-significant first)
        let mut remaining = byte;
        let mut chunk_trits = [1u8; 5]; // 1 = Zero (no-op)
        for j in (0..5).rev() {
            chunk_trits[j] = remaining % 3;
            remaining /= 3;
        }

        for j in 0..5 {
            if trit_idx >= num_elements {
                break;
            }

            let trit = chunk_trits[j];
            if trit == 2 {
                // Trit::Pos (+1): acc += widen(value)
                let widened = upscale_to_compute(values[trit_idx]);
                acc = compute_add(acc, widened);
            } else if trit == 0 {
                // Trit::Neg (-1): acc -= widen(value)
                let widened = upscale_to_compute(values[trit_idx]);
                acc = compute_sub(acc, widened);
            }
            // trit == 1 → Trit::Zero: skip (zero multiply eliminated)

            trit_idx += 1;
        }
    }

    // Single downscale of the accumulated dot product
    let dot_storage = round_to_storage(acc);

    // Apply per-block scale: result = dot * scale, at compute tier
    compute_tier_mul_pair(dot_storage, scale)
}

/// Row-wise trit-weighted matrix-vector product at compute tier.
///
/// Computes `result[row] = sum_j(trit[row][j] * values[j]) * scales[row]`
/// for each row, where the inner sum is a zero-multiply trit dot product.
///
/// # Arguments
/// - `packed_trits`: row-major packed trit matrix (each row = ceil(cols/5) bytes)
/// - `rows`: number of matrix rows
/// - `cols`: number of columns (= length of values vector)
/// - `values`: input activation vector
/// - `scales`: per-row scale factors (one per row)
///
/// # Returns
/// Output vector of length `rows`.
#[allow(dead_code)]
pub fn compute_tier_trit_matvec_raw(
    packed_trits: &[u8],
    rows: usize,
    cols: usize,
    values: &[BinaryStorage],
    scales: &[BinaryStorage],
) -> Vec<BinaryStorage> {
    assert!(values.len() >= cols, "compute_tier_trit_matvec_raw: values shorter than cols");
    assert!(scales.len() >= rows, "compute_tier_trit_matvec_raw: scales shorter than rows");

    let bytes_per_row = (cols + 4) / 5;
    let mut result = Vec::with_capacity(rows);

    for row in 0..rows {
        let row_start = row * bytes_per_row;
        let row_end = row_start + bytes_per_row;
        let row_trits = &packed_trits[row_start..row_end];

        let dot = compute_tier_trit_dot_raw(row_trits, cols, values, scales[row]);
        result.push(dot);
    }

    result
}

// Compute-tier helpers for trit operations
#[allow(dead_code)]
#[inline]
fn compute_zero() -> ComputeStorage {
    #[cfg(table_format = "q64_64")]
    { I256::zero() }
    #[cfg(table_format = "q32_32")]
    { 0i128 }
    #[cfg(table_format = "q16_16")]
    { 0i64 }
    #[cfg(table_format = "q128_128")]
    { I512::zero() }
    #[cfg(table_format = "q256_256")]
    { I1024::zero() }
}

#[allow(dead_code)]
#[inline]
fn compute_add(a: ComputeStorage, b: ComputeStorage) -> ComputeStorage {
    a + b
}

#[allow(dead_code)]
#[inline]
fn compute_sub(a: ComputeStorage, b: ComputeStorage) -> ComputeStorage {
    a - b
}

/// Multiply two BinaryStorage values at compute tier with single downscale.
#[allow(dead_code)]
#[inline]
fn compute_tier_mul_pair(a: BinaryStorage, b: BinaryStorage) -> BinaryStorage {
    #[cfg(table_format = "q64_64")]
    {
        let a_wide = I256::from_i128(a);
        let b_wide = I256::from_i128(b);
        let product = a_wide * b_wide;
        round_to_storage(product)
    }
    #[cfg(table_format = "q32_32")]
    {
        // i64 × i64 → i128 (native widening, no I256 needed)
        let product = (a as i128) * (b as i128);
        round_to_storage(product)
    }
    #[cfg(table_format = "q16_16")]
    {
        // i32 × i32 → i64 (native widening, no I128 needed)
        let product = (a as i64) * (b as i64);
        round_to_storage(product)
    }
    #[cfg(table_format = "q128_128")]
    {
        let a_neg = a.is_negative();
        let b_neg = b.is_negative();
        let result_neg = a_neg != b_neg;
        let abs_a = if a_neg { -a } else { a };
        let abs_b = if b_neg { -b } else { b };
        let product = abs_a.mul_to_i512(abs_b);
        let product = if result_neg { -product } else { product };
        round_to_storage(product)
    }
    #[cfg(table_format = "q256_256")]
    {
        let a_neg = a.is_negative();
        let b_neg = b.is_negative();
        let result_neg = a_neg != b_neg;
        let abs_a = if a_neg { -a } else { a };
        let abs_b = if b_neg { -b } else { b };
        let product = abs_a.mul_to_i1024(abs_b);
        let product = if result_neg { -product } else { product };
        round_to_storage(product)
    }
}