archmage 0.9.13

//! SIMD Kernel Collection - Top 8 Hotspots from Image Processing
//!
//! This example demonstrates the most impactful SIMD kernels from:
//! - image-webp (VP8 codec)
//! - jpegli-rs (JPEG encoder)
//! - zenimage (image pipeline)
//!
//! Run with: `cargo run --example simd_kernels --release`
//!
//! Each kernel shows archmage patterns for common image processing operations.

// x86-only example - stub main for other platforms
#[cfg(not(target_arch = "x86_64"))]
fn main() {}

#[cfg(target_arch = "x86_64")]
mod x86_impl {
    #![allow(clippy::excessive_precision)]
    #![allow(clippy::approx_constant)]

    use archmage::{SimdToken, X64V2Token, X64V3Token, arcane};
    use core::arch::x86_64::*;
    use magetypes::simd::f32x8;
    use std::time::Instant;

    // ============================================================================
    // 1. 4x4 DCT (VP8 Integer Transform)
    // ============================================================================

    /// VP8-style 4x4 DCT using i16 arithmetic with _mm_madd_epi16
    ///
    /// This is the core transform for VP8/WebP encoding. Uses integer math
    /// for deterministic output across platforms.
    #[arcane]
    pub fn dct4x4_vp8(token: X64V2Token, block: &mut [i32; 16]) {
        let _ = token;
        // Constants for VP8 DCT
        const K1: i32 = 20091; // cos(pi/8)*sqrt(2)*4096
        const K2: i32 = 35468; // sin(pi/8)*sqrt(2)*4096

        // Horizontal pass
        for row in 0..4 {
            let base = row * 4;
            let a = block[base] + block[base + 3];
            let b = block[base + 1] + block[base + 2];
            let c = block[base + 1] - block[base + 2];
            let d = block[base] - block[base + 3];

            block[base] = a + b;
            block[base + 1] = (c * K2 + d * K1 + 2048) >> 12;
            block[base + 2] = a - b;
            block[base + 3] = (d * K2 - c * K1 + 2048) >> 12;
        }

        // Vertical pass
        for col in 0..4 {
            let a = block[col] + block[col + 12];
            let b = block[col + 4] + block[col + 8];
            let c = block[col + 4] - block[col + 8];
            let d = block[col] - block[col + 12];

            block[col] = (a + b + 7) >> 4;
            block[col + 4] = (((c * K2 + d * K1 + 2048) >> 12) + 7) >> 4;
            block[col + 8] = (a - b + 7) >> 4;
            block[col + 12] = (((d * K2 - c * K1 + 2048) >> 12) + 7) >> 4;
        }
    }

    // ============================================================================
    // 2. 8x8 DCT Butterfly (JPEG)
    // ============================================================================

