//! Inverse Discrete Cosine Transform (IDCT) implementation.
//!
//! This module provides 8x8 IDCT-III transform used for JPEG decoding.
//! Uses the recursive splitting algorithm from jpegli for exact compatibility.
//!
//! SIMD optimization via the `wide` crate is always enabled.

#![allow(dead_code)]

use crate::foundation::consts::DCT_BLOCK_SIZE;

use wide::f32x8;

/// WcMultipliers constants for IDCT (same as DCT).
/// Generated by: 1.0 / (2 * cos((i + 0.5) * pi / N))
const WC4: [f32; 2] = [0.541196100146197, 1.3065629648763764];

const WC8: [f32; 4] = [
    0.5097955791041592,
    0.6013448869350453,
    0.8999762231364156,
    2.5629154477415055,
];

const SQRT2: f32 = 1.41421356237;

/// Transpose an 8x8 block (scalar version, used in tests).
#[cfg(test)]
#[inline]
fn transpose_8x8(input: &[f32; 64], output: &mut [f32; 64]) {
    for row in 0..8 {
        for col in 0..8 {
            output[col * 8 + row] = input[row * 8 + col];
        }
    }
}

// SIMD-optimized implementations (always available via `wide` crate)
mod simd {
    use super::*;
    use multiversed::multiversed;

    /// SIMD-optimized 8x8 transpose using wide's built-in transpose.
    /// Accelerated on AVX.
    #[multiversed]
    pub fn transpose_8x8_simd(input: &[f32; 64], output: &mut [f32; 64]) {
        // Load all 8 rows
        let rows = [
            f32x8::from(<[f32; 8]>::try_from(&input[0..8]).unwrap()),
            f32x8::from(<[f32; 8]>::try_from(&input[8..16]).unwrap()),
            f32x8::from(<[f32; 8]>::try_from(&input[16..24]).unwrap()),
            f32x8::from(<[f32; 8]>::try_from(&input[24..32]).unwrap()),
            f32x8::from(<[f32; 8]>::try_from(&input[32..40]).unwrap()),
            f32x8::from(<[f32; 8]>::try_from(&input[40..48]).unwrap()),
            f32x8::from(<[f32; 8]>::try_from(&input[48..56]).unwrap()),
            f32x8::from(<[f32; 8]>::try_from(&input[56..64]).unwrap()),
        ];

        // Use wide's built-in transpose (AVX-accelerated)
        let transposed = f32x8::transpose(rows);

        // Store results
        output[0..8].copy_from_slice(&transposed[0].to_array());
        output[8..16].copy_from_slice(&transposed[1].to_array());
        output[16..24].copy_from_slice(&transposed[2].to_array());
        output[24..32].copy_from_slice(&transposed[3].to_array());
        output[32..40].copy_from_slice(&transposed[4].to_array());
        output[40..48].copy_from_slice(&transposed[5].to_array());
        output[48..56].copy_from_slice(&transposed[6].to_array());
        output[56..64].copy_from_slice(&transposed[7].to_array());
    }

    /// Transpose 8 f32x8 vectors in-place (no load/store overhead).
    /// Each f32x8 represents a row; returns transposed rows.
    #[inline(always)]
    fn transpose_vec(rows: [f32x8; 8]) -> [f32x8; 8] {
        f32x8::transpose(rows)
    }

