rusty_mp3 1.0.4

//! Hybrid synthesis stage 2: the 32-band polyphase synthesis filterbank.
//!
//! Each of the 18 passes takes 32 subband samples and produces 32 PCM samples.
//! A 32→64 cosine matrixing feeds the 1024-sample FIFO `V[]`; a windowing step
//! gathers `U[]` from `V[]`, multiplies by the 512-tap `D[]` window, and sums to
//! 32 outputs. The matrixing is computed here; `D[]` is tabulated in [`tables`](crate::tables).

use std::f64::consts::PI;
use std::sync::OnceLock;

use crate::frame::{GRANULE_LINES, SUBBANDS, SUBBAND_LINES};

use super::synth_window::SYNTH_D;

/// Matrixing coefficients `N[i][k] = cos((16+i)·(2k+1)·π/64)`, 64×32. The dense
/// reference — still the source of truth that [`matrixing_fast`] is gated against.
fn matrix() -> &'static [[f32; SUBBANDS]; 64] {
    static T: OnceLock<[[f32; SUBBANDS]; 64]> = OnceLock::new();
    T.get_or_init(|| {
        let mut n = [[0f32; SUBBANDS]; 64];
        for (i, row) in n.iter_mut().enumerate() {
            for (k, c) in row.iter_mut().enumerate() {
                *c = (PI / 64.0 * (16 + i) as f64 * (2 * k + 1) as f64).cos() as f32;
            }
        }
        n
    })
}

/// Half-width DCT kernel `C[m][k] = cos(m·(2k+1)·π/64)`, `m=0..31`, `k=0..15` —
/// the 32 distinct matrixing values `G[m]` after both symmetries are folded out.
fn half_dct() -> &'static [[f32; 16]; 32] {
    static T: OnceLock<[[f32; 16]; 32]> = OnceLock::new();
    T.get_or_init(|| {
        let mut c = [[0f32; 16]; 32];
        for (m, row) in c.iter_mut().enumerate() {
            for (k, e) in row.iter_mut().enumerate() {
                *e = (PI / 64.0 * m as f64 * (2 * k + 1) as f64).cos() as f32;
            }
        }
        c
    })
}

/// **B3** — the 64 matrixing outputs `V[i] = Σ_k cos((16+i)(2k+1)π/64)·s[k]`,
/// computed via two exact cosine symmetries instead of a dense 64×32 product:
///
/// * within `k`: `cos(m·(63−2k)π/64) = (−1)^m cos(m·(2k+1)π/64)` folds the 32
///   inputs into 16 sum/difference terms, so each `G[m]` is 16 mults not 32;
/// * across `i`: every `V[i]` is a signed copy of one of just 32 values
///   `G[m] = Σ_k cos(m(2k+1)π/64)·s[k]` (`m=0..31`) — `V[16]` is identically 0.
///
/// 512 mults vs 2048. Equal to [`matrix`] up to float reassociation (the symmetry
/// is exact; only ULP rounding differs) — pinned by `fast_matrixing_matches_dense`.
fn matrixing_fast(s: &[f32; SUBBANDS]) -> [f32; 64] {
    let c = half_dct();
    // Fold k↔31−k: even-m terms use the sum, odd-m terms the difference.
    let mut plus = [0f32; 16];
    let mut minus = [0f32; 16];
    for k in 0..16 {
        plus[k] = s[k] + s[31 - k];
        minus[k] = s[k] - s[31 - k];
    }
    // The 32 distinct DCT outputs G[0..31].
    let mut g = [0f32; 32];
    for (m, gm) in g.iter_mut().enumerate() {
        let src = if m & 1 == 0 { &plus } else { &minus };
        let mut acc = 0f32;
        for k in 0..16 {
            acc += c[m][k] * src[k];
        }
        *gm = acc;
    }
    // Map G → the 64 V outputs by sign/index (V[16] = 0).
    expand_g(&g)
}

/// `half_dct` transposed to `[k][m]`, so the 32 outputs `G[m]` are CONTIGUOUS
/// for a fixed `k`. Built once.
///
/// This is what makes the matrixing vectorisable without reassociating anything:
/// the natural loop reduces over `k` (a horizontal sum, which would need
/// reassociation and stop being bit-identical), but with `m` as the lane, each
/// `G[m]` accumulates over `k` in the ORIGINAL order, independently.
fn half_dct_t() -> &'static [[f32; 32]; 16] {
    static T: OnceLock<[[f32; 32]; 16]> = OnceLock::new();
    T.get_or_init(|| {
        let c = half_dct();
        let mut t = [[0f32; 32]; 16];
        for (m, row) in c.iter().enumerate() {
            for (k, v) in row.iter().enumerate() {
                t[k][m] = *v;
            }
        }
        t
    })
}

/// **D5** — AVX twin of [`matrixing_fast`], 8 of the 32 `G[m]` per lane.
///
/// Bit-identical for the same reason as [`window_avx`]: `m` is the lane, so each
/// output accumulates its 16 terms in the original `k` order, with separate
/// mul+add rather than FMA. The even/odd `m` split (even takes the sum-folded
/// input, odd the difference-folded one) becomes a single blend, because a
/// lane group always starts on an even `m`.
///
/// # Safety
/// Caller must have verified AVX is available.
#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx")]
unsafe fn matrixing_avx(s: &[f32; SUBBANDS], dst: &mut [f32; 64]) {
    use std::arch::x86_64::*;
    let ct = half_dct_t();
    let mut plus = [0f32; 16];
    let mut minus = [0f32; 16];
    for k in 0..16 {
        plus[k] = s[k] + s[31 - k];
        minus[k] = s[k] - s[31 - k];
    }
    // SAFETY: a register constant -- no memory access; the ISA is this fn's contract. Redundant from Rust 1.86 (safe target-feature calls), required at the crate's MSRV (1.85).
    #[allow(unused_unsafe)]
    let mut acc = [unsafe { _mm256_setzero_ps() }; 4];
    for k in 0..16 {
        // SAFETY: `ct[k]` is `[f32; 32]` (`half_dct_t`); `v * 8 + 8 <= 32` for `v < 4`.
        unsafe {
            // Lane m takes plus[k] when m is even, minus[k] when odd. Each group
            // of 8 starts at a multiple of 8, hence always even, so one mask does.
            let src = _mm256_blend_ps(_mm256_set1_ps(plus[k]), _mm256_set1_ps(minus[k]), 0xAA);
            for (v, accv) in acc.iter_mut().enumerate() {
                let c = _mm256_loadu_ps(ct[k].as_ptr().add(v * 8));
                *accv = _mm256_add_ps(*accv, _mm256_mul_ps(c, src));
            }
        }
    }
    // SAFETY: AVX per this fn's contract.
    unsafe { expand_g_avx(&acc, dst) };
}