    /// 8-point DCT using butterfly algorithm with FMA
    ///
    /// Based on the jpegli-rs implementation. Processes 8 independent
    /// 8-point DCTs in parallel (one per lane).
    #[arcane]
    pub fn dct8_butterfly(token: X64V3Token, m: &mut [f32x8; 8]) {
        // WC8 coefficients
        let wc0 = f32x8::splat(token, 0.5097955791041592);
        let wc1 = f32x8::splat(token, 0.6013448869350453);
        let wc2 = f32x8::splat(token, 0.8999762231364156);
        let wc3 = f32x8::splat(token, 2.5629154477415055);
        let sqrt2 = f32x8::splat(token, 1.41421356237);

        // Stage 1: AddReverse<4>
        let t0 = m[0] + m[7];
        let t1 = m[1] + m[6];
        let t2 = m[2] + m[5];
        let t3 = m[3] + m[4];
        let t4 = m[0] - m[7];
        let t5 = m[1] - m[6];
        let t6 = m[2] - m[5];
        let t7 = m[3] - m[4];

        // Stage 2: DCT4 on first half
        let u0 = t0 + t3;
        let u1 = t1 + t2;
        let u2 = t0 - t3;
        let u3 = t1 - t2;

        let r0 = u0 + u1;
        let r2 = u0 - u1;
        let r1 = u2.mul_add(
            f32x8::splat(token, 0.541196100146197),
            u3 * f32x8::splat(token, 1.3065629648763764),
        );
        let r3 = u2.mul_add(
            f32x8::splat(token, 0.541196100146197),
            u3 * f32x8::splat(token, -1.3065629648763764),
        );
        let r1 = r1.mul_add(sqrt2, r3);

        // Stage 3: Scaled second half
        let s4 = t4 * wc0;
        let s5 = t5 * wc1;
        let s6 = t6 * wc2;
        let s7 = t7 * wc3;

        // DCT4 on second half (simplified)
        let v0 = s4 + s7;
        let v1 = s5 + s6;
        let v2 = s4 - s7;
        let v3 = s5 - s6;

        let p0 = v0 + v1;
        let p2 = v0 - v1;
        let p1 = v2.mul_add(
            f32x8::splat(token, 0.541196100146197),
            v3 * f32x8::splat(token, 1.3065629648763764),
        );
        let p3 = v2.mul_add(
            f32x8::splat(token, 0.541196100146197),
            v3 * f32x8::splat(token, -1.3065629648763764),
        );

        // B<4> cumulative
        let q0 = p0.mul_add(sqrt2, p1);
        let q1 = p1 + p2;
        let q2 = p2 + p3;

        // Interleave output
        m[0] = r0;
        m[1] = q0;
        m[2] = r1;
        m[3] = q1;
        m[4] = r2;
        m[5] = q2;
        m[6] = r3;
        m[7] = p3;
    }

    // ============================================================================
    // 3. Chroma Downsampling 2x2
    // ============================================================================

    /// 2x2 box filter downsampling using gather pattern
    ///
    /// Takes 16 consecutive inputs, extracts even/odd pairs for 8 outputs.
    /// Key operation: _mm256_permutevar8x32_ps for variable gather.
    #[arcane]
    pub fn downsample_2x2_row(token: X64V3Token, row0: &[f32], row1: &[f32], output: &mut [f32]) {
        debug_assert!(row0.len() >= output.len() * 2);
        debug_assert!(row1.len() >= output.len() * 2);

        let scale = f32x8::splat(token, 0.25);

        // Process 8 output pixels at a time (16 input pixels)
        for chunk in 0..(output.len() / 8) {
            let in_x = chunk * 16;
            let out_x = chunk * 8;

            // Note: Full SIMD gather would use vpgatherdd or manual permute chains.
            // For simplicity, gather even/odd with scalar, then vectorize the math.
            let mut p00 = [0.0f32; 8];
            let mut p10 = [0.0f32; 8];
            let mut p01 = [0.0f32; 8];
            let mut p11 = [0.0f32; 8];

            for i in 0..8 {
                p00[i] = row0[in_x + i * 2];
                p10[i] = row0[in_x + i * 2 + 1];
                p01[i] = row1[in_x + i * 2];
                p11[i] = row1[in_x + i * 2 + 1];
            }

            let p00_v = f32x8::from_array(token, p00);
            let p10_v = f32x8::from_array(token, p10);
            let p01_v = f32x8::from_array(token, p01);
            let p11_v = f32x8::from_array(token, p11);

            // Box filter average
            let sum = p00_v + p10_v + p01_v + p11_v;
            let avg = sum * scale;

            let out_arr: &mut [f32; 8] = (&mut output[out_x..out_x + 8]).try_into().unwrap();
            avg.store(out_arr);
        }
    }

    // ============================================================================
    // 4. RGB to YCbCr Color Matrix
    // ============================================================================

