wacore 0.7.0

//! MLow LPC ANALYSIS FRONT-END (encoder side): `smpl_window` / `smpl_lpc` / `smpl_bwe_expand` plus
//! `smpl_A2NLSF_16` (wrapping silk `silk_A2NLSF`). It turns a 20 ms windowed LPC buffer into the
//! post-bandwidth-expansion LPC coefficients `A[0..16]` and their analysis NLSF, which feed the
//! (bit-exact) LSF quantizer.
//!
//! The autocorrelation is derived NOT by time-domain lags but by an FFT of the windowed signal
//! (512-pt real FFT) -> power spectrum -> a cosine transform (`brute_dct`) -> reflection coefficients
//! (`smpl_ac2rc`) -> `A` (`smpl_rc2a`), then bandwidth expansion (`A[i] *= 0.9999^i`). A time-domain
//! autocorr + Gaussian lag window would produce wrong formants and mute the outbound audio.
//!
//! Float tolerance: the reference uses PFFFT; our portable mixed-radix FFT (shared with `smpl_perc`)
//! differs only in FFT-internal rounding, so the resulting `A`/NLSF match to a tight float tolerance,
//! not bit-for-bit (documented like the decoder postfilters).
#![allow(clippy::needless_range_loop)]

use super::smpl_perc::{FftScratch, rfft_forward_ordered_sc};

pub(crate) const SMPL_LPC_ORDER: usize = 16;
const SMPL_LPC_NFFT: usize = 512;
// The reference uses this exact truncated literal (`SMPL_PI 3.1415926535897f`); keep it verbatim for
// bit-faithful window/DCT generation rather than `std::f32::consts::PI`.
#[allow(clippy::approx_constant, clippy::excessive_precision)]
const SMPL_PI: f32 = 3.1415926535897;
const SMPL_PI_F64: f64 = SMPL_PI as f64;
const SMPL_LPC_REG: f32 = 5e-7;
const SMPL_LPC_BWE: f32 = 0.9999;

const SMPL_LPC_WIN1_20MS_LEN: usize = 264;
const SMPL_WIN3_LONG_LEN: usize = 64;
const SMPL_WIN3_SHORT_LEN: usize = 32;
/// 20 ms LPC analysis buffer length (`lpcbuf_len` for 60 ms packets).
pub(crate) const SMPL_LPC_BUF_LEN: usize = 448;

// window generation (gen_sin_win / gen_cos_win)

/// `gen_sin_win`: `win[i] = sinf((i+1)/(N+1) * PI/2)`.
fn gen_sin_win(n: usize) -> Vec<f32> {
    (0..n)
        .map(|i| ((i as f32 + 1.0) / (n as f32 + 1.0) * SMPL_PI / 2.0).sin())
        .collect()
}

/// `gen_cos_win`: `win[i] = cosf((i+1)/(N+1) * PI/2)`.
fn gen_cos_win(n: usize) -> Vec<f32> {
    (0..n)
        .map(|i| ((i as f32 + 1.0) / (n as f32 + 1.0) * SMPL_PI / 2.0).cos())
        .collect()
}

/// `smpl_window` for the LPC 20 ms path (`use_lpc_win=true`, `frame_ms=20`, `len=448`).
/// `use_long_win` selects the 64-tap vs 32-tap trailing cosine window; the short variant zeros the
/// trailing 32 samples. Produces the windowed buffer the FFT consumes.
pub(crate) fn smpl_window_lpc20(
    input: &[f32; SMPL_LPC_BUF_LEN],
    use_long_win: bool,
) -> [f32; SMPL_LPC_BUF_LEN] {
    let win1 = gen_sin_win(SMPL_LPC_WIN1_20MS_LEN);
    let (win3, win3len) = if use_long_win {
        (gen_cos_win(SMPL_WIN3_LONG_LEN), SMPL_WIN3_LONG_LEN)
    } else {
        (gen_cos_win(SMPL_WIN3_SHORT_LEN), SMPL_WIN3_SHORT_LEN)
    };
    let mut out = [0.0f32; SMPL_LPC_BUF_LEN];
    for i in 0..SMPL_LPC_WIN1_20MS_LEN {
        out[i] = input[i] * win1[i];
    }
    // Untouched middle section copies straight through (no window taper).
    let mid = SMPL_LPC_BUF_LEN - SMPL_LPC_WIN1_20MS_LEN - SMPL_WIN3_LONG_LEN;
    out[SMPL_LPC_WIN1_20MS_LEN..SMPL_LPC_WIN1_20MS_LEN + mid]
        .copy_from_slice(&input[SMPL_LPC_WIN1_20MS_LEN..SMPL_LPC_WIN1_20MS_LEN + mid]);
    let base = SMPL_LPC_BUF_LEN - SMPL_WIN3_LONG_LEN;
    for i in 0..win3len {
        out[base + i] = input[base + i] * win3[i];
    }
    if !use_long_win {
        for s in out
            .iter_mut()
            .take(base + SMPL_WIN3_LONG_LEN)
            .skip(base + SMPL_WIN3_SHORT_LEN)
        {
            *s = 0.0;
        }
    }
    out
}