    /// Vectorized 1D IDCT on 8 rows in parallel.
    /// Input: 8 f32x8 vectors where each represents corresponding values from 8 rows.
    /// (i.e., cols[0] contains element 0 from each of 8 rows)
    /// Returns: transformed values in the same layout.
    #[inline(always)]
    fn idct_1d_vec(cols: [f32x8; 8]) -> [f32x8; 8] {
        let [c0, c1, c2, c3, c4, c5, c6, c7] = cols;

        // ForwardEvenOdd: de-interleave
        let t0 = c0;
        let t1 = c2;
        let t2 = c4;
        let t3 = c6;
        let t4 = c1;
        let t5 = c3;
        let t6 = c5;
        let t7 = c7;

        // IDCT1D<4> on even part [t0, t1, t2, t3]
        let e0 = t0;
        let e1 = t2;
        let o0 = t1;
        let o1 = t3;

        // IDCT1D<2> on [e0, e1]
        let e00 = e0 + e1;
        let e01 = e0 - e1;

        // BTranspose<2> and IDCT1D<2> on [o0, o1]
        let o1b = o1 + o0;
        let o0b = o0 * f32x8::splat(SQRT2);
        let o00 = o0b + o1b;
        let o01 = o0b - o1b;

        // MultiplyAndAdd<4>
        let wc4_0 = f32x8::splat(WC4[0]);
        let wc4_1 = f32x8::splat(WC4[1]);
        let prod0 = wc4_0 * o00;
        let prod1 = wc4_1 * o01;
        let r0 = e00 + prod0;
        let r3 = e00 - prod0;
        let r1 = e01 + prod1;
        let r2 = e01 - prod1;

        // BTranspose<4> on odd part [t4, t5, t6, t7]
        let t7b = t7 + t6;
        let t6b = t6 + t5;
        let t5b = t5 + t4;
        let t4b = t4 * f32x8::splat(SQRT2);

        // IDCT1D<4> on [t4b, t5b, t6b, t7b]
        let e0o = t4b;
        let e1o = t6b;
        let o0o = t5b;
        let o1o = t7b;

        let e00o = e0o + e1o;
        let e01o = e0o - e1o;

        let o1bo = o1o + o0o;
        let o0bo = o0o * f32x8::splat(SQRT2);
        let o00o = o0bo + o1bo;
        let o01o = o0bo - o1bo;

        let prod0o = wc4_0 * o00o;
        let prod1o = wc4_1 * o01o;
        let r0o = e00o + prod0o;
        let r3o = e00o - prod0o;
        let r1o = e01o + prod1o;
        let r2o = e01o - prod1o;

        // MultiplyAndAdd<8>
        let wc8 = [
            f32x8::splat(WC8[0]),
            f32x8::splat(WC8[1]),
            f32x8::splat(WC8[2]),
            f32x8::splat(WC8[3]),
        ];

        let prod8_0 = wc8[0] * r0o;
        let prod8_1 = wc8[1] * r1o;
        let prod8_2 = wc8[2] * r2o;
        let prod8_3 = wc8[3] * r3o;

        [
            r0 + prod8_0, // c0
            r1 + prod8_1, // c1
            r2 + prod8_2, // c2
            r3 + prod8_3, // c3
            r3 - prod8_3, // c4
            r2 - prod8_2, // c5
            r1 - prod8_1, // c6
            r0 - prod8_0, // c7
        ]
    }