/// **D14** — [`expand_g`] in registers: the 32 `G[m]` (four accumulators,
/// `G[8v..8v + 8]` in `acc[v]`) become the 64 `V` values by reversal, a
/// one-lane shift and sign flips, stored as eight full vectors.
///
/// Pure data movement: every output is one `G[m]`, `-G[m]` (an XOR of the sign
/// bit -- exactly what f32 negation is) or the `+0.0` of `V[16]`, the same as
/// the scalar map. What it replaces cost far more than its 64 moves: the scalar
/// map stored `V` one float at a time, and the next readers (the by-value
/// return, the FIFO copy, the window's loads of the newest slot) read it back
/// eight floats at a time, which the core cannot forward from eight separate
/// stores -- each such load waited for the stores to retire. Measured as a
/// ceiling probe (map stubbed out) it was ~25% of the whole decode.
/// Pinned by `expand_g_avx_matches_scalar`.
///
/// # Safety
/// Caller must have verified AVX is available.
#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx")]
unsafe fn expand_g_avx(acc: &[std::arch::x86_64::__m256; 4], dst: &mut [f32; 64]) {
    use std::arch::x86_64::*;
    // SAFETY: register-only operations; AVX per this fn's contract. The stores
    // write `dst[8t..8t + 8]` for `t < 8`, inside the `[f32; 64]`.
    unsafe {
        let sign = _mm256_set1_ps(-0.0);
        let zero = _mm256_setzero_ps();
        // Reverse eight lanes: within each 128-bit half, then swap the halves.
        let rev = |x: __m256| {
            let r = _mm256_permute_ps(x, 0x1B);
            _mm256_permute2f128_ps(r, r, 0x01)
        };
        // `[p7, x0, x1, .. x6]`: x moved up one lane, p's top lane shifted in.
        let shift1 = |x: __m256, p: __m256| {
            let t = _mm256_permute2f128_ps(x, p, 0x03); // [p.hi, x.lo]
            let rx = _mm256_permute_ps(x, 0x93); // per half [x3, x0, x1, x2]
            let rt = _mm256_permute_ps(t, 0x93); // per half lane 0 = t[3]
            _mm256_blend_ps(rx, rt, 0x11)
        };
        // N[v] = -reverse(G[8v..8v + 8]).
        let n: [__m256; 4] = [
            _mm256_xor_ps(rev(acc[0]), sign),
            _mm256_xor_ps(rev(acc[1]), sign),
            _mm256_xor_ps(rev(acc[2]), sign),
            _mm256_xor_ps(rev(acc[3]), sign),
        ];
        let d = dst.as_mut_ptr();
        _mm256_storeu_ps(d, acc[2]); // V[0..8]   = G[16..24]
        _mm256_storeu_ps(d.add(8), acc[3]); // V[8..16]  = G[24..32]
        _mm256_storeu_ps(d.add(16), shift1(n[3], zero)); // [+0, -G31 .. -G25]
        _mm256_storeu_ps(d.add(24), shift1(n[2], n[3])); // [-G24 .. -G17]
        _mm256_storeu_ps(d.add(32), shift1(n[1], n[2])); // [-G16 .. -G9]
        _mm256_storeu_ps(d.add(40), shift1(n[0], n[1])); // [-G8 .. -G1]
        _mm256_storeu_ps(d.add(48), _mm256_xor_ps(acc[0], sign)); // -G[0..8]
        _mm256_storeu_ps(d.add(56), _mm256_xor_ps(acc[1], sign)); // -G[8..16]
    }
}

/// NEON twin of [`matrixing_avx`], 4 of the 32 `G[m]` per lane. Each 4-lane
/// group starts on an even `m`, so the even/odd input split is the fixed lane
/// pattern `[plus, minus, plus, minus]`. Separate `vmulq` + `vaddq` keep it
/// bit-identical to `matrixing_fast` (pinned by `matrixing_simd_matches_scalar`,
/// run under qemu-aarch64).
///
/// # Safety
/// NEON is baseline on AArch64; all accesses are into fixed-size arrays.
#[cfg(target_arch = "aarch64")]
unsafe fn matrixing_neon(s: &[f32; SUBBANDS]) -> [f32; 64] {
    use std::arch::aarch64::*;
    let ct = half_dct_t();
    // SAFETY: a register constant -- no memory access; the ISA is this fn's contract. Redundant from Rust 1.86 (safe target-feature calls), required at the crate's MSRV (1.85).
    #[allow(unused_unsafe)]
    let mut acc = [unsafe { vdupq_n_f32(0.0) }; 8];
    for k in 0..16 {
        let (p, m) = (s[k] + s[31 - k], s[k] - s[31 - k]);
        let pm = [p, m, p, m];
        // SAFETY: `pm` is `[f32; 4]`; `ct[k]` is `[f32; 32]` and `v * 4 + 4 <= 32` for `v < 8`.
        unsafe {
            let src = vld1q_f32(pm.as_ptr());
            for (v, accv) in acc.iter_mut().enumerate() {
                let c = vld1q_f32(ct[k].as_ptr().add(v * 4));
                *accv = vaddq_f32(*accv, vmulq_f32(c, src));
            }
        }
    }
    let mut g = [0f32; 32];
    for (v, accv) in acc.iter().enumerate() {
        // SAFETY: `g` is `[f32; 32]`; `v * 4 + 4 <= 32` for `v < 8`.
        unsafe { vst1q_f32(g.as_mut_ptr().add(v * 4), *accv) };
    }
    expand_g(&g)
}