// brute_dct cosine tables (gen_cos_row)

const NFFT4: usize = SMPL_LPC_NFFT / 4; // 128

/// `gen_cos_row`: `row[k] = cos(omega)*scale` with `omega` accumulated via `fmod(omega+domega, 2pi)`
/// in f64 (NOT `cos(k*domega)`; the running fmod differs in float rounding).
fn gen_cos_row(domega: f64, scale: f64) -> [f64; NFFT4] {
    let mut row = [0.0f64; NFFT4];
    let mut omega = 0.0f64;
    let two_pi = 2.0 * SMPL_PI_F64;
    for k in 0..NFFT4 {
        row[k] = omega.cos() * scale;
        omega = (omega + domega).rem_euclid(two_pi);
    }
    row
}

struct DctTables {
    cdif: [[f64; NFFT4]; SMPL_LPC_ORDER / 2],     // 8 rows
    csumdiff: [[f64; NFFT4]; SMPL_LPC_ORDER / 4], // 4 rows
    csumsum: [[f64; NFFT4]; SMPL_LPC_ORDER / 4],  // 4 rows
}

fn build_dct_tables() -> DctTables {
    let two_pi = 2.0 * SMPL_PI_F64;
    let mut cdif = [[0.0f64; NFFT4]; SMPL_LPC_ORDER / 2];
    for j in 0..SMPL_LPC_ORDER / 2 {
        cdif[j] = gen_cos_row(
            (1 + j * 2) as f64 * two_pi / SMPL_LPC_NFFT as f64,
            2.0 / SMPL_LPC_NFFT as f64,
        );
    }
    let mut csumdiff = [[0.0f64; NFFT4]; SMPL_LPC_ORDER / 4];
    for j in 0..SMPL_LPC_ORDER / 4 {
        csumdiff[j] = gen_cos_row(
            (2 + j * 4) as f64 * two_pi / SMPL_LPC_NFFT as f64,
            1.0 / SMPL_LPC_NFFT as f64,
        );
    }
    let mut csumsum = [[0.0f64; NFFT4]; SMPL_LPC_ORDER / 4];
    for j in 0..SMPL_LPC_ORDER / 4 {
        csumsum[j] = gen_cos_row(
            (4 + j * 4) as f64 * two_pi / SMPL_LPC_NFFT as f64,
            1.0 / SMPL_LPC_NFFT as f64,
        );
    }
    DctTables {
        cdif,
        csumdiff,
        csumsum,
    }
}

/// `brute_dct`: autocorrelation `R[0..16]` from the power spectrum `F2[0..256]` via the precomputed
/// cosine sums. All accumulation in f64 (F2 is cast to double first).
fn brute_dct(t: &DctTables, f2: &[f64], order: usize, r: &mut [f64]) {
    let half = SMPL_LPC_NFFT / 2; // 256
    let mut f2sum = 0.0f64;
    let mut f2_dif = [0.0f64; NFFT4];
    let mut f2_sumsum = [0.0f64; NFFT4];
    let mut f2_sumdif = [0.0f64; NFFT4];
    for n in 0..NFFT4 {
        f2sum += f2[n] + f2[NFFT4 + n];
        f2_dif[n] = f2[n] - f2[half - n];
        f2_sumsum[n] = f2[n] + f2[half - n] + f2[NFFT4 + n] + f2[NFFT4 - n];
        f2_sumdif[n] = f2[n] + f2[half - n] - f2[NFFT4 + n] - f2[NFFT4 - n];
    }
    f2_dif[0] *= 0.5;
    r[0] = (2.0 * f2sum - f2[0] + f2[half]) / SMPL_LPC_NFFT as f64;

    for j in 0..order / 2 {
        let mut rtmp = 0.0f64;
        let row = &t.cdif[j];
        for k in 0..NFFT4 {
            rtmp += row[k] * f2_dif[k];
        }
        r[1 + j * 2] = rtmp;
    }
    for j in 0..order / 4 {
        let mut rtmp = 0.0f64;
        let row = &t.csumdiff[j];
        for k in 0..NFFT4 {
            rtmp += row[k] * f2_sumdif[k];
        }
        r[2 + j * 4] = rtmp;
    }
    for j in 0..order / 4 {
        let mut rtmp = 0.0f64;
        let row = &t.csumsum[j];
        for k in 0..NFFT4 {
            rtmp += row[k] * f2_sumsum[k];
        }
        r[4 + j * 4] = rtmp;
    }
}