    /// Fast 1D IDCT on a single row (8 values).
    /// Uses the AAN (Arai, Agui, Nakajima) scaled IDCT algorithm.
    #[allow(dead_code)]
    #[inline]
    pub fn idct1d_fast(input: &[f32], output: &mut [f32]) {
        // Scaled IDCT constants (AAN algorithm)
        // These are cos(k*pi/16) * sqrt(2) factors
        const A1: f32 = 0.707_106_77; // cos(4*pi/16) * sqrt(2) = 1/sqrt(2)
        const A2: f32 = 0.541_196_1; // cos(6*pi/16) * sqrt(2)
        const A3: f32 = 0.707_106_77; // cos(4*pi/16) * sqrt(2) = 1/sqrt(2)
        const A4: f32 = 1.306_562_96; // cos(2*pi/16) * sqrt(2)
        const A5: f32 = 0.382_683_43; // cos(6*pi/16)

        let v0 = input[0];
        let v1 = input[1];
        let v2 = input[2];
        let v3 = input[3];
        let v4 = input[4];
        let v5 = input[5];
        let v6 = input[6];
        let v7 = input[7];

        // Stage 1: Reorder and prepare
        let t0 = v0;
        let t1 = v4;
        let t2 = v2;
        let t3 = v6;
        let t4 = v1;
        let t5 = v5;
        let t6 = v3;
        let t7 = v7;

        // Stage 2: First butterfly operations
        let s0 = t0 + t1;
        let s1 = t0 - t1;
        let s2 = t2 * A1 - t3 * A1;
        let s3 = t2 * A1 + t3 * A1;

        // Stage 3: Even part completion
        let e0 = s0 + s3;
        let e1 = s1 + s2;
        let e2 = s1 - s2;
        let e3 = s0 - s3;

        // Stage 2: Odd part
        let p4 = t4 + t7;
        let p5 = t5 + t6;
        let p6 = t4 + t6;
        let p7 = t5 + t7;
        let z5 = (p6 - p7) * A5;

        let o4 = t4 * A2 + z5 + p4 * (-A2 - A5);
        let o5 = t5 * A3 + p5 * (-A3);
        let o6 = t6 * A4 + z5 + p6 * (-A4 + A5);
        let o7 = t7 * A1;

        let q4 = o4 + o7;
        let q5 = o5 + o6;
        let q6 = o5 - o6;
        let q7 = -o4 + o7;

        // Stage 4: Final combination
        output[0] = e0 + q4;
        output[1] = e1 + q5;
        output[2] = e2 + q6;
        output[3] = e3 + q7;
        output[4] = e3 - q7;
        output[5] = e2 - q6;
        output[6] = e1 - q5;
        output[7] = e0 - q4;
    }

    /// Full SIMD-optimized 2D IDCT with f32x8 chaining.
    /// Keeps data in SIMD registers throughout the pipeline to avoid gather/scatter overhead.
    /// Uses only 2 transposes instead of the previous 4 gather/scatter cycles.
    #[multiversed]
    pub fn inverse_dct_8x8_simd(input: &[f32; 64]) -> [f32; 64] {
        // Load input as rows (each f32x8 is one row of the 8x8 block)
        let rows = [
            f32x8::from(<[f32; 8]>::try_from(&input[0..8]).unwrap()),
            f32x8::from(<[f32; 8]>::try_from(&input[8..16]).unwrap()),
            f32x8::from(<[f32; 8]>::try_from(&input[16..24]).unwrap()),
            f32x8::from(<[f32; 8]>::try_from(&input[24..32]).unwrap()),
            f32x8::from(<[f32; 8]>::try_from(&input[32..40]).unwrap()),
            f32x8::from(<[f32; 8]>::try_from(&input[40..48]).unwrap()),
            f32x8::from(<[f32; 8]>::try_from(&input[48..56]).unwrap()),
            f32x8::from(<[f32; 8]>::try_from(&input[56..64]).unwrap()),
        ];

        // Transpose: rows[j] -> cols[i] where cols[i] = [row0[i], row1[i], ..., row7[i]]
        // This layout allows idct_1d_vec to process 8 rows in parallel
        let cols = transpose_vec(rows);

        // Row IDCT: processes all 8 rows in parallel
        // After this, cols[i] = [IDCT(row0)[i], IDCT(row1)[i], ...]
        let cols_after_row = idct_1d_vec(cols);

        // Transpose back to row layout
        // rows[j] = complete IDCT'd row j = [IDCT(row_j)[0], IDCT(row_j)[1], ...]
        // This is also the correct layout for column IDCT:
        // idct_1d_vec needs input[i] = coeff i from each parallel column
        // coeff i of column k = row_i[k], so input[i] = rows[i]
        let rows_for_col = transpose_vec(cols_after_row);

        // Column IDCT: processes all 8 columns in parallel using row layout
        let final_rows = idct_1d_vec(rows_for_col);

        // Apply 1/8 scaling and store (already in row layout for storage)
        let scale = f32x8::splat(1.0 / 8.0);
        let mut output = [0.0f32; 64];

        output[0..8].copy_from_slice(&(final_rows[0] * scale).to_array());
        output[8..16].copy_from_slice(&(final_rows[1] * scale).to_array());
        output[16..24].copy_from_slice(&(final_rows[2] * scale).to_array());
        output[24..32].copy_from_slice(&(final_rows[3] * scale).to_array());
        output[32..40].copy_from_slice(&(final_rows[4] * scale).to_array());
        output[40..48].copy_from_slice(&(final_rows[5] * scale).to_array());
        output[48..56].copy_from_slice(&(final_rows[6] * scale).to_array());
        output[56..64].copy_from_slice(&(final_rows[7] * scale).to_array());

        output
    }
}