/// This arch's SIMD matrixing (AVX / NEON / scalar).
///
/// # Safety
/// Caller must have checked [`super::isa::simd_available`].
#[inline]
unsafe fn matrixing_simd(s: &[f32; SUBBANDS], dst: &mut [f32; 64]) {
    #[cfg(target_arch = "x86_64")]
    // SAFETY: forwards this fn's contract (the caller checked `simd_available()`).
    unsafe {
        matrixing_avx(s, dst);
    }
    #[cfg(target_arch = "aarch64")]
    // SAFETY: forwards this fn's contract (the caller checked `simd_available()`).
    unsafe {
        *dst = matrixing_neon(s);
    }
    #[cfg(not(any(target_arch = "x86_64", target_arch = "aarch64")))]
    {
        *dst = matrixing_fast(s);
    }
}

/// Map the 32 distinct DCT outputs onto the 64 `V` values by sign and index
/// (`V[16]` is identically zero). Shared by the scalar and AVX matrixing.
#[inline]
fn expand_g(g: &[f32; 32]) -> [f32; 64] {
    let mut vv = [0f32; 64];
    for (i, vi) in vv.iter_mut().enumerate() {
        *vi = if i < 16 {
            g[16 + i]
        } else if i == 16 {
            0.0
        } else if i < 48 {
            -g[48 - i]
        } else {
            -g[i - 48]
        };
    }
    vv
}

/// **D4** — AVX twin of the 16-tap windowing sum, 8 outputs per lane.
///
/// BIT-IDENTICAL to the scalar loop, not merely close: each output `j` stays in
/// its own lane and accumulates its 16 taps in the original order (i ascending,
/// the `a` term before the `b` term), and separate `mul`+`add` are used rather
/// than FMA — an FMA would round once where the scalar rounds twice. Pinned by
/// `window_simd_matches_scalar`.
///
/// Worth doing because auto-vectorization demonstrably did NOT happen here:
/// `--emit asm` on the scalar version counts 0 packed ops against 1698 scalar
/// ones, even after the operands were made contiguous.
///
/// # Safety
/// Caller must have verified AVX is available. All memory bounds are established
/// inside (fixed-size arrays; `head` is masked to a multiple of 64).
#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx")]
unsafe fn window_avx(fifo: &[f32; 1024], d: &[f32; 512], head: usize, out: &mut [f32; 32]) {
    use std::arch::x86_64::*;
    // `polyphase` only ever passes a multiple of 64; masking makes the bounds
    // below hold for ANY `head`, so the kernel's soundness is local (one `and`).
    let head = head & !63;
    // SAFETY: a register constant -- no memory access; the ISA is this fn's contract. Redundant from Rust 1.86 (safe target-feature calls), required at the crate's MSRV (1.85).
    #[allow(unused_unsafe)]
    let mut acc = [unsafe { _mm256_setzero_ps() }; 4];
    for i in 0..8 {
        // `a` is masked with 960, `b` with 1023, and the asymmetry is the point.
        // `head` starts at 0 and advances by `(head + 1024 - 64) & 1023`, so it is
        // always 64-aligned and `head + i * 128` has its low six bits clear -- the
        // narrower mask is the identity for `a` AND proves `a <= 960`, hence
        // `a + 32 <= 1024`, removing the span's bounds check. It is NOT valid for
        // `b`: `+ 96` is not 64-aligned (96 mod 64 = 32), so masking bit 5 away
        // there destroys signal. Tried; it moved the decode hash and failed two
        // reconstruction tests.
        let a = (head + i * 128) & 0x3C0; // 960
        let b = (head + i * 128 + 96) & 1023;
        for (v, accv) in acc.iter_mut().enumerate() {
            // SAFETY: `head` is masked to a multiple of 64 above, so `a <= 960` and `b` is 32 mod 64 with `b <= 992`: both 32-float spans (`+ v * 8 + 8 <= + 32`) end inside `fifo` (1024). `d` offsets reach `7 * 64 + 32 + 24 + 8 = 512`, its length.
            unsafe {
                let fa = _mm256_loadu_ps(fifo.as_ptr().add(a + v * 8));
                let da = _mm256_loadu_ps(d.as_ptr().add(i * 64 + v * 8));
                *accv = _mm256_add_ps(*accv, _mm256_mul_ps(fa, da));
                let fb = _mm256_loadu_ps(fifo.as_ptr().add(b + v * 8));
                let db = _mm256_loadu_ps(d.as_ptr().add(i * 64 + 32 + v * 8));
                *accv = _mm256_add_ps(*accv, _mm256_mul_ps(fb, db));
            }
        }
    }
    for (v, accv) in acc.iter().enumerate() {
        // SAFETY: `out` is `[f32; 32]`; `v * 8 + 8 <= 32` for `v < 4`.
        unsafe { _mm256_storeu_ps(out.as_mut_ptr().add(v * 8), *accv) };
    }
}