/// `smpl_ac2rc_dbl`: autocorrelation `R[0..order]` -> reflection coefficients (Schur), with
/// `C0[0] *= (1 + reg)`. `rc[k]` is truncated to f32 each step (load-bearing for bit-faithfulness).
fn ac2rc_dbl(corr: &[f64], order: usize, reg: f32, rc: &mut [f32]) {
    let mut c0 = vec![0.0f64; order + 1];
    let mut c1 = vec![0.0f64; order + 1];
    c0[..order + 1].copy_from_slice(&corr[..order + 1]);
    c0[0] *= (1.0f32 + reg) as f64;
    c1[..order + 1].copy_from_slice(&c0[..order + 1]);
    for r in rc.iter_mut().take(order) {
        *r = 0.0;
    }
    for k in 0..order {
        if c0[k + 1] > c1[0] {
            rc[k] = -1.0;
            break;
        }
        if c0[k + 1] < -c1[0] {
            rc[k] = 1.0;
            break;
        }
        if c1[0] == 0.0 {
            break;
        }
        let rc_tmp = -c0[k + 1] / c1[0];
        rc[k] = rc_tmp as f32;
        for n in 0..(order - k) {
            let ctmp1 = c0[n + k + 1];
            let ctmp2 = c1[n];
            c0[n + k + 1] = ctmp1 + ctmp2 * rc_tmp;
            c1[n] = ctmp2 + ctmp1 * rc_tmp;
        }
    }
}

/// `smpl_rc2a`: reflection coefficients -> monic LPC `A[0..order]` (`A[0]=1`). All f32.
fn rc2a(rc: &[f32], order: usize, a: &mut [f32]) {
    for v in a.iter_mut().take(order + 1).skip(1) {
        *v = 0.0;
    }
    a[0] = 1.0;
    for k in 0..order {
        let rc_tmp = rc[k];
        for n in 0..k.div_ceil(2) {
            let tmp1 = a[n + 1];
            let tmp2 = a[k - n];
            a[n + 1] = tmp1 + tmp2 * rc_tmp;
            a[k - n] = tmp2 + tmp1 * rc_tmp;
        }
        a[k + 1] = rc_tmp;
    }
}

/// `smpl_bwe_expand` (the `bwe>0` path): `A[i] *= bwe^i` for i in 1..=order.
fn bwe_expand(a: &mut [f32], order: usize, bwe: f32) {
    let mut c = bwe;
    for i in 1..order + 1 {
        a[i] *= c;
        c *= bwe;
    }
}

/// Number of power-spectrum bins `smpl_pitch` / `smpl_get_signal_mode` read (`SMPL_F_LEN`).
pub(crate) const SMPL_F_LEN: usize = SMPL_LPC_NFFT / 2 + 1;

/// Reusable scratch for the fixed 512-pt LPC-analysis FFT. Held in the encoder state so the twiddle
/// tables are built once per stream instead of rebuilt on every internal frame.
pub(crate) fn new_lpc_fft_scratch() -> FftScratch {
    FftScratch::new(SMPL_LPC_NFFT)
}

/// Cached DCT cos tables: deterministic in the fixed NFFT/order, built once instead of recomputing
/// ~2k cos() per LPC analysis.
fn dct_tables() -> &'static DctTables {
    static TABLES: std::sync::OnceLock<DctTables> = std::sync::OnceLock::new();
    TABLES.get_or_init(build_dct_tables)
}