/// ForwardEvenOdd: De-interleave even/odd parts (opposite of InverseEvenOdd)
/// out[i] = in[2*i] for i in 0..N/2
/// out[N/2+i] = in[2*i+1] for i in 0..N/2
#[cfg(test)]
#[inline]
fn forward_even_odd<const N: usize>(input: &[f32], output: &mut [f32]) {
    let half = N / 2;
    for i in 0..half {
        output[i] = input[2 * i];
        output[half + i] = input[2 * i + 1];
    }
}

/// BTranspose: Reverse of B transform
/// First cumulative sum from end to start, then coeff[0] *= sqrt2
#[cfg(test)]
#[inline]
fn b_transpose<const N: usize>(coeff: &mut [f32]) {
    // Cumulative sum from end to start
    for i in (1..N).rev() {
        coeff[i] += coeff[i - 1];
    }
    coeff[0] *= SQRT2;
}

/// MultiplyAndAdd: Combine even and odd parts with Wc multipliers
/// out[i] = in[i] + Wc[i] * in[N/2+i]
/// out[N-1-i] = in[i] - Wc[i] * in[N/2+i]
#[cfg(test)]
#[inline]
fn multiply_and_add_8(input: &[f32], output: &mut [f32]) {
    for i in 0..4 {
        let even = input[i];
        let odd = input[4 + i];
        let prod = WC8[i] * odd;
        output[i] = even + prod;
        output[7 - i] = even - prod;
    }
}

#[cfg(test)]
#[inline]
fn multiply_and_add_4(input: &[f32], output: &mut [f32]) {
    for i in 0..2 {
        let even = input[i];
        let odd = input[2 + i];
        let prod = WC4[i] * odd;
        output[i] = even + prod;
        output[3 - i] = even - prod;
    }
}

/// IDCT base case for N=2
#[cfg(test)]
#[inline]
fn idct1d_2(input: &[f32], output: &mut [f32]) {
    let in1 = input[0];
    let in2 = input[1];
    output[0] = in1 + in2;
    output[1] = in1 - in2;
}

/// IDCT for N=4 (recursive)
#[cfg(test)]
fn idct1d_4(input: &[f32], output: &mut [f32]) {
    let mut tmp = [0.0f32; 4];

    // ForwardEvenOdd<4>: de-interleave
    forward_even_odd::<4>(input, &mut tmp);

    // IDCT1D<2> on first half (in-place)
    let (first, second) = tmp.split_at_mut(2);
    let mut first_out = [0.0f32; 2];
    idct1d_2(first, &mut first_out);
    first[0] = first_out[0];
    first[1] = first_out[1];

    // BTranspose<2>: just coeff[0] *= sqrt2 (after sum which doesn't exist for N=2)
    second[1] += second[0];
    second[0] *= SQRT2;

    // IDCT1D<2> on second half
    let mut second_out = [0.0f32; 2];
    idct1d_2(second, &mut second_out);
    second[0] = second_out[0];
    second[1] = second_out[1];

    // Recombine tmp for MultiplyAndAdd
    let tmp_combined = [first[0], first[1], second[0], second[1]];

    // MultiplyAndAdd<4>
    multiply_and_add_4(&tmp_combined, output);
}