/// NEON twin of [`window_avx`]: 4 outputs per lane, 8 accumulators, the same tap
/// order (i ascending, `a` before `b`) with separate `vmulq` + `vaddq`. Pinned by
/// `window_simd_matches_scalar` under qemu-aarch64.
///
/// # Safety
/// NEON is baseline on AArch64; all memory bounds are established inside
/// (fixed-size arrays; `head` is masked to a multiple of 64).
#[cfg(target_arch = "aarch64")]
unsafe fn window_neon(fifo: &[f32; 1024], d: &[f32; 512], head: usize, out: &mut [f32; 32]) {
    use std::arch::aarch64::*;
    // See `window_avx`: masking makes the spans' bounds independent of the caller.
    let head = head & !63;
    // SAFETY: a register constant -- no memory access; the ISA is this fn's contract. Redundant from Rust 1.86 (safe target-feature calls), required at the crate's MSRV (1.85).
    #[allow(unused_unsafe)]
    let mut acc = [unsafe { vdupq_n_f32(0.0) }; 8];
    for i in 0..8 {
        let a = (head + i * 128) & 0x3C0; // 960
        let b = (head + i * 128 + 96) & 1023;
        for (v, accv) in acc.iter_mut().enumerate() {
            // SAFETY: `head` is masked to a multiple of 64 above, so `a <= 960` and `b <= 992`: both 32-float spans (`+ v * 4 + 4 <= + 32`) end inside `fifo` (1024). `d` offsets reach `7 * 64 + 32 + 28 + 4 = 512`, its length.
            unsafe {
                let fa = vld1q_f32(fifo.as_ptr().add(a + v * 4));
                let da = vld1q_f32(d.as_ptr().add(i * 64 + v * 4));
                *accv = vaddq_f32(*accv, vmulq_f32(fa, da));
                let fb = vld1q_f32(fifo.as_ptr().add(b + v * 4));
                let db = vld1q_f32(d.as_ptr().add(i * 64 + 32 + v * 4));
                *accv = vaddq_f32(*accv, vmulq_f32(fb, db));
            }
        }
    }
    for (v, accv) in acc.iter().enumerate() {
        // SAFETY: `out` is `[f32; 32]`; `v * 4 + 4 <= 32` for `v < 8`.
        unsafe { vst1q_f32(out.as_mut_ptr().add(v * 4), *accv) };
    }
}

/// This arch's SIMD windowing (AVX / NEON / scalar).
///
/// # Safety
/// Caller must have checked [`super::isa::simd_available`], and `head` must be a
/// multiple of 64.
#[inline]
unsafe fn window_simd(fifo: &[f32; 1024], d: &[f32; 512], head: usize, out: &mut [f32; 32]) {
    #[cfg(target_arch = "x86_64")]
    // SAFETY: forwards this fn's contract (the caller checked `simd_available()`).
    unsafe {
        window_avx(fifo, d, head, out);
    }
    #[cfg(target_arch = "aarch64")]
    // SAFETY: forwards this fn's contract (the caller checked `simd_available()`).
    unsafe {
        window_neon(fifo, d, head, out)
    }
    #[cfg(not(any(target_arch = "x86_64", target_arch = "aarch64")))]
    window_scalar(fifo, d, head, out)
}

/// The scalar windowing sum — the oracle the SIMD twin is gated against, and the
/// fallback on machines without AVX.
#[inline]
/// One 32-sample output block of the polyphase window.
///
/// `out` is a fixed-size block, not a slice: as `&mut [f32]` the compiler cannot
/// know it holds 32 elements, so each of the two accumulations per tap carried a
/// bounds check -- 512 of them per call, in the stage that is ~30% of decode.
///
/// REFUTED, kept as a warning: masking the FIFO offsets with 960 instead of 1023
/// looked equivalent -- `head` is 64-aligned, so `head + i * 128` has its low six
/// bits clear -- and would have proven `a + 32 <= 1024`, removing the span checks.
/// It is wrong for the SECOND offset: `+ 96` is not 64-aligned (96 mod 64 = 32), so
/// the narrower mask clears a bit that carries signal. Caught immediately by the
/// decode hash and two reconstruction tests.
fn window_scalar(fifo: &[f32; 1024], d: &[f32; 512], head: usize, out: &mut [f32; 32]) {
    out.fill(0.0);
    for i in 0..8 {
        // `a` is masked with 960, `b` with 1023, and the asymmetry is the point.
        // `head` starts at 0 and advances by `(head + 1024 - 64) & 1023`, so it is
        // always 64-aligned and `head + i * 128` has its low six bits clear -- the
        // narrower mask is the identity for `a` AND proves `a <= 960`, hence
        // `a + 32 <= 1024`, removing the span's bounds check. It is NOT valid for
        // `b`: `+ 96` is not 64-aligned (96 mod 64 = 32), so masking bit 5 away
        // there destroys signal. Tried; it moved the decode hash and failed two
        // reconstruction tests.
        let a = (head + i * 128) & 0x3C0; // 960
        let b = (head + i * 128 + 96) & 1023;
        let (fa, fb) = (&fifo[a..a + 32], &fifo[b..b + 32]);
        let (da, db) = (&d[i * 64..i * 64 + 32], &d[i * 64 + 32..i * 64 + 64]);
        for j in 0..32 {
            out[j] += fa[j] * da[j];
            out[j] += fb[j] * db[j];
        }
    }
}

/// Raw SIMD availability, for the oracle tests (they must run the kernel whatever
/// `MP3_ISA` says). The shipping path asks [`super::isa::use_simd`], resolved
/// once per call rather than per pass.
#[cfg(test)]
fn have_simd() -> bool {
    super::isa::simd_available()
}