/// Full `smpl_lpc` + `smpl_bwe_expand` over a 20 ms windowed LPC buffer (`SMPL_LPC_BUF_LEN`),
/// returning the post-BWE monic LPC `A[0..16]` (`A[0]=1`) and the power spectrum `F2[0..256]`
/// (`lpcbuf_F2`) that `smpl_pitch` (harm_strength) and `smpl_get_signal_mode` (spectral tilt) consume.
pub(crate) fn smpl_lpc_analyze_with_f2(
    windowed: &[f32; SMPL_LPC_BUF_LEN],
    fft: &mut FftScratch,
) -> ([f32; SMPL_LPC_ORDER + 1], [f32; SMPL_F_LEN]) {
    // Zero-pad to NFFT and forward real FFT (pffft ordered layout).
    let mut xbuf = [0.0f32; SMPL_LPC_NFFT];
    xbuf[..SMPL_LPC_BUF_LEN].copy_from_slice(windowed);
    let mut f = [0.0f32; SMPL_LPC_NFFT];
    rfft_forward_ordered_sc(&xbuf, &mut f, fft);

    // Power spectrum F2[0..256] (float), pffft packs DC in F[0], Nyquist in F[1].
    let mut f2 = [0.0f32; SMPL_F_LEN];
    f2[0] = f[0] * f[0];
    f2[SMPL_LPC_NFFT / 2] = f[1] * f[1];
    for i in 1..SMPL_LPC_NFFT / 2 {
        f2[i] = f[2 * i] * f[2 * i] + f[2 * i + 1] * f[2 * i + 1];
    }
    let f2d: Vec<f64> = f2.iter().map(|&v| v as f64).collect();

    let mut r = [0.0f64; SMPL_LPC_ORDER + 1];
    brute_dct(dct_tables(), &f2d, SMPL_LPC_ORDER, &mut r);

    let mut rc = [0.0f32; SMPL_LPC_ORDER];
    ac2rc_dbl(&r, SMPL_LPC_ORDER, SMPL_LPC_REG, &mut rc);
    let mut a = [0.0f32; SMPL_LPC_ORDER + 1];
    rc2a(&rc, SMPL_LPC_ORDER, &mut a);
    bwe_expand(&mut a, SMPL_LPC_ORDER, SMPL_LPC_BWE);
    (a, f2)
}

// LPC interpolation (smpl_lpc_interpol)

/// `smpl_lsf_interpol_4[idx]`: per-subframe interpolation weights (idx 0 and the alternative idx 1
/// that `lsf_interpol_search` may pick when it lowers the residual energy).
const SMPL_LSF_INTERPOL_4_TBL: [[f32; SMPL_SUBFRS]; 2] =
    [[0.55, 0.88, 1.0, 1.0], [0.3, 0.65, 0.95, 1.0]];
const SMPL_SUBFRS: usize = 4;
const MAX_RC_STABLE: f32 = 0.9995;

/// `smpl_lpc_is_stable`: true if the monic LPC `A[0..16]` (`A[0]=1`) is a stable all-pole filter.
fn lpc_is_stable(a: &[f32]) -> bool {
    let order = SMPL_LPC_ORDER;
    if a[order] * a[order] > MAX_RC_STABLE {
        return false;
    }
    let mut a0 = [0.0f64; SMPL_LPC_ORDER];
    for i in 0..order {
        a0[i] = a[i + 1] as f64;
    }
    let mut a1 = [0.0f64; SMPL_LPC_ORDER];
    let mut m = order - 1;
    loop {
        let den = 1.0 - a0[m] * a0[m];
        if den == 0.0 {
            return false;
        }
        let inv = 1.0 / den;
        for k in 0..m {
            a1[k] = (a0[k] - a0[m] * a0[m - k - 1]) * inv;
        }
        if a1[m - 1] * a1[m - 1] > MAX_RC_STABLE as f64 {
            return false;
        }
        if m == 1 {
            return true;
        }
        m -= 1;
        let den = 1.0 - a1[m] * a1[m];
        if den == 0.0 {
            return false;
        }
        let inv = 1.0 / den;
        for k in 0..m {
            a0[k] = (a1[k] - a1[m] * a1[m - k - 1]) * inv;
        }
        if a0[m - 1] * a0[m - 1] > MAX_RC_STABLE as f64 {
            return false;
        }
        if m == 1 {
            return true;
        }
        m -= 1;
    }
}

/// `smpl_lpc_stabilize`: bandwidth-expand until stable.
fn lpc_stabilize(a: &mut [f32]) {
    if lpc_is_stable(a) {
        return;
    }
    let mut iter = 0;
    loop {
        iter += 1;
        bwe_expand(a, SMPL_LPC_ORDER, 1.0 - iter as f32 * 0.001);
        if lpc_is_stable(a) {
            return;
        }
    }
}