/// IDCT for N=8 (recursive)
#[cfg(test)]
fn idct1d_8(input: &[f32], output: &mut [f32]) {
    let mut tmp = [0.0f32; 8];

    // ForwardEvenOdd<8>: de-interleave
    forward_even_odd::<8>(input, &mut tmp);

    // IDCT1D<4> on first half
    let mut first_out = [0.0f32; 4];
    idct1d_4(&tmp[0..4], &mut first_out);
    tmp[0] = first_out[0];
    tmp[1] = first_out[1];
    tmp[2] = first_out[2];
    tmp[3] = first_out[3];

    // BTranspose<4> on second half
    b_transpose::<4>(&mut tmp[4..8]);

    // IDCT1D<4> on second half
    let mut second_out = [0.0f32; 4];
    idct1d_4(&tmp[4..8], &mut second_out);
    tmp[4] = second_out[0];
    tmp[5] = second_out[1];
    tmp[6] = second_out[2];
    tmp[7] = second_out[3];

    // MultiplyAndAdd<8>
    multiply_and_add_8(&tmp, output);
}

/// 1D IDCT on all 8 rows (no scaling - scaling handled in main function)
#[cfg(test)]
fn idct_rows(input: &[f32; 64], output: &mut [f32; 64]) {
    for row in 0..8 {
        let in_slice = &input[row * 8..(row + 1) * 8];
        let mut out_row = [0.0f32; 8];
        idct1d_8(in_slice, &mut out_row);
        for i in 0..8 {
            output[row * 8 + i] = out_row[i];
        }
    }
}

/// Check if all AC coefficients are zero (DC-only block)
#[inline]
fn is_dc_only(input: &[f32; DCT_BLOCK_SIZE]) -> bool {
    // Check if all coefficients after DC are zero
    // Use a simple threshold check for floating point
    input[1..].iter().all(|&v| v.abs() < 1e-10)
}

/// Performs an 8x8 inverse DCT on the input coefficients.
///
/// Uses SIMD optimization when the `simd` feature is enabled.
/// Includes 1/8 scaling to match forward DCT.
///
/// # Arguments
/// * `input` - 8x8 block of DCT coefficients (after dequantization)
///
/// # Returns
/// 8x8 block of pixel values (before level shift adjustment)
#[must_use]
pub fn inverse_dct_8x8(input: &[f32; DCT_BLOCK_SIZE]) -> [f32; DCT_BLOCK_SIZE] {
    // Fast path: DC-only block (all AC coefficients are zero)
    // This is very common in JPEGs and can be computed directly
    if is_dc_only(input) {
        // For DC-only input, all output pixels have the same value
        // DC value after IDCT = DC_coeff / 8 (due to 1/8 scaling)
        let dc_value = input[0] / 8.0;
        return [dc_value; DCT_BLOCK_SIZE];
    }

    // Use SIMD-optimized version (always available via wide crate)
    simd::inverse_dct_8x8_simd(input)
}

/// Performs inverse DCT with level shift and clamping to u8.
///
/// # Arguments
/// * `input` - 8x8 block of dequantized DCT coefficients
///
/// # Returns
/// 8x8 block of pixel values (0-255)
#[must_use]
pub fn inverse_dct_8x8_u8(input: &[f32; DCT_BLOCK_SIZE]) -> [u8; DCT_BLOCK_SIZE] {
    use wide::f32x8;

    let output = inverse_dct_8x8(input);
    let mut result = [0u8; DCT_BLOCK_SIZE];

    let level_shift = f32x8::splat(128.0);
    let zero = f32x8::splat(0.0);
    let max_val = f32x8::splat(255.0);

    // Process 8 values at a time
    for chunk in 0..8 {
        let k = chunk * 8;

        // Load contiguous values directly from slice
        let vals = f32x8::from(<[f32; 8]>::try_from(&output[k..k + 8]).unwrap());

        // Add level shift, round, and clamp
        let shifted = vals + level_shift;
        let rounded = shifted.round();
        let clamped = rounded.max(zero).min(max_val);

        let arr: [f32; 8] = clamped.into();
        result[k] = arr[0] as u8;
        result[k + 1] = arr[1] as u8;
        result[k + 2] = arr[2] as u8;
        result[k + 3] = arr[3] as u8;
        result[k + 4] = arr[4] as u8;
        result[k + 5] = arr[5] as u8;
        result[k + 6] = arr[6] as u8;
        result[k + 7] = arr[7] as u8;
    }

    result
}