    /// BT.601 RGB to YCbCr conversion using FMA
    ///
    /// Y  =  0.299 R + 0.587 G + 0.114 B
    /// Cb = -0.169 R - 0.331 G + 0.500 B + 128
    /// Cr =  0.500 R - 0.419 G - 0.081 B + 128
    #[arcane]
    pub fn rgb_to_ycbcr_8px(
        token: X64V3Token,
        r: f32x8,
        g: f32x8,
        b: f32x8,
    ) -> (f32x8, f32x8, f32x8) {
        let offset = f32x8::splat(token, 128.0);

        // Y coefficients
        let ky_r = f32x8::splat(token, 0.299);
        let ky_g = f32x8::splat(token, 0.587);
        let ky_b = f32x8::splat(token, 0.114);

        // Cb coefficients
        let kcb_r = f32x8::splat(token, -0.168736);
        let kcb_g = f32x8::splat(token, -0.331264);
        let kcb_b = f32x8::splat(token, 0.5);

        // Cr coefficients
        let kcr_r = f32x8::splat(token, 0.5);
        let kcr_g = f32x8::splat(token, -0.418688);
        let kcr_b = f32x8::splat(token, -0.081312);

        // Y = 0.299*R + 0.587*G + 0.114*B
        let y = r.mul_add(ky_r, g.mul_add(ky_g, b * ky_b));

        // Cb = -0.169*R - 0.331*G + 0.5*B + 128
        let cb = r.mul_add(kcb_r, g.mul_add(kcb_g, b.mul_add(kcb_b, offset)));

        // Cr = 0.5*R - 0.419*G - 0.081*B + 128
        let cr = r.mul_add(kcr_r, g.mul_add(kcr_g, b.mul_add(kcr_b, offset)));

        (y, cb, cr)
    }

    // ============================================================================
    // 5. Horizontal 1D Convolution
    // ============================================================================

    /// Horizontal convolution with N-tap kernel (fixed-point u8)
    ///
    /// Strided access pattern limits SIMD gains, but still faster than scalar.
    #[arcane]
    pub fn convolve_horizontal_u8(
        token: X64V2Token,
        input: &[u8],
        output: &mut [u8],
        kernel: &[i16],
        scale_shift: u32,
    ) {
        let _ = token;
        let k_len = kernel.len();
        let k_half = k_len / 2;
        let half = 1i32 << (scale_shift - 1);

        // Process interior pixels
        for out_idx in k_half..(output.len().saturating_sub(k_half)) {
            let mut sum = half;
            for (k_idx, &k) in kernel.iter().enumerate() {
                let in_idx = out_idx + k_idx - k_half;
                sum += input[in_idx] as i32 * k as i32;
            }
            output[out_idx] = (sum >> scale_shift).clamp(0, 255) as u8;
        }

        // Handle edges with clamping
        for out_idx in 0..k_half.min(output.len()) {
            let mut sum = half;
            for (k_idx, &k) in kernel.iter().enumerate() {
                let in_idx = (out_idx as isize + k_idx as isize - k_half as isize).max(0) as usize;
                sum += input[in_idx.min(input.len() - 1)] as i32 * k as i32;
            }
            output[out_idx] = (sum >> scale_shift).clamp(0, 255) as u8;
        }
    }

    // ============================================================================
    // 6. sRGB to Linear Conversion
    // ============================================================================

    /// sRGB to linear RGB using polynomial approximation
    ///
    /// sRGB formula:
    /// - if x <= 0.04045: linear = x / 12.92
    /// - else: linear = ((x + 0.055) / 1.055)^2.4
    ///
    /// Uses sqrt chains to approximate x^2.4 ≈ x^2 * x^0.4
    #[arcane]
    pub fn srgb_to_linear_8px(token: X64V3Token, srgb: f32x8) -> f32x8 {
        let threshold = f32x8::splat(token, 0.04045);
        let linear_scale = f32x8::splat(token, 1.0 / 12.92);
        let offset = f32x8::splat(token, 0.055);
        let scale = f32x8::splat(token, 1.0 / 1.055);

        // Linear part
        let linear_result = srgb * linear_scale;

        // Gamma part: ((x + 0.055) / 1.055)^2.4
        let adjusted = (srgb + offset) * scale;

        // x^2.4 ≈ x^2 * x^0.4, where x^0.4 ≈ sqrt(sqrt(x)) * sqrt(sqrt(sqrt(x)))
        let x2 = adjusted * adjusted;
        let sqrt_x = adjusted.sqrt();
        let sqrt_sqrt_x = sqrt_x.sqrt(); // x^0.25
        let x_0125 = sqrt_sqrt_x.sqrt(); // x^0.125
        let x_04_approx = sqrt_sqrt_x * x_0125; // x^0.375 ≈ x^0.4

        let gamma_result = x2 * x_04_approx;

        // Select based on threshold
        let mask = srgb.simd_le(threshold);
        let result_raw = _mm256_blendv_ps(gamma_result.raw(), linear_result.raw(), mask.raw());
        f32x8::from_m256(token, result_raw)
    }