/// Run the synthesis filterbank for one channel's granule (subband-major `time`),
/// returning 576 PCM samples. `fifo` is the persistent `V[]` state.
pub fn polyphase(time: &[f32; GRANULE_LINES], fifo: &mut [f32; 1024]) -> [f32; GRANULE_LINES] {
    // Resolved once per call, not per pass (codec-measurement: the A/B switch
    // itself is measurement overhead if it sits in the hot loop).
    let simd = super::isa::use_simd();
    super::prof::kernel_tally(
        if simd {
            &super::prof::K_SYNTH_SIMD
        } else {
            &super::prof::K_SYNTH_SCALAR
        },
        SUBBAND_LINES as u64,
    );
    // One `#[target_feature]` boundary per granule rather than two per pass: the
    // AVX entry compiles the whole body with AVX, so `matrixing_avx` and
    // `window_avx` inline instead of being called 36 times (see `imdct::hybrid`).
    #[cfg(target_arch = "x86_64")]
    if simd {
        // SAFETY: on x86_64 `use_simd()` is true only when AVX was detected.
        return unsafe { polyphase_avx(time, fifo) };
    }
    polyphase_impl(time, fifo, simd)
}

/// [`polyphase_impl`] compiled with AVX enabled, so the kernel calls inside inline.
///
/// # Safety
/// AVX must be available.
#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx")]
unsafe fn polyphase_avx(
    time: &[f32; GRANULE_LINES],
    fifo: &mut [f32; 1024],
) -> [f32; GRANULE_LINES] {
    polyphase_impl(time, fifo, true)
}

/// The synthesis body; `simd` selects the kernel twins.
// Load-bearing: the body must inline INTO the `#[target_feature]` entry so the
// kernels inline with it (docs/plans/mp3-kernel-ledger.md, brick 4).
#[allow(clippy::inline_always)]
#[inline(always)]
fn polyphase_impl(
    time: &[f32; GRANULE_LINES],
    fifo: &mut [f32; 1024],
    simd: bool,
) -> [f32; GRANULE_LINES] {
    let d = &SYNTH_D;
    let mut pcm = [0f32; GRANULE_LINES];

    // **D2** — the V FIFO is addressed circularly rather than shifted. Logical
    // index `L` (0 = newest) lives at physical `(head + L) & 1023`, so advancing
    // the FIFO by 64 is a subtraction on `head` instead of a 960-float memmove.
    //
    // The old form moved 960 floats on every one of the 18 passes — 17,280 float
    // moves per granule per channel, ~17.5 GB across a 7-minute track. Here the
    // only movement is one 1024-float rotation at the end, 16.9x less.
    //
    // Every span stays contiguous, so nothing is given up in exchange: `head` is
    // always a multiple of 64, `i*128` a multiple of 128, and 1024 a multiple of
    // 32 — so each 32-long run sits inside one 32-aligned block and never wraps
    // mid-run.
    let mut head = 0usize;

    for v in 0..SUBBAND_LINES {
        // Gather this pass's 32 subband samples.
        let mut s = [0f32; SUBBANDS];
        for (k, sv) in s.iter_mut().enumerate() {
            *sv = time[k * SUBBAND_LINES + v];
        }
        // Advance V by 64 and matrix the new 64 values into the front.
        head = (head + 1024 - 64) & 1023;
        let slot: &mut [f32; 64] = (&mut fifo[head..head + 64])
            .try_into()
            .expect("head is a multiple of 64 below 1024");
        if simd {
            // SAFETY: `simd` came from `use_simd()` above.
            unsafe { matrixing_simd(&s, slot) };
        } else {
            *slot = matrixing_fast(&s);
        }
        // Build U from V, window with D, sum 16 taps → one PCM sample per j.
        //
        // `i` outer, `j` inner: each of the four operands is then a CONTIGUOUS
        // 32-float run, which auto-vectorizes. The old `j`-outer form walked the
        // FIFO with a stride of 128 floats in its inner loop, which cannot
        // vectorize, and it recomputed the two circular indices 32 times over.
        //
        // Each `out[j]` still accumulates its 16 taps in the original order
        // (i ascending, the `a` term before the `b` term), so this is
        // bit-identical, not merely close.
        let out: &mut [f32; 32] = pcm[v * 32..]
            .first_chunk_mut()
            .expect("granule is 18 blocks of 32");
        if simd {
            // SAFETY: `simd` was resolved from `use_simd()` before the loop;
            // `head` is a multiple of 64 so both spans leave 32 floats in bounds.
            unsafe { window_simd(fifo, d, head, out) };
        } else {
            window_scalar(fifo, d, head, out);
        }
    }
    // Restore the canonical `head == 0` layout the caller's state is defined in:
    // 18 passes advanced head by 18*64 = 1152 ≡ 128 (mod 1024), leaving it at 896.
    fifo.rotate_left(head);
    pcm
}

/// **D12** — the synthesis filterbank for BOTH channels of a stereo granule at
/// once, writing interleaved PCM (`out[2 * n + ch]`).
///
/// Per channel this is exactly [`polyphase`]: the same passes, the same
/// matrixing and windowing sums in the same order. What changes is that the
/// two channels' passes run side by side in one AVX kernel, so each step keeps
/// eight independent accumulators live instead of four (the 16-deep add chains
/// were the bound), the shared cosine and window rows are loaded once for both
/// channels, and the interleave is two register shuffles per eight samples
/// instead of a separate pass over the granule. Pinned by
/// `polyphase_stereo_matches_two_mono_passes`.
pub fn polyphase_stereo(
    time: [&[f32; GRANULE_LINES]; 2],
    fifo: [&mut [f32; 1024]; 2],
    out: &mut [f32; 2 * GRANULE_LINES],
) {
    let simd = super::isa::use_simd();
    #[cfg(target_arch = "x86_64")]
    if simd {
        super::prof::kernel_tally(&super::prof::K_SYNTH_SIMD, 2 * SUBBAND_LINES as u64);
        // SAFETY: on x86_64 `use_simd()` is true only when AVX was detected.
        unsafe { polyphase_stereo_avx(time, fifo, out) };
        return;
    }
    // Elsewhere: two mono passes and an interleave (NEON keeps its kernels).
    let [f0, f1] = fifo;
    let a = polyphase(time[0], f0);
    let b = polyphase(time[1], f1);
    for ((o, &l), &r) in out.chunks_exact_mut(2).zip(&a).zip(&b) {
        o[0] = l;
        o[1] = r;
    }
}