/// Per-subframe interpolated LPC predcoefs (`smpl_lpc_interpol`, lsf_interpol_idx 0), and the updated
/// `prev_lsf` (the last subframe's interpolated NLSF, carried to the next internal frame). `nlsf2a`
/// is the decoder's NLSF->A. On a reset (`prev_lsf` all zero) `lsf` is copied into `prev_lsf`.
pub(crate) fn smpl_lpc_interpol(
    lsf: &[f32],
    prev_lsf: &[f32],
    nlsf2a: impl Fn(&[f32]) -> Vec<f32>,
) -> (
    [[f32; SMPL_LPC_ORDER + 1]; SMPL_SUBFRS],
    [f32; SMPL_LPC_ORDER],
) {
    smpl_lpc_interpol_idx(lsf, prev_lsf, 0, nlsf2a)
}

/// As [`smpl_lpc_interpol`] but for an explicit `lsf_interpol_idx` (the `smpl_lsf_interpol_4` row).
pub(crate) fn smpl_lpc_interpol_idx(
    lsf: &[f32],
    prev_lsf: &[f32],
    interpol_idx: usize,
    nlsf2a: impl Fn(&[f32]) -> Vec<f32>,
) -> (
    [[f32; SMPL_LPC_ORDER + 1]; SMPL_SUBFRS],
    [f32; SMPL_LPC_ORDER],
) {
    let interp = &SMPL_LSF_INTERPOL_4_TBL[interpol_idx.min(1)];
    let mut prev = [0.0f32; SMPL_LPC_ORDER];
    if prev_lsf.len() == SMPL_LPC_ORDER && prev_lsf[SMPL_LPC_ORDER - 1] != 0.0 {
        prev.copy_from_slice(prev_lsf);
    } else {
        prev.copy_from_slice(&lsf[..SMPL_LPC_ORDER]);
    }
    let mut predcoefs = [[0.0f32; SMPL_LPC_ORDER + 1]; SMPL_SUBFRS];
    let mut ilsf = [0.0f32; SMPL_LPC_ORDER];
    for j in 0..SMPL_SUBFRS {
        let w = interp[j];
        if w == 1.0 {
            ilsf.copy_from_slice(&lsf[..SMPL_LPC_ORDER]);
        } else {
            for k in 0..SMPL_LPC_ORDER {
                ilsf[k] = (1.0 - w) * prev[k] + w * lsf[k];
            }
        }
        let a = nlsf2a(&ilsf);
        let mut pc = [0.0f32; SMPL_LPC_ORDER + 1];
        for (i, &c) in a.iter().enumerate().take(SMPL_LPC_ORDER + 1) {
            pc[i] = c;
        }
        pc[0] = 1.0;
        lpc_stabilize(&mut pc);
        predcoefs[j] = pc;
    }
    (predcoefs, ilsf)
}

// silk_A2NLSF (fixed-point forward A -> NLSF)

const LSF_COS_TAB_SZ_FIX: usize = 128;
const BIN_DIV_STEPS_A2NLSF_FIX: i32 = 3;
const MAX_ITERATIONS_A2NLSF_FIX: i32 = 16;
const SILK_INT16_MAX: i32 = 32767;

include!("silk_lsf_cos_tab.rs");

#[inline]
fn silk_rshift_round(a: i32, shift: i32) -> i32 {
    if shift == 1 {
        (a >> 1) + (a & 1)
    } else {
        ((a >> (shift - 1)) + 1) >> 1
    }
}

#[inline]
fn silk_smlaww(a32: i32, b32: i32, c32: i32) -> i32 {
    // silk_SMLAWW in its canonical 64-bit form: a32 + ((i64)b32 * c32 >> 16).
    (a32 as i64 + (((b32 as i64) * (c32 as i64)) >> 16)) as i32
}

#[inline]
fn silk_div32(a: i32, b: i32) -> i32 {
    a / b
}

#[inline]
fn silk_div32_16(a: i32, b: i32) -> i32 {
    a / b
}

/// `silk_bwexpander_32`: chirp (bandwidth expand) LPC coefficients in Q16 in place.
/// `chirp_Q16 += RSHIFT_ROUND(MUL(chirp, chirp-65536), 16)` with `MUL` a wrapping i32 multiply.
fn silk_bwexpander_32(ar: &mut [i32], d: usize, chirp_q16: i32) {
    let mut chirp = chirp_q16;
    let chirp_minus_one = chirp_q16 - 65536;
    for i in 0..d - 1 {
        ar[i] = (((chirp as i64) * (ar[i] as i64)) >> 16) as i32; // silk_SMULWW
        let mul = chirp.wrapping_mul(chirp_minus_one); // silk_MUL (i32 wrap)
        chirp = chirp.wrapping_add(silk_rshift_round(mul, 16));
    }
    ar[d - 1] = (((chirp as i64) * (ar[d - 1] as i64)) >> 16) as i32;
}