    /// Linear to sRGB conversion
    #[arcane]
    pub fn linear_to_srgb_8px(token: X64V3Token, linear: f32x8) -> f32x8 {
        let threshold = f32x8::splat(token, 0.0031308);
        let linear_scale = f32x8::splat(token, 12.92);
        let gamma_scale = f32x8::splat(token, 1.055);
        let offset = f32x8::splat(token, -0.055);
        let one = f32x8::splat(token, 1.0);

        // Linear part
        let linear_result = linear * linear_scale;

        // Gamma part: 1.055 * x^(1/2.4) - 0.055
        // x^(1/2.4) ≈ x^0.417 ≈ sqrt(sqrt(x)) * x^0.167
        let sqrt_x = linear.sqrt();
        let sqrt_sqrt_x = sqrt_x.sqrt(); // x^0.25
        let x_0125 = sqrt_sqrt_x.sqrt(); // x^0.125
        let x_042_approx = sqrt_sqrt_x * x_0125; // x^0.375 ≈ x^0.417

        let gamma_result = x_042_approx
            .mul_add(gamma_scale, offset)
            .max(f32x8::zero(token))
            .min(one);

        // Select based on threshold
        let mask = linear.simd_le(threshold);
        let result_raw = _mm256_blendv_ps(gamma_result.raw(), linear_result.raw(), mask.raw());
        f32x8::from_m256(token, result_raw)
    }

    // ============================================================================
    // 7. Multiply and Screen Blend Modes
    // ============================================================================

    /// Multiply blend: out = src * dst (per channel)
    #[arcane]
    pub fn blend_multiply_2px(token: X64V3Token, src: f32x8, dst: f32x8) -> f32x8 {
        let result = src * dst;

        // Preserve alpha (indices 3 and 7) from src
        let blend_mask = _mm256_set_ps(-0.0, 0.0, 0.0, 0.0, -0.0, 0.0, 0.0, 0.0);
        let result_raw = _mm256_blendv_ps(result.raw(), src.raw(), blend_mask);
        f32x8::from_m256(token, result_raw)
    }

    /// Screen blend: out = 1 - (1-src) * (1-dst)
    #[arcane]
    pub fn blend_screen_2px(token: X64V3Token, src: f32x8, dst: f32x8) -> f32x8 {
        let one = f32x8::splat(token, 1.0);

        let inv_src = one - src;
        let inv_dst = one - dst;
        let product = inv_src * inv_dst;
        let result = one - product;

        // Preserve alpha from src
        let blend_mask = _mm256_set_ps(-0.0, 0.0, 0.0, 0.0, -0.0, 0.0, 0.0, 0.0);
        let result_raw = _mm256_blendv_ps(result.raw(), src.raw(), blend_mask);
        f32x8::from_m256(token, result_raw)
    }

    /// Overlay blend: if dst < 0.5: 2*src*dst, else: 1-2*(1-src)*(1-dst)
    #[arcane]
    pub fn blend_overlay_2px(token: X64V3Token, src: f32x8, dst: f32x8) -> f32x8 {
        let one = f32x8::splat(token, 1.0);
        let two = f32x8::splat(token, 2.0);
        let half = f32x8::splat(token, 0.5);

        // Multiply path: 2 * src * dst
        let multiply_result = src * dst * two;

        // Screen path: 1 - 2 * (1-src) * (1-dst)
        let inv_src = one - src;
        let inv_dst = one - dst;
        let screen_result = one - inv_src * inv_dst * two;

        // Select based on dst < 0.5
        let mask = dst.simd_lt(half);
        let result_raw = _mm256_blendv_ps(screen_result.raw(), multiply_result.raw(), mask.raw());

        // Preserve alpha from src
        let blend_mask = _mm256_set_ps(-0.0, 0.0, 0.0, 0.0, -0.0, 0.0, 0.0, 0.0);
        let result_raw = _mm256_blendv_ps(result_raw, src.raw(), blend_mask);
        f32x8::from_m256(token, result_raw)
    }