/// The AVX body of [`polyphase_stereo`].
///
/// # Safety
/// AVX must be available.
#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx")]
unsafe fn polyphase_stereo_avx(
    time: [&[f32; GRANULE_LINES]; 2],
    fifo: [&mut [f32; 1024]; 2],
    out: &mut [f32; 2 * GRANULE_LINES],
) {
    let [f0, f1] = fifo;
    let d = &SYNTH_D;
    // See `polyphase_impl` (D2): circular FIFO, both channels advance together.
    let mut head = 0usize;
    for v in 0..SUBBAND_LINES {
        let mut s0 = [0f32; SUBBANDS];
        let mut s1 = [0f32; SUBBANDS];
        for k in 0..SUBBANDS {
            s0[k] = time[0][k * SUBBAND_LINES + v];
            s1[k] = time[1][k * SUBBAND_LINES + v];
        }
        head = (head + 1024 - 64) & 1023;
        let d0: &mut [f32; 64] = (&mut f0[head..head + 64])
            .try_into()
            .expect("64-aligned slot");
        let d1: &mut [f32; 64] = (&mut f1[head..head + 64])
            .try_into()
            .expect("64-aligned slot");
        // SAFETY: AVX per this fn's contract.
        unsafe { matrixing2_avx(&s0, &s1, d0, d1) };
        let o: &mut [f32; 64] = out[v * 64..]
            .first_chunk_mut()
            .expect("granule is 18 blocks of 2 x 32");
        // SAFETY: AVX per this fn's contract; `head` is a multiple of 64.
        unsafe { window2_avx(f0, f1, d, head, o) };
    }
    f0.rotate_left(head);
    f1.rotate_left(head);
}

/// [`matrixing_avx`] for two channels in one pass: per channel the identical
/// lane-wise sums (same `k` order, separate mul and add), with the cosine rows
/// loaded once for both.
///
/// # Safety
/// Caller must have verified AVX is available.
#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx")]
unsafe fn matrixing2_avx(
    s0: &[f32; SUBBANDS],
    s1: &[f32; SUBBANDS],
    d0: &mut [f32; 64],
    d1: &mut [f32; 64],
) {
    use std::arch::x86_64::*;
    let ct = half_dct_t();
    let (mut p0, mut m0, mut p1, mut m1) = ([0f32; 16], [0f32; 16], [0f32; 16], [0f32; 16]);
    for k in 0..16 {
        p0[k] = s0[k] + s0[31 - k];
        m0[k] = s0[k] - s0[31 - k];
        p1[k] = s1[k] + s1[31 - k];
        m1[k] = s1[k] - s1[31 - k];
    }
    // SAFETY: a register constant -- no memory access; the ISA is this fn's contract. Redundant from Rust 1.86 (safe target-feature calls), required at the crate's MSRV (1.85).
    #[allow(unused_unsafe)]
    let (mut a0, mut a1) = unsafe { ([_mm256_setzero_ps(); 4], [_mm256_setzero_ps(); 4]) };
    for k in 0..16 {
        // SAFETY: `ct[k]` is `[f32; 32]` (`half_dct_t`); `v * 8 + 8 <= 32` for `v < 4`.
        unsafe {
            // Even lanes take the sum-folded input, odd lanes the difference.
            let x0 = _mm256_blend_ps(_mm256_set1_ps(p0[k]), _mm256_set1_ps(m0[k]), 0xAA);
            let x1 = _mm256_blend_ps(_mm256_set1_ps(p1[k]), _mm256_set1_ps(m1[k]), 0xAA);
            for v in 0..4 {
                let c = _mm256_loadu_ps(ct[k].as_ptr().add(v * 8));
                a0[v] = _mm256_add_ps(a0[v], _mm256_mul_ps(c, x0));
                a1[v] = _mm256_add_ps(a1[v], _mm256_mul_ps(c, x1));
            }
        }
    }
    // SAFETY: AVX per this fn's contract.
    unsafe {
        expand_g_avx(&a0, d0);
        expand_g_avx(&a1, d1);
    }
}