/// Performs inverse DCT on multiple blocks.
pub fn inverse_dct_blocks(blocks: &[[f32; DCT_BLOCK_SIZE]]) -> Vec<[f32; DCT_BLOCK_SIZE]> {
    blocks.iter().map(inverse_dct_8x8).collect()
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::encode::dct::forward_dct_8x8;

    #[test]
    fn test_dct_idct_roundtrip() {
        // Test with a pattern
        let mut input = [0.0f32; DCT_BLOCK_SIZE];
        for i in 0..DCT_BLOCK_SIZE {
            input[i] = (i as f32 * 3.7).sin() * 100.0;
        }

        let dct = forward_dct_8x8(&input);
        let recovered = inverse_dct_8x8(&dct);

        for i in 0..DCT_BLOCK_SIZE {
            assert!(
                (input[i] - recovered[i]).abs() < 0.1,
                "Mismatch at {}: {} vs {}",
                i,
                input[i],
                recovered[i]
            );
        }
    }

    #[test]
    fn test_idct_dc_only() {
        // DC-only input should produce constant output
        let mut input = [0.0f32; DCT_BLOCK_SIZE];
        input[0] = 64.0; // DC coefficient only

        let output = inverse_dct_8x8(&input);

        // All values should be the same (or very close)
        let first = output[0];
        for (i, &v) in output.iter().enumerate() {
            assert!(
                (v - first).abs() < 0.01,
                "Value at {} differs: {} vs {}",
                i,
                v,
                first
            );
        }
    }

    #[test]
    fn test_idct_u8_clamping() {
        // Test that output is properly clamped to [0, 255]
        let mut input = [0.0f32; DCT_BLOCK_SIZE];
        input[0] = 1000.0; // Very large DC - should produce high values

        let output = inverse_dct_8x8_u8(&input);

        // All values must be in valid u8 range (implicitly true for u8)
        // With a large positive DC, most values should be high
        let max_val = output.iter().copied().max().unwrap();

        // The output should have high values due to large DC
        assert!(
            max_val > 200,
            "Large DC should produce high values, got max={}",
            max_val
        );

        // Test negative clamping - with large negative DC, min should be 0 (clamped)
        input[0] = -1000.0;
        let output_neg = inverse_dct_8x8_u8(&input);
        let min_val_neg = output_neg.iter().copied().min().unwrap();
        // Should be clamped to low values (near 0)
        assert!(
            min_val_neg < 10,
            "Large negative DC should produce low values clamped near 0, got min={}",
            min_val_neg
        );
    }

    #[test]
    fn test_roundtrip_random_patterns() {
        // Test several different patterns
        for seed in 0..10 {
            let mut input = [0.0f32; DCT_BLOCK_SIZE];
            for i in 0..DCT_BLOCK_SIZE {
                input[i] = ((i + seed * 7) as f32 * 1.23).sin() * 127.0;
            }

            let dct = forward_dct_8x8(&input);
            let recovered = inverse_dct_8x8(&dct);

            for i in 0..DCT_BLOCK_SIZE {
                assert!(
                    (input[i] - recovered[i]).abs() < 0.1,
                    "Pattern {}: Mismatch at {}: {} vs {}",
                    seed,
                    i,
                    input[i],
                    recovered[i]
                );
            }
        }
    }

    #[test]
    fn test_dc_value_preservation() {
        // A constant block should roundtrip with the correct DC
        let value = 50.0f32;
        let input = [value; DCT_BLOCK_SIZE];

        let dct = forward_dct_8x8(&input);
        let recovered = inverse_dct_8x8(&dct);

        for (i, &v) in recovered.iter().enumerate() {
            assert!(
                (v - value).abs() < 0.01,
                "DC preservation failed at {}: {} vs {}",
                i,
                v,
                value
            );
        }
    }

    #[test]
    fn test_inverse_dct_blocks() {
        // Create DC-only blocks (only first coefficient non-zero)
        let mut block1 = [0.0f32; DCT_BLOCK_SIZE];
        block1[0] = 64.0;
        let mut block2 = [0.0f32; DCT_BLOCK_SIZE];
        block2[0] = 32.0;
        let blocks = vec![block1, block2];

        let results = inverse_dct_blocks(&blocks);
        assert_eq!(results.len(), 2);

        // DC-only blocks should produce constant outputs
        let first_val = results[0][0];
        for &v in &results[0] {
            assert!((v - first_val).abs() < 0.01);
        }
    }

    #[test]
    fn test_inverse_dct_8x8_u8_negative() {
        // Test negative DC that would underflow
        let mut input = [0.0f32; DCT_BLOCK_SIZE];
        input[0] = -2000.0; // Very negative DC

        let output = inverse_dct_8x8_u8(&input);
        // All values should be clamped to 0
        for &v in &output {
            assert_eq!(v, 0);
        }
    }

    #[test]
    fn test_is_dc_only() {
        let dc_only = [
            100.0f32, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
        ];
        assert!(is_dc_only(&dc_only));

        let not_dc_only = [
            100.0f32, 50.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
        ];
        assert!(!is_dc_only(&not_dc_only));
    }

    #[test]
    fn test_transpose_8x8() {
        let mut input = [0.0f32; 64];
        for i in 0..64 {
            input[i] = i as f32;
        }
        let mut output = [0.0f32; 64];
        transpose_8x8(&input, &mut output);

        // Check that output[col*8 + row] = input[row*8 + col]
        for row in 0..8 {
            for col in 0..8 {
                assert_eq!(output[col * 8 + row], input[row * 8 + col]);
            }
        }
    }

    #[test]
    fn test_idct1d_2() {
        let input = [3.0f32, 1.0];
        let mut output = [0.0f32; 2];
        idct1d_2(&input, &mut output);
        assert_eq!(output[0], 4.0); // 3 + 1
        assert_eq!(output[1], 2.0); // 3 - 1
    }

    #[test]
    fn test_forward_even_odd() {
        let input = [0.0f32, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0];
        let mut output = [0.0f32; 8];
        forward_even_odd::<8>(&input, &mut output);
        // Even: 0, 2, 4, 6; Odd: 1, 3, 5, 7
        assert_eq!(output[0], 0.0);
        assert_eq!(output[1], 2.0);
        assert_eq!(output[2], 4.0);
        assert_eq!(output[3], 6.0);
        assert_eq!(output[4], 1.0);
        assert_eq!(output[5], 3.0);
        assert_eq!(output[6], 5.0);
        assert_eq!(output[7], 7.0);
    }

    #[test]
    fn test_b_transpose() {
        let mut coeff = [1.0f32, 2.0, 3.0, 4.0];
        b_transpose::<4>(&mut coeff);
        // After cumulative sum from end: [1, 2+1=3, 3+2=5, 4+3=7] -> [1, 3, 5, 7]
        // But wait, it's from end to start: [1+2, 2+3, 3+4, 4] = [3, 5, 7, 4]
        // Actually: for i in (1..N).rev() means i = 3,2,1
        // coeff[3] += coeff[2] -> 4+3=7
        // coeff[2] += coeff[1] -> 3+2=5
        // coeff[1] += coeff[0] -> 2+1=3
        // Then coeff[0] *= sqrt2 -> 1 * sqrt2 ≈ 1.414
        assert!((coeff[0] - SQRT2).abs() < 0.001);
        assert_eq!(coeff[1], 3.0);
        // coeff[2] depends on intermediate state
    }

    #[test]
    fn test_multiply_and_add_4() {
        let input = [10.0f32, 20.0, 5.0, 8.0]; // even: 10, 20; odd: 5, 8
        let mut output = [0.0f32; 4];
        multiply_and_add_4(&input, &mut output);
        // output[0] = 10 + WC4[0] * 5
        // output[3] = 10 - WC4[0] * 5
        // output[1] = 20 + WC4[1] * 8
        // output[2] = 20 - WC4[1] * 8
        assert!((output[0] - (10.0 + WC4[0] * 5.0)).abs() < 0.001);
        assert!((output[3] - (10.0 - WC4[0] * 5.0)).abs() < 0.001);
    }

    #[test]
    fn test_multiply_and_add_8() {
        let input = [1.0f32, 2.0, 3.0, 4.0, 0.5, 0.6, 0.7, 0.8];
        let mut output = [0.0f32; 8];
        multiply_and_add_8(&input, &mut output);
        // Check a few values
        assert!((output[0] - (1.0 + WC8[0] * 0.5)).abs() < 0.001);
        assert!((output[7] - (1.0 - WC8[0] * 0.5)).abs() < 0.001);
    }

    #[test]
    fn test_idct1d_4() {
        let input = [10.0f32, 5.0, 2.0, 1.0];
        let mut output = [0.0f32; 4];
        idct1d_4(&input, &mut output);
        // Just verify it doesn't panic and produces reasonable values
        for v in &output {
            assert!(v.is_finite());
        }
    }

    #[test]
    fn test_idct1d_8() {
        let input = [64.0f32, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
        let mut output = [0.0f32; 8];
        idct1d_8(&input, &mut output);
        // DC-only input should produce constant output
        let first = output[0];
        for v in &output {
            assert!((v - first).abs() < 0.01);
        }
    }

    #[test]
    fn test_idct_rows() {
        let mut input = [0.0f32; 64];
        input[0] = 64.0; // DC only
        let mut output = [0.0f32; 64];
        idct_rows(&input, &mut output);
        // First row should have uniform values
        let first = output[0];
        for i in 0..8 {
            assert!((output[i] - first).abs() < 0.01);
        }
    }

    #[test]
    fn test_simd_transpose_matches_scalar() {
        let mut input = [0.0f32; 64];
        for i in 0..64 {
            input[i] = (i as f32 * 1.5).sin() * 100.0;
        }

        let mut scalar_output = [0.0f32; 64];
        let mut simd_output = [0.0f32; 64];

        transpose_8x8(&input, &mut scalar_output);
        simd::transpose_8x8_simd(&input, &mut simd_output);

        for i in 0..64 {
            assert!(
                (scalar_output[i] - simd_output[i]).abs() < 0.001,
                "Mismatch at {}: {} vs {}",
                i,
                scalar_output[i],
                simd_output[i]
            );
        }
    }

    #[test]
    fn test_simd_idct_matches_scalar() {
        // Test with random input
        let mut input = [0.0f32; 64];
        for i in 0..64 {
            input[i] = ((i * 17) as f32).sin() * 50.0;
        }

        let simd_result = simd::inverse_dct_8x8_simd(&input);

        // Compare with scalar implementation
        let mut block0 = input;
        let mut block1 = [0.0f32; 64];

        transpose_8x8(&block0, &mut block1);
        idct_rows(&block1, &mut block0);
        transpose_8x8(&block0, &mut block1);
        let mut scalar_result = [0.0f32; 64];
        idct_rows(&block1, &mut scalar_result);
        for v in &mut scalar_result {
            *v /= 8.0;
        }

        for i in 0..64 {
            assert!(
                (scalar_result[i] - simd_result[i]).abs() < 0.1,
                "SIMD/scalar mismatch at {}: {} vs {}",
                i,
                scalar_result[i],
                simd_result[i]
            );
        }
    }
}