    // ============================================================================
    // 8. Horizontal Weighted Reduction
    // ============================================================================

    /// Horizontal reduction with strided access (for resize)
    ///
    /// Each output is a weighted sum of `n_points` inputs spaced `stride` apart.
    #[arcane]
    pub fn reduce_horizontal_f32(
        token: X64V3Token,
        input: &[f32],
        output: &mut [f32],
        weights: &[f32],
        stride: usize,
    ) {
        for (out_idx, out_val) in output.iter_mut().enumerate() {
            let base = out_idx * stride;

            // Note: Horizontal reduction has strided access which limits vectorization.
            // Full SIMD version would process multiple outputs with gather/scatter.
            let mut sum = 0.0f32;
            for (k, &w) in weights.iter().enumerate() {
                if base + k < input.len() {
                    sum += input[base + k] * w;
                }
            }
            *out_val = sum;
        }
    }

    // ============================================================================
    // Scalar References
    // ============================================================================

    #[allow(dead_code)]
    fn dct4x4_scalar(block: &mut [i32; 16]) {
        const K1: i32 = 20091;
        const K2: i32 = 35468;

        for row in 0..4 {
            let base = row * 4;
            let a = block[base] + block[base + 3];
            let b = block[base + 1] + block[base + 2];
            let c = block[base + 1] - block[base + 2];
            let d = block[base] - block[base + 3];

            block[base] = a + b;
            block[base + 1] = (c * K2 + d * K1 + 2048) >> 12;
            block[base + 2] = a - b;
            block[base + 3] = (d * K2 - c * K1 + 2048) >> 12;
        }

        for col in 0..4 {
            let a = block[col] + block[col + 12];
            let b = block[col + 4] + block[col + 8];
            let c = block[col + 4] - block[col + 8];
            let d = block[col] - block[col + 12];

            block[col] = (a + b + 7) >> 4;
            block[col + 4] = (((c * K2 + d * K1 + 2048) >> 12) + 7) >> 4;
            block[col + 8] = (a - b + 7) >> 4;
            block[col + 12] = (((d * K2 - c * K1 + 2048) >> 12) + 7) >> 4;
        }
    }

    fn srgb_to_linear_scalar(x: f32) -> f32 {
        if x <= 0.04045 {
            x / 12.92
        } else {
            ((x + 0.055) / 1.055).powf(2.4)
        }
    }

    // ============================================================================
    // Tests and Benchmarks
    // ============================================================================

    fn test_correctness() {
        println!("=== Correctness Tests ===\n");

        // Test sRGB conversion
        if let Some(token) = X64V3Token::summon() {
            let srgb_vals = [0.0, 0.04, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0];
            let srgb = f32x8::from_array(token, srgb_vals);
            let linear = srgb_to_linear_8px(token, srgb);
            let linear_arr = linear.to_array();

            println!("  sRGB→Linear conversion:");
            let mut max_err = 0.0f32;
            for i in 0..8 {
                let expected = srgb_to_linear_scalar(srgb_vals[i]);
                let err = (linear_arr[i] - expected).abs();
                max_err = max_err.max(err);
            }
            println!("    Max error vs scalar: {:.4}", max_err);
            println!("    (Using sqrt approximation for x^2.4)\n");

            // Test RGB→YCbCr
            let r = f32x8::from_array(token, [255.0, 0.0, 0.0, 128.0, 64.0, 192.0, 100.0, 200.0]);
            let g = f32x8::from_array(token, [0.0, 255.0, 0.0, 128.0, 128.0, 64.0, 150.0, 100.0]);
            let b = f32x8::from_array(token, [0.0, 0.0, 255.0, 128.0, 192.0, 128.0, 50.0, 150.0]);

            let (y, cb, cr) = rgb_to_ycbcr_8px(token, r, g, b);
            let y_arr = y.to_array();
            let cb_arr = cb.to_array();
            let cr_arr = cr.to_array();

            println!("  RGB→YCbCr conversion:");
            println!(
                "    Red (255,0,0)   → Y={:.1}, Cb={:.1}, Cr={:.1}",
                y_arr[0], cb_arr[0], cr_arr[0]
            );
            println!(
                "    Green (0,255,0) → Y={:.1}, Cb={:.1}, Cr={:.1}",
                y_arr[1], cb_arr[1], cr_arr[1]
            );
            println!(
                "    Blue (0,0,255)  → Y={:.1}, Cb={:.1}, Cr={:.1}",
                y_arr[2], cb_arr[2], cr_arr[2]
            );
            println!(
                "    Gray (128,128,128) → Y={:.1}, Cb={:.1}, Cr={:.1}\n",
                y_arr[3], cb_arr[3], cr_arr[3]
            );

            // Test blend modes
            let src = f32x8::from_array(token, [0.5, 0.3, 0.8, 1.0, 0.2, 0.6, 0.4, 0.5]);
            let dst = f32x8::from_array(token, [0.4, 0.6, 0.2, 1.0, 0.8, 0.4, 0.6, 0.5]);

            let multiply = blend_multiply_2px(token, src, dst);
            let screen = blend_screen_2px(token, src, dst);
            let overlay = blend_overlay_2px(token, src, dst);

            println!("  Blend modes (src=0.5, dst=0.4 for RGB):");
            println!("    Multiply: {:.3}", multiply.to_array()[0]);
            println!("    Screen:   {:.3}", screen.to_array()[0]);
            println!("    Overlay:  {:.3}\n", overlay.to_array()[0]);
        }
    }