fn silk_a2nlsf_trans_poly(p: &mut [i32], dd: usize) {
    for k in 2..=dd {
        for n in (k + 1..=dd).rev() {
            p[n - 2] -= p[n];
        }
        p[k - 2] -= p[k] << 1;
    }
}

fn silk_a2nlsf_eval_poly(p: &[i32], x: i32, dd: usize) -> i32 {
    let x_q16 = x << 4;
    let mut y32 = p[dd];
    for n in (0..dd).rev() {
        y32 = silk_smlaww(p[n], y32, x_q16);
    }
    y32
}

fn silk_a2nlsf_init(a_q16: &[i32], p: &mut [i32], q: &mut [i32], dd: usize) {
    p[dd] = 1 << 16;
    q[dd] = 1 << 16;
    for k in 0..dd {
        p[k] = -a_q16[dd - k - 1] - a_q16[dd + k];
        q[k] = -a_q16[dd - k - 1] + a_q16[dd + k];
    }
    for k in (1..=dd).rev() {
        p[k - 1] -= p[k];
        q[k - 1] += q[k];
    }
    silk_a2nlsf_trans_poly(p, dd);
    silk_a2nlsf_trans_poly(q, dd);
}

/// `silk_A2NLSF`: monic whitening coefficients (Q16) -> NLSF (Q15, 0..2^15-1). Mutates `a_q16`
/// (bandwidth expansion on non-convergence). `d` is the filter order (even).
fn silk_a2nlsf(nlsf: &mut [i32], a_q16: &mut [i32], d: usize) {
    let dd = d >> 1;
    let mut p = vec![0i32; dd + 1];
    let mut q = vec![0i32; dd + 1];

    silk_a2nlsf_init(a_q16, &mut p, &mut q, dd);

    // p_sel: false => P, true => Q
    let mut use_q = false;
    let mut xlo = SILK_LSF_COS_TAB_FIX_Q12[0];
    let mut ylo = silk_a2nlsf_eval_poly(if use_q { &q } else { &p }, xlo, dd);

    let mut root_ix;
    if ylo < 0 {
        nlsf[0] = 0;
        use_q = true;
        ylo = silk_a2nlsf_eval_poly(&q, xlo, dd);
        root_ix = 1;
    } else {
        root_ix = 0;
    }
    let mut k = 1usize;
    let mut i = 0i32;
    let mut thr = 0i32;
    loop {
        let xhi = SILK_LSF_COS_TAB_FIX_Q12[k];
        let mut yhi = silk_a2nlsf_eval_poly(if use_q { &q } else { &p }, xhi, dd);

        if (ylo <= 0 && yhi >= thr) || (ylo >= 0 && yhi <= -thr) {
            thr = if yhi == 0 { 1 } else { 0 };
            let mut xlo_l = xlo;
            let mut ylo_l = ylo;
            let mut xhi_l = xhi;
            let mut ffrac = -256i32;
            for m in 0..BIN_DIV_STEPS_A2NLSF_FIX {
                let xmid = silk_rshift_round(xlo_l + xhi_l, 1);
                let ymid = silk_a2nlsf_eval_poly(if use_q { &q } else { &p }, xmid, dd);
                if (ylo_l <= 0 && ymid >= 0) || (ylo_l >= 0 && ymid <= 0) {
                    xhi_l = xmid;
                    yhi = ymid;
                } else {
                    xlo_l = xmid;
                    ylo_l = ymid;
                    ffrac += 128 >> m;
                }
            }
            if ylo_l.abs() < 65536 {
                let den = ylo_l - yhi;
                let nom = (ylo_l << (8 - BIN_DIV_STEPS_A2NLSF_FIX)) + (den >> 1);
                if den != 0 {
                    ffrac += silk_div32(nom, den);
                }
            } else {
                ffrac += silk_div32(ylo_l, (ylo_l - yhi) >> (8 - BIN_DIV_STEPS_A2NLSF_FIX));
            }
            nlsf[root_ix] = ((k as i32) << 8).wrapping_add(ffrac).min(SILK_INT16_MAX);

            root_ix += 1;
            if root_ix >= d {
                break;
            }
            use_q = (root_ix & 1) != 0;
            xlo = SILK_LSF_COS_TAB_FIX_Q12[k - 1];
            ylo = (1 - (root_ix as i32 & 2)) << 12;
        } else {
            k += 1;
            xlo = xhi;
            ylo = yhi;
            thr = 0;
            if k > LSF_COS_TAB_SZ_FIX {
                i += 1;
                if i > MAX_ITERATIONS_A2NLSF_FIX {
                    nlsf[0] = silk_div32_16(1 << 15, d as i32 + 1);
                    for kk in 1..d {
                        nlsf[kk] = nlsf[kk - 1] + nlsf[0];
                    }
                    return;
                }
                silk_bwexpander_32(a_q16, d, 65536 - (1 << i));
                silk_a2nlsf_init(a_q16, &mut p, &mut q, dd);
                use_q = false;
                xlo = SILK_LSF_COS_TAB_FIX_Q12[0];
                ylo = silk_a2nlsf_eval_poly(&p, xlo, dd);
                if ylo < 0 {
                    nlsf[0] = 0;
                    use_q = true;
                    ylo = silk_a2nlsf_eval_poly(&q, xlo, dd);
                    root_ix = 1;
                } else {
                    root_ix = 0;
                }
                k = 1;
            }
        }
    }
}