/// [`window_avx`] for two channels, writing the 2 x 32 outputs interleaved:
/// per channel the identical 16-tap sums (`i` ascending, `a` before `b`,
/// separate mul and add), with the window rows loaded once for both.
///
/// # Safety
/// Caller must have verified AVX is available. Bounds are local, as in
/// [`window_avx`] (`head` is masked to a multiple of 64).
#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx")]
unsafe fn window2_avx(
    f0: &[f32; 1024],
    f1: &[f32; 1024],
    d: &[f32; 512],
    head: usize,
    out: &mut [f32; 64],
) {
    use std::arch::x86_64::*;
    let head = head & !63;
    // SAFETY: a register constant -- no memory access; the ISA is this fn's contract. Redundant from Rust 1.86 (safe target-feature calls), required at the crate's MSRV (1.85).
    #[allow(unused_unsafe)]
    let (mut a0, mut a1) = unsafe { ([_mm256_setzero_ps(); 4], [_mm256_setzero_ps(); 4]) };
    for i in 0..8 {
        // The offsets and masks of `window_avx` (`a <= 960`, `b <= 992`).
        let a = (head + i * 128) & 0x3C0;
        let b = (head + i * 128 + 96) & 1023;
        for v in 0..4 {
            // SAFETY: as in `window_avx`: both 32-float spans end inside each FIFO (1024), and the `d` offsets reach at most 512, its length.
            unsafe {
                let da = _mm256_loadu_ps(d.as_ptr().add(i * 64 + v * 8));
                let db = _mm256_loadu_ps(d.as_ptr().add(i * 64 + 32 + v * 8));
                let fa0 = _mm256_loadu_ps(f0.as_ptr().add(a + v * 8));
                let fa1 = _mm256_loadu_ps(f1.as_ptr().add(a + v * 8));
                a0[v] = _mm256_add_ps(a0[v], _mm256_mul_ps(fa0, da));
                a1[v] = _mm256_add_ps(a1[v], _mm256_mul_ps(fa1, da));
                let fb0 = _mm256_loadu_ps(f0.as_ptr().add(b + v * 8));
                let fb1 = _mm256_loadu_ps(f1.as_ptr().add(b + v * 8));
                a0[v] = _mm256_add_ps(a0[v], _mm256_mul_ps(fb0, db));
                a1[v] = _mm256_add_ps(a1[v], _mm256_mul_ps(fb1, db));
            }
        }
    }
    for v in 0..4 {
        // Interleave L/R: [l0 r0 l1 r1 .. l7 r7] for this lane group of eight.
        // SAFETY: `out` is `[f32; 64]`; `v * 16 + 16 <= 64` for `v < 4`.
        unsafe {
            let lo = _mm256_unpacklo_ps(a0[v], a1[v]);
            let hi = _mm256_unpackhi_ps(a0[v], a1[v]);
            _mm256_storeu_ps(
                out.as_mut_ptr().add(v * 16),
                _mm256_permute2f128_ps(lo, hi, 0x20),
            );
            _mm256_storeu_ps(
                out.as_mut_ptr().add(v * 16 + 8),
                _mm256_permute2f128_ps(lo, hi, 0x31),
            );
        }
    }
}

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

    /// **D14 gate.** The register expand is BIT-identical to the scalar
    /// `expand_g` -- every lane, signed zeros, infinities and NaN payloads
    /// included (negation must be a pure sign flip).
    #[test]
    #[cfg(target_arch = "x86_64")]
    fn expand_g_avx_matches_scalar() {
        if !have_simd() {
            eprintln!("no SIMD twin on this host - scalar path only, gate skipped");
            return;
        }
        let mut st = 0x510E_527Fu32;
        for trial in 0..512 {
            let mut g = [0f32; 32];
            for (m, x) in g.iter_mut().enumerate() {
                st ^= st << 13;
                st ^= st >> 17;
                st ^= st << 5;
                *x = match (st >> 28) & 7 {
                    0 => 0.0,
                    1 => -0.0,
                    2 => f32::from_bits(0x7FC0_0000 | (st & 0x3F_FFFF)),
                    3 => f32::NEG_INFINITY,
                    _ => (st >> 8) as f32 / (1u32 << 24) as f32 - 0.5 + m as f32,
                };
            }
            let want = expand_g(&g);
            let mut got = [0f32; 64];
            // SAFETY: gated on `have_simd()` above; loads are from `g` (32 floats).
            unsafe {
                use std::arch::x86_64::*;
                let acc = [
                    _mm256_loadu_ps(g.as_ptr()),
                    _mm256_loadu_ps(g.as_ptr().add(8)),
                    _mm256_loadu_ps(g.as_ptr().add(16)),
                    _mm256_loadu_ps(g.as_ptr().add(24)),
                ];
                expand_g_avx(&acc, &mut got);
            }
            for i in 0..64 {
                assert_eq!(got[i].to_bits(), want[i].to_bits(), "V[{i}] trial {trial}");
            }
        }
    }

    /// **D12 gate.** The two-channel synthesis is BIT-identical to two mono
    /// passes on the scalar oracle, interleaved -- over consecutive granules (the
    /// FIFOs carry), with signed zeros and silent subbands in the input.
    #[test]
    fn polyphase_stereo_matches_two_mono_passes() {
        let mut st = 0x6A09_E667u32;
        let mut next = move || {
            st ^= st << 13;
            st ^= st >> 17;
            st ^= st << 5;
            st
        };
        let (mut fa0, mut fa1) = ([0f32; 1024], [0f32; 1024]);
        let (mut fb0, mut fb1) = ([0f32; 1024], [0f32; 1024]);
        for g in 0..40 {
            let mut t = [[0f32; GRANULE_LINES]; 2];
            for (c, tc) in t.iter_mut().enumerate() {
                for (i, x) in tc.iter_mut().enumerate() {
                    let r = next();
                    *x = if (i / 18 + c + g) % 7 == 0 || r % 11 == 0 {
                        if r & 1 == 0 {
                            0.0
                        } else {
                            -0.0
                        }
                    } else {
                        (r >> 8) as f32 / (1u32 << 24) as f32 - 0.5
                    };
                }
            }
            let mut got = [0f32; 2 * GRANULE_LINES];
            polyphase_stereo([&t[0], &t[1]], [&mut fa0, &mut fa1], &mut got);
            let l = polyphase_impl(&t[0], &mut fb0, false);
            let r = polyphase_impl(&t[1], &mut fb1, false);
            for n in 0..GRANULE_LINES {
                assert_eq!(got[2 * n].to_bits(), l[n].to_bits(), "L[{n}] granule {g}");
                assert_eq!(
                    got[2 * n + 1].to_bits(),
                    r[n].to_bits(),
                    "R[{n}] granule {g}"
                );
            }
            assert!(fa0
                .iter()
                .zip(&fb0)
                .all(|(x, y)| x.to_bits() == y.to_bits()));
            assert!(fa1
                .iter()
                .zip(&fb1)
                .all(|(x, y)| x.to_bits() == y.to_bits()));
        }
    }

    /// **D5 gate.** The AVX matrixing must be BIT-identical to the scalar
    /// `matrixing_fast`. It can be, because `m` is the lane: each of the 32
    /// outputs accumulates its 16 terms in the original `k` order, so nothing is
    /// reassociated. (Reducing over `k` instead — the natural SIMD shape — would
    /// be a horizontal sum and would NOT be bit-identical.)
    #[test]
    fn matrixing_simd_matches_scalar() {
        if !have_simd() {
            eprintln!("no SIMD twin on this host - scalar path only, gate skipped");
            return;
        }
        let mut st = 0x3C6E_F372u32;
        let mut rng = || {
            st = st.wrapping_mul(1_664_525).wrapping_add(1_013_904_223);
            (st >> 8) as f32 / (1u32 << 24) as f32 - 0.5
        };
        for trial in 0..64 {
            let mut s = [0f32; SUBBANDS];
            for v in &mut s {
                *v = rng();
            }
            let a = matrixing_fast(&s);
            // SAFETY: gated on `have_simd()` above.
            let mut b = [0f32; 64];
            unsafe { matrixing_simd(&s, &mut b) };
            assert_eq!(a, b, "matrixing SIMD/scalar mismatch on trial {trial}");
        }
    }

    /// **D4 gate.** The AVX windowing twin must be BIT-identical to the scalar
    /// oracle — not within a tolerance. It can be, because each output stays in
    /// its own lane and accumulates in the original order, and the kernel uses
    /// separate mul+add rather than FMA (an FMA rounds once where the scalar
    /// rounds twice, which would change the result).
    ///
    /// Swept over every legal `head` (multiples of 64) so the circular wrap is
    /// covered, with random FIFO contents.
    #[test]
    fn window_simd_matches_scalar() {
        if !have_simd() {
            eprintln!("no SIMD twin on this host - scalar path only, gate skipped");
            return;
        }
        let mut st = 0x7F4A_7C15u32;
        let mut rng = || {
            st = st.wrapping_mul(1_664_525).wrapping_add(1_013_904_223);
            (st >> 8) as f32 / (1u32 << 24) as f32 - 0.5
        };
        let mut fifo = [0f32; 1024];
        for f in &mut fifo {
            *f = rng();
        }
        let d = &SYNTH_D;
        for head in (0..1024).step_by(64) {
            let mut a = [0f32; 32];
            let mut b = [0f32; 32];
            window_scalar(&fifo, d, head, &mut a);
            // SAFETY: gated on `have_simd()`; `head` steps by 64.
            unsafe { window_simd(&fifo, d, head, &mut b) };
            assert_eq!(a, b, "SIMD/scalar mismatch at head={head}");
        }
    }

    #[test]
    fn matrix_matches_cosine_formula() {
        let n = matrix();
        // Spot-check a couple of entries against N[i][k] = cos((16+i)(2k+1)π/64).
        let e =
            |i: usize, k: usize| (PI / 64.0 * (16 + i) as f64 * (2 * k + 1) as f64).cos() as f32;
        assert!((n[0][0] - e(0, 0)).abs() < 1e-6);
        assert!((n[33][7] - e(33, 7)).abs() < 1e-6);
        assert!((n[63][31] - e(63, 31)).abs() < 1e-6);
    }

    /// The dense 64×32 matrixing — the reference `matrixing_fast` must reproduce.
    fn matrixing_dense(s: &[f32; SUBBANDS]) -> [f32; 64] {
        let n = matrix();
        let mut v = [0f32; 64];
        for (i, vi) in v.iter_mut().enumerate() {
            let mut acc = 0f32;
            for k in 0..SUBBANDS {
                acc += n[i][k] * s[k];
            }
            *vi = acc;
        }
        v
    }

    /// **B3 gate.** The fast matrixing must equal the dense reference to float
    /// precision (the symmetries are exact; only ULP reassociation differs) — for
    /// many random subband inputs. A bug in the index/sign map shows up as a large
    /// error here, long before it could reach the FFmpeg-conformance check.
    #[test]
    #[cfg_attr(
        miri,
        ignore = "end-to-end / exhaustive: too slow to interpret under Miri (> 5 min)"
    )]
    fn fast_matrixing_matches_dense() {
        let mut st = 0x1234_9ABCu32;
        let mut rng = || {
            st = st.wrapping_mul(1_664_525).wrapping_add(1_013_904_223);
            (st >> 8) as f32 / (1u32 << 24) as f32 * 2.0 - 1.0
        };
        let mut worst = 0f32;
        for _ in 0..2000 {
            let mut s = [0f32; SUBBANDS];
            for sv in &mut s {
                *sv = rng() * 100.0;
            }
            let dense = matrixing_dense(&s);
            let fast = matrixing_fast(&s);
            for i in 0..64 {
                // relative error vs the row's scale (sum of |s|).
                let scale = s.iter().map(|x| x.abs()).sum::<f32>().max(1.0);
                worst = worst.max((dense[i] - fast[i]).abs() / scale);
            }
        }
        eprintln!("[B3] worst relative matrixing error fast vs dense: {worst:.2e}");
        assert!(
            worst < 1e-5,
            "fast matrixing diverges from dense: {worst:.2e}"
        );
    }

    #[test]
    fn fifo_advances_and_output_is_finite() {
        // Feed every pass (all-ones) so the FIFO stays populated across the 18
        // shifts. With the placeholder D (zeros) the PCM is zero but finite; the
        // matrixing/FIFO must run cleanly and leave the FIFO non-empty.
        let time = [1f32; GRANULE_LINES];
        let mut fifo = [0f32; 1024];
        let pcm = polyphase(&time, &mut fifo);
        assert!(pcm.iter().all(|v| v.is_finite()));
        assert!(
            fifo.iter().any(|&v| v != 0.0),
            "matrixing must fill the FIFO"
        );
    }
}