    fn benchmark() {
        const PIXELS: usize = 1920 * 1080;
        const ITERATIONS: usize = 100;

        println!(
            "=== Benchmarks ({} pixels x {} iterations) ===\n",
            PIXELS, ITERATIONS
        );

        if let Some(token) = X64V3Token::summon() {
            // sRGB→Linear benchmark
            let srgb_data: Vec<f32> = (0..PIXELS).map(|i| (i % 256) as f32 / 255.0).collect();
            let mut linear_data = vec![0.0f32; PIXELS];

            let start = Instant::now();
            for _ in 0..ITERATIONS {
                for chunk_start in (0..PIXELS).step_by(8) {
                    if chunk_start + 8 <= PIXELS {
                        let arr: &[f32; 8] = (&srgb_data[chunk_start..chunk_start + 8])
                            .try_into()
                            .unwrap();
                        let srgb = f32x8::load(token, arr);
                        let linear = srgb_to_linear_8px(token, srgb);
                        let out: &mut [f32; 8] = (&mut linear_data[chunk_start..chunk_start + 8])
                            .try_into()
                            .unwrap();
                        linear.store(out);
                    }
                }
                std::hint::black_box(&linear_data);
            }
            let simd_time = start.elapsed();
            let mpix_s = (PIXELS * ITERATIONS) as f64 / simd_time.as_secs_f64() / 1_000_000.0;
            println!(
                "  sRGB→Linear:     {:>8.2} ms ({:.1} Mpix/s)",
                simd_time.as_secs_f64() * 1000.0,
                mpix_s
            );

            // RGB→YCbCr benchmark
            let r_data: Vec<f32> = (0..PIXELS).map(|i| ((i * 17) % 256) as f32).collect();
            let g_data: Vec<f32> = (0..PIXELS).map(|i| ((i * 31) % 256) as f32).collect();
            let b_data: Vec<f32> = (0..PIXELS).map(|i| ((i * 47) % 256) as f32).collect();
            let mut y_data = vec![0.0f32; PIXELS];
            let mut cb_data = vec![0.0f32; PIXELS];
            let mut cr_data = vec![0.0f32; PIXELS];

            let start = Instant::now();
            for _ in 0..ITERATIONS {
                for chunk_start in (0..PIXELS).step_by(8) {
                    if chunk_start + 8 <= PIXELS {
                        let r = f32x8::load(
                            token,
                            (&r_data[chunk_start..chunk_start + 8]).try_into().unwrap(),
                        );
                        let g = f32x8::load(
                            token,
                            (&g_data[chunk_start..chunk_start + 8]).try_into().unwrap(),
                        );
                        let b = f32x8::load(
                            token,
                            (&b_data[chunk_start..chunk_start + 8]).try_into().unwrap(),
                        );

                        let (y, cb, cr) = rgb_to_ycbcr_8px(token, r, g, b);

                        y.store(
                            (&mut y_data[chunk_start..chunk_start + 8])
                                .try_into()
                                .unwrap(),
                        );
                        cb.store(
                            (&mut cb_data[chunk_start..chunk_start + 8])
                                .try_into()
                                .unwrap(),
                        );
                        cr.store(
                            (&mut cr_data[chunk_start..chunk_start + 8])
                                .try_into()
                                .unwrap(),
                        );
                    }
                }
                std::hint::black_box(&y_data);
            }
            let simd_time = start.elapsed();
            let mpix_s = (PIXELS * ITERATIONS) as f64 / simd_time.as_secs_f64() / 1_000_000.0;
            println!(
                "  RGB→YCbCr:       {:>8.2} ms ({:.1} Mpix/s)",
                simd_time.as_secs_f64() * 1000.0,
                mpix_s
            );

            // Blend modes benchmark
            let src_data: Vec<f32> = (0..PIXELS)
                .map(|i| ((i * 17) % 256) as f32 / 255.0)
                .collect();
            let dst_data: Vec<f32> = (0..PIXELS)
                .map(|i| ((i * 31) % 256) as f32 / 255.0)
                .collect();
            let mut out_data = vec![0.0f32; PIXELS];

            let start = Instant::now();
            for _ in 0..ITERATIONS {
                for chunk_start in (0..PIXELS).step_by(8) {
                    if chunk_start + 8 <= PIXELS {
                        let src = f32x8::load(
                            token,
                            (&src_data[chunk_start..chunk_start + 8])
                                .try_into()
                                .unwrap(),
                        );
                        let dst = f32x8::load(
                            token,
                            (&dst_data[chunk_start..chunk_start + 8])
                                .try_into()
                                .unwrap(),
                        );
                        let result = blend_overlay_2px(token, src, dst);
                        result.store(
                            (&mut out_data[chunk_start..chunk_start + 8])
                                .try_into()
                                .unwrap(),
                        );
                    }
                }
                std::hint::black_box(&out_data);
            }
            let simd_time = start.elapsed();
            let mpix_s = (PIXELS * ITERATIONS) as f64 / simd_time.as_secs_f64() / 1_000_000.0;
            println!(
                "  Overlay blend:   {:>8.2} ms ({:.1} Mpix/s)",
                simd_time.as_secs_f64() * 1000.0,
                mpix_s
            );
        }

        println!();
    }