/// `smpl_A2NLSF_16`: post-BWE float `A[0..16]` (`A[0]=1`) -> analysis NLSF in radians (0..pi).
pub(crate) fn smpl_a2nlsf_16(a: &[f32]) -> [f32; SMPL_LPC_ORDER] {
    let mut a_q16 = [0i32; SMPL_LPC_ORDER];
    for i in 0..SMPL_LPC_ORDER {
        a_q16[i] = (-a[i + 1] * 65536.0).round() as i32;
    }
    let mut lsf_q15 = [0i32; SMPL_LPC_ORDER];
    silk_a2nlsf(&mut lsf_q15, &mut a_q16, SMPL_LPC_ORDER);
    let mut nlsf = [0.0f32; SMPL_LPC_ORDER];
    for i in 0..SMPL_LPC_ORDER {
        nlsf[i] = (lsf_q15[i] as f32) / 32768.0 * SMPL_PI;
    }
    nlsf
}

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

    fn fvec(v: &Value) -> Vec<f32> {
        v.as_array()
            .unwrap()
            .iter()
            .map(|x| x.as_f64().unwrap() as f32)
            .collect()
    }

    // The forward A2NLSF (silk_A2NLSF, fixed-point) must reproduce the reference `lsf` from the
    // post-BWE `A` for every captured record. Both use the same fixed-point arithmetic, so this is
    // exact (NLSF in radians within float rounding of the Q15->radians scale).
    #[test]
    fn a2nlsf_matches_c() {
        let recs: Value = serde_json::from_str(include_str!("testdata/lsf_quant_io.json")).unwrap();
        let arr = recs.as_array().unwrap();
        let mut worst = 0.0f32;
        for (n, r) in arr.iter().enumerate() {
            let a = fvec(&r["A"]);
            let want = fvec(&r["lsf"]);
            let got = smpl_a2nlsf_16(&a);
            for k in 0..SMPL_LPC_ORDER {
                let d = (got[k] - want[k]).abs();
                worst = worst.max(d);
                assert!(
                    d < 1e-4,
                    "rec {n} nlsf[{k}] {:.7} != ref {:.7} (d={d:.2e})",
                    got[k],
                    want[k]
                );
            }
        }
        eprintln!("a2nlsf worst abs error = {worst:.3e}");
    }

    // The full front-end (re-window from the raw lpcbuf, then FFT-autocorr -> reflection -> A -> BWE)
    // must reproduce the reference post-BWE `A` from the raw lpcbuf. The FFT differs only in internal
    // rounding (portable mixed-radix vs PFFFT), so A matches to a tight float tolerance (documented,
    // like the decoder postfilters).
    #[test]
    fn front_end_a_matches_c() {
        let recs: Value = serde_json::from_str(include_str!("testdata/fe_dump.json")).unwrap();
        let arr = recs.as_array().unwrap();
        assert!(arr.len() >= 12, "need front-end vectors");
        let mut worst = 0.0f32;
        let mut worst_win = 0.0f32;
        for (n, r) in arr.iter().enumerate() {
            let lpcbuf = fvec(&r["lpcbuf"]);
            let numframe = r["numframe"].as_i64().unwrap();
            assert_eq!(lpcbuf.len(), SMPL_LPC_BUF_LEN);
            let mut buf = [0.0f32; SMPL_LPC_BUF_LEN];
            buf.copy_from_slice(&lpcbuf);
            // use_long_win = numframe < 2 (frames_per_packet-1 == 2)
            let win = smpl_window_lpc20(&buf, numframe < 2);
            // Cross-check the windowing exactly against the reference `windowed`.
            let want_win = fvec(&r["windowed"]);
            for k in 0..SMPL_LPC_BUF_LEN {
                worst_win = worst_win.max((win[k] - want_win[k]).abs());
            }
            let a = smpl_lpc_analyze_with_f2(&win, &mut new_lpc_fft_scratch()).0;
            let want_a = fvec(&r["A"]);
            let r0: f64 = r["R"][0].as_f64().unwrap();
            let mut rd = 0.0f32;
            for k in 0..=SMPL_LPC_ORDER {
                rd = rd.max((a[k] - want_a[k]).abs());
            }
            worst = worst.max(rd);
            // Frames with audible energy must match `A` tightly (FFT-internal rounding only). Below a
            // very low energy floor the autocorrelation is dominated by `reg` + numerical noise so the
            // LPC is ill-conditioned and the reconstructed frame is silent regardless of `A`; skip the
            // `A` assert there but still verify the autocorrelation `R` itself (the FFT+DCT output).
            if r0 > 1e-7 {
                assert!(
                    rd < 5e-3,
                    "rec {n} (nf {numframe}, R0={r0:.2e}) |dA|={rd:.2e} too large",
                );
            }
        }
        // The windowing is exact; the autocorrelation tracks the reference FFT/DCT to a tight tolerance.
        assert!(worst_win < 1e-6, "windowing |dwin|={worst_win:.2e}");
    }

    // DIAGNOSTIC: the wire round-trip; feed the captured `qi` (grid+stage2) + threaded prev_nlsf to
    // the decoder's NLSF reconstruction and compare to the captured `qlsf`. Proves grid/stage2 map
    // directly onto the decoder wire (grid=qi[0], cond centroid=16) and that the decoder rebuilds the
    // same envelope.
    #[test]
    fn decoder_reconstructs_c_qlsf() {
        use super::super::smpl_lsf_quant::{lsf_quant, lsf_quant_cond};
        use super::super::smpl_synth::{load_smpl_synth_tables, smpl_reconstruct_nlsf};
        let recs: Value = serde_json::from_str(include_str!("testdata/lsf_quant_io.json")).unwrap();
        let arr = recs.as_array().unwrap();
        let st = load_smpl_synth_tables();
        let mut prev_nlsf: Vec<f32> = Vec::new();
        let mut worst = 0.0f32;
        for (n, r) in arr.iter().enumerate() {
            let a = fvec(&r["A"]);
            let lsf = fvec(&r["lsf"]);
            let voiced = r["voiced"].as_i64().unwrap() as usize;
            let low_rate = r["lowRate"].as_i64().unwrap() as usize;
            let surv = r["surv"].as_i64().unwrap() as usize;
            let rd_w = r["RDw_adj"].as_f64().unwrap() as f32;
            let cond = r["cond_coding"].as_i64().unwrap() != 0;
            let prev_lsf = fvec(&r["prev_lsf"]);
            let res = if cond {
                lsf_quant_cond(&a, &lsf, &prev_lsf, voiced, low_rate, rd_w, surv)
            } else {
                lsf_quant(&a, &lsf, voiced, low_rate, rd_w, surv)
            };
            let grid = res.qi[0] as usize;
            let mut stage2 = [0i32; 16];
            stage2.copy_from_slice(&res.qi[1..=16]);
            let rec = smpl_reconstruct_nlsf(st, voiced, 0, grid, &stage2, &prev_nlsf);
            // The decoder reconstruction must match the captured qlsf (the synthesis envelope),
            // proving the grid/stage2 (incl. cond centroid 16) round-trip and the decoder rebuilds the
            // same envelope, so no synthesis-time stabilization step is needed in the decoder.
            let want = fvec(&r["qlsf"]);
            let rd = (0..SMPL_LPC_ORDER)
                .map(|k| (rec[k] - want[k]).abs())
                .fold(0.0f32, f32::max);
            assert!(
                rd < 1e-3,
                "rec {n} cond={cond} grid={grid} reconstruct vs qlsf {rd:.2e}"
            );
            worst = worst.max(rd);
            prev_nlsf = rec;
        }
        assert!(
            worst < 2e-3,
            "decoder reconstruct vs C qlsf worst {worst:.3e}"
        );
    }
}