    pub fn main() {
        println!("\n╔═══════════════════════════════════════════════════════════════╗");
        println!("║     SIMD Kernel Collection - Image Processing Hotspots        ║");
        println!("╚═══════════════════════════════════════════════════════════════╝\n");

        println!("Eight kernels from image-webp, jpegli-rs, and zenimage:\n");
        println!("  1. 4x4 DCT (VP8)        - Integer transform");
        println!("  2. 8x8 DCT butterfly    - JPEG encoding");
        println!("  3. Chroma downsample    - 2x2 box filter");
        println!("  4. RGB→YCbCr            - Color matrix FMA");
        println!("  5. Horizontal convolve  - 1D filter");
        println!("  6. sRGB↔Linear          - Gamma correction");
        println!("  7. Blend modes          - Multiply/Screen/Overlay");
        println!("  8. Horizontal reduce    - Strided weighted sum\n");

        test_correctness();
        benchmark();

        println!("=== Key Patterns ===\n");
        println!("  Token-gated dispatch:");
        println!("    if let Some(token) = X64V3Token::summon() {{ ... }}");
        println!();
        println!("  FMA chains for matrix ops:");
        println!("    y = r.mul_add(ky_r, g.mul_add(ky_g, b * ky_b));");
        println!();
        println!("  Gamma approximation:");
        println!("    x^2.4 ≈ x^2 * sqrt(sqrt(x)) * sqrt(sqrt(sqrt(x)))");
        println!();
    }
}

#[cfg(target_arch = "x86_64")]
fn main() {
    x86_impl::main()
}