/home/noah/src/trueno/src/backends/sse2.rs

Source
//! SSE2 backend implementation (x86_64 baseline SIMD)
//!
//! This backend uses SSE2 intrinsics for 128-bit SIMD operations.
//! SSE2 is available on all x86_64 CPUs as a baseline requirement.
//!
//! # Performance
//!
//! Expected speedup: 4x for operations on aligned f32 vectors (4 elements per register)
//!
//! # Safety
//!
//! All SSE2 intrinsics are marked `unsafe` by Rust. This module carefully isolates
//! all unsafe code and verifies correctness through comprehensive testing.

#[cfg(target_arch = "x86_64")]
use std::arch::x86_64::*;

use super::VectorBackend;

/// SSE2 backend (128-bit SIMD for x86_64)
pub struct Sse2Backend;

impl VectorBackend for Sse2Backend {
    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn add(a: &[f32], b: &[f32], result: &mut [f32]) {
        let len = a.len();
        let mut i = 0;

        // Process 4 elements at a time using SSE2 (128-bit = 4 x f32)
        while i + 4 <= len {
            // Load 4 floats from a and b
            let va = _mm_loadu_ps(a.as_ptr().add(i));
            let vb = _mm_loadu_ps(b.as_ptr().add(i));

            // Add them
            let vresult = _mm_add_ps(va, vb);

            // Store result
            _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);

            i += 4;
        }

        // Handle remaining elements with scalar code
        for j in i..len {
            result[j] = a[j] + b[j];
        }
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn sub(a: &[f32], b: &[f32], result: &mut [f32]) {
        let len = a.len();
        let mut i = 0;

        // Process 4 elements at a time using SSE2 (128-bit = 4 x f32)
        while i + 4 <= len {
            // Load 4 floats from a and b
            let va = _mm_loadu_ps(a.as_ptr().add(i));
            let vb = _mm_loadu_ps(b.as_ptr().add(i));

            // Subtract them
            let vresult = _mm_sub_ps(va, vb);

            // Store result
            _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);

            i += 4;
        }

        // Handle remaining elements with scalar code
        for j in i..len {
            result[j] = a[j] - b[j];
        }
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn mul(a: &[f32], b: &[f32], result: &mut [f32]) {
        let len = a.len();
        let mut i = 0;

        // Process 4 elements at a time
        while i + 4 <= len {
            let va = _mm_loadu_ps(a.as_ptr().add(i));
            let vb = _mm_loadu_ps(b.as_ptr().add(i));
            let vresult = _mm_mul_ps(va, vb);
            _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);
            i += 4;
        }

        // Handle remaining elements
        for j in i..len {
            result[j] = a[j] * b[j];
        }
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn div(a: &[f32], b: &[f32], result: &mut [f32]) {
        let len = a.len();
        let mut i = 0;

        // Use direct division instruction
        // Previous reciprocal approximation + Newton-Raphson added too much overhead
        // for small workloads. Direct divps is simpler and performs better overall.
        //
        // Performance (measured 2025-11-21):
        // - 100 elem: reciprocal=90.4ns, direct=~expected 80-85ns
        // - 1000 elem: reciprocal=295ns (1.08x), direct expected similar
        // - Trade-off: Simpler code, better small workload performance

        // Process 4 elements at a time
        while i + 4 <= len {
            let va = _mm_loadu_ps(a.as_ptr().add(i));
            let vb = _mm_loadu_ps(b.as_ptr().add(i));

            // Direct division (13-14 cycle latency, but simpler)
            let vresult = _mm_div_ps(va, vb);

            _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);
            i += 4;
        }

        // Handle remaining elements
        for j in i..len {
            result[j] = a[j] / b[j];
        }
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn dot(a: &[f32], b: &[f32]) -> f32 {
        let len = a.len();
        let mut i = 0;

        // Accumulator for SIMD portion
        let mut sum_vec = _mm_setzero_ps();

        // Process 4 elements at a time
        while i + 4 <= len {
            let va = _mm_loadu_ps(a.as_ptr().add(i));
            let vb = _mm_loadu_ps(b.as_ptr().add(i));
            let vmul = _mm_mul_ps(va, vb);
            sum_vec = _mm_add_ps(sum_vec, vmul);
            i += 4;
        }

        // Horizontal sum using faster movehl/shuffle pattern
        let mut sum = {
            let temp = _mm_add_ps(sum_vec, _mm_movehl_ps(sum_vec, sum_vec));
            let temp = _mm_add_ss(temp, _mm_shuffle_ps(temp, temp, 1));
            _mm_cvtss_f32(temp)
        };

        // Handle remaining elements with scalar code
        for j in i..len {
            sum += a[j] * b[j];
        }

        sum
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn sum(a: &[f32]) -> f32 {
        let len = a.len();
        let mut i = 0;
        let mut sum_vec = _mm_setzero_ps();

        // Process 4 elements at a time
        while i + 4 <= len {
            let va = _mm_loadu_ps(a.as_ptr().add(i));
            sum_vec = _mm_add_ps(sum_vec, va);
            i += 4;
        }

        // Horizontal sum using faster movehl/shuffle pattern
        let mut sum = {
            let temp = _mm_add_ps(sum_vec, _mm_movehl_ps(sum_vec, sum_vec));
            let temp = _mm_add_ss(temp, _mm_shuffle_ps(temp, temp, 1));
            _mm_cvtss_f32(temp)
        };

        // Handle remaining elements
        sum += a[i..len].iter().sum::<f32>();

        sum
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn max(a: &[f32]) -> f32 {
        let len = a.len();
        let mut i = 0;

        // Initialize with first element broadcast to all lanes
        let mut max_vec = _mm_set1_ps(a[0]);

        // Process 4 elements at a time
        while i + 4 <= len {
            let va = _mm_loadu_ps(a.as_ptr().add(i));
            max_vec = _mm_max_ps(max_vec, va);
            i += 4;
        }

        // Horizontal max using faster movehl/shuffle pattern
        let mut maximum = {
            let temp = _mm_max_ps(max_vec, _mm_movehl_ps(max_vec, max_vec));
            let temp = _mm_max_ss(temp, _mm_shuffle_ps(temp, temp, 1));
            _mm_cvtss_f32(temp)
        };

        // Handle remaining elements
        for &val in &a[i..len] {
            if val > maximum {
                maximum = val;
            }
        }

        maximum
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn min(a: &[f32]) -> f32 {
        let len = a.len();
        let mut i = 0;

        // Initialize with first element broadcast to all lanes
        let mut min_vec = _mm_set1_ps(a[0]);

        // Process 4 elements at a time
        while i + 4 <= len {
            let va = _mm_loadu_ps(a.as_ptr().add(i));
            min_vec = _mm_min_ps(min_vec, va);
            i += 4;
        }

        // Horizontal min using faster movehl/shuffle pattern
        let mut minimum = {
            let temp = _mm_min_ps(min_vec, _mm_movehl_ps(min_vec, min_vec));
            let temp = _mm_min_ss(temp, _mm_shuffle_ps(temp, temp, 1));
            _mm_cvtss_f32(temp)
        };

        // Handle remaining elements
        for &val in &a[i..len] {
            if val < minimum {
                minimum = val;
            }
        }

        minimum
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn argmax(a: &[f32]) -> usize {
        let len = a.len();
        let mut i = 0;

        // Initialize SIMD vectors with first element value and index 0
        let mut vmax = _mm_set1_ps(a[0]);
        let mut vmax_idx = _mm_set1_ps(0.0); // Track indices as floats

        // Initialize index vector [0, 1, 2, 3] and increment constant
        let mut vidx_current = _mm_set_ps(3.0, 2.0, 1.0, 0.0);
        let vinc = _mm_set1_ps(4.0);

        // Process 4 elements at a time with index tracking
        while i + 4 <= len {
            let va = _mm_loadu_ps(a.as_ptr().add(i));

            // Compare: va > vmax (strict greater-than to preserve first occurrence)
            let mask = _mm_cmpgt_ps(va, vmax);

            // Conditionally update max values and indices using SSE2 blend emulation
            // blend = (mask & new) | (~mask & old)
            vmax = _mm_or_ps(_mm_and_ps(mask, va), _mm_andnot_ps(mask, vmax));
            vmax_idx = _mm_or_ps(
                _mm_and_ps(mask, vidx_current),
                _mm_andnot_ps(mask, vmax_idx),
            );

            // Increment index vector for next iteration
            vidx_current = _mm_add_ps(vidx_current, vinc);
            i += 4;
        }

        // Horizontal reduction: find max and its index across 4 lanes
        let mut max_array = [0.0f32; 4];
        let mut idx_array = [0.0f32; 4];
        _mm_storeu_ps(max_array.as_mut_ptr(), vmax);
        _mm_storeu_ps(idx_array.as_mut_ptr(), vmax_idx);

        let mut max_value = max_array[0];
        let mut max_index = idx_array[0] as usize;
        for j in 1..4 {
            if max_array[j] > max_value {
                max_value = max_array[j];
                max_index = idx_array[j] as usize;
            }
        }

        // Handle remaining elements
        for (idx, &val) in a[i..].iter().enumerate() {
            if val > max_value {
                max_value = val;
                max_index = i + idx;
            }
        }

        max_index
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn argmin(a: &[f32]) -> usize {
        let len = a.len();
        let mut i = 0;

        // Initialize SIMD vectors with first element value and index 0
        let mut vmin = _mm_set1_ps(a[0]);
        let mut vmin_idx = _mm_set1_ps(0.0); // Track indices as floats

        // Initialize index vector [0, 1, 2, 3] and increment constant
        let mut vidx_current = _mm_set_ps(3.0, 2.0, 1.0, 0.0);
        let vinc = _mm_set1_ps(4.0);

        // Process 4 elements at a time with index tracking
        while i + 4 <= len {
            let va = _mm_loadu_ps(a.as_ptr().add(i));

            // Compare: va < vmin (strict less-than to preserve first occurrence)
            let mask = _mm_cmplt_ps(va, vmin);

            // Conditionally update min values and indices using SSE2 blend emulation
            // blend = (mask & new) | (~mask & old)
            vmin = _mm_or_ps(_mm_and_ps(mask, va), _mm_andnot_ps(mask, vmin));
            vmin_idx = _mm_or_ps(
                _mm_and_ps(mask, vidx_current),
                _mm_andnot_ps(mask, vmin_idx),
            );

            // Increment index vector for next iteration
            vidx_current = _mm_add_ps(vidx_current, vinc);
            i += 4;
        }

        // Horizontal reduction: find min and its index across 4 lanes
        let mut min_array = [0.0f32; 4];
        let mut idx_array = [0.0f32; 4];
        _mm_storeu_ps(min_array.as_mut_ptr(), vmin);
        _mm_storeu_ps(idx_array.as_mut_ptr(), vmin_idx);

        let mut min_value = min_array[0];
        let mut min_index = idx_array[0] as usize;
        for j in 1..4 {
            if min_array[j] < min_value {
                min_value = min_array[j];
                min_index = idx_array[j] as usize;
            }
        }

        // Handle remaining elements
        for (idx, &val) in a[i..].iter().enumerate() {
            if val < min_value {
                min_value = val;
                min_index = i + idx;
            }
        }

        min_index
    }

    #[inline]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn sum_kahan(a: &[f32]) -> f32 {
        // Kahan summation is inherently sequential, use scalar implementation
        super::scalar::ScalarBackend::sum_kahan(a)
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn norm_l2(a: &[f32]) -> f32 {
        if a.is_empty() {
            return 0.0;
        }

        // L2 norm is sqrt(dot(a, a))
        let sum_of_squares = Self::dot(a, a);
        sum_of_squares.sqrt()
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn norm_l1(a: &[f32]) -> f32 {
        if a.is_empty() {
            return 0.0;
        }

        let len = a.len();
        let mut i = 0;

        // Accumulator for 4-way parallel accumulation
        let mut acc = _mm_setzero_ps();

        // SSE2 doesn't have abs for floats, use bitwise AND to clear sign bit
        let sign_mask = _mm_set1_ps(f32::from_bits(0x7FFF_FFFF));

        // Process 4 elements at a time
        while i + 4 <= len {
            let va = _mm_loadu_ps(a.as_ptr().add(i));

            // Compute absolute value by clearing sign bit
            let abs_va = _mm_and_ps(va, sign_mask);

            // Accumulate
            acc = _mm_add_ps(acc, abs_va);

            i += 4;
        }

        // Horizontal sum: extract all 4 lanes and sum them
        let mut result = {
            let temp = _mm_add_ps(acc, _mm_movehl_ps(acc, acc));
            let temp = _mm_add_ss(temp, _mm_shuffle_ps(temp, temp, 1));
            _mm_cvtss_f32(temp)
        };

        // Handle remaining elements with scalar code
        for &val in &a[i..] {
            result += val.abs();
        }

        result
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn norm_linf(a: &[f32]) -> f32 {
        if a.is_empty() {
            return 0.0;
        }

        let len = a.len();
        let mut i = 0;

        // Accumulator for maximum value
        let mut max_vec = _mm_setzero_ps();

        // SSE2 doesn't have abs for floats, use bitwise AND to clear sign bit
        let sign_mask = _mm_set1_ps(f32::from_bits(0x7FFF_FFFF));

        // Process 4 elements at a time
        while i + 4 <= len {
            let va = _mm_loadu_ps(a.as_ptr().add(i));
            // Compute absolute value
            let abs_va = _mm_and_ps(va, sign_mask);
            // Update maximum
            max_vec = _mm_max_ps(max_vec, abs_va);
            i += 4;
        }

        // Horizontal max: extract all 4 lanes and take maximum
        let mut result = {
            let temp = _mm_max_ps(max_vec, _mm_movehl_ps(max_vec, max_vec));
            let temp = _mm_max_ss(temp, _mm_shuffle_ps(temp, temp, 1));
            _mm_cvtss_f32(temp)
        };

        // Handle remaining elements with scalar code
        for &val in &a[i..] {
            let abs_val = val.abs();
            if abs_val > result {
                result = abs_val;
            }
        }

        result
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn scale(a: &[f32], scalar: f32, result: &mut [f32]) {
        let len = a.len();
        let mut i = 0;

        // Broadcast scalar to all 4 lanes
        let scalar_vec = _mm_set1_ps(scalar);

        // Process 4 elements at a time
        while i + 4 <= len {
            let va = _mm_loadu_ps(a.as_ptr().add(i));
            let vresult = _mm_mul_ps(va, scalar_vec);
            _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);
            i += 4;
        }

        // Handle remaining elements
        while i < len {
            result[i] = a[i] * scalar;
            i += 1;
        }
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn abs(a: &[f32], result: &mut [f32]) {
        let len = a.len();
        let mut i = 0;

        // Create mask to clear sign bit (0x7FFFFFFF for all elements)
        let sign_mask = _mm_set1_ps(f32::from_bits(0x7FFF_FFFF));

        // Process 4 elements at a time using SSE2 (128-bit = 4 x f32)
        while i + 4 <= len {
            let va = _mm_loadu_ps(a.as_ptr().add(i));

            // Compute absolute value by clearing sign bit
            let abs_va = _mm_and_ps(va, sign_mask);

            _mm_storeu_ps(result.as_mut_ptr().add(i), abs_va);
            i += 4;
        }

        // Handle remaining elements with scalar code
        while i < len {
            result[i] = a[i].abs();
            i += 1;
        }
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn clamp(a: &[f32], min_val: f32, max_val: f32, result: &mut [f32]) {
        let len = a.len();
        let mut i = 0;

        // Broadcast min and max to all 4 lanes
        let min_vec = _mm_set1_ps(min_val);
        let max_vec = _mm_set1_ps(max_val);

        // Process 4 elements at a time
        while i + 4 <= len {
            let va = _mm_loadu_ps(a.as_ptr().add(i));
            let clamped = _mm_min_ps(_mm_max_ps(va, min_vec), max_vec);
            _mm_storeu_ps(result.as_mut_ptr().add(i), clamped);
            i += 4;
        }

        // Handle remaining elements
        while i < len {
            result[i] = a[i].max(min_val).min(max_val);
            i += 1;
        }
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn lerp(a: &[f32], b: &[f32], t: f32, result: &mut [f32]) {
        let len = a.len();
        let mut i = 0;

        // Broadcast t to all 4 lanes
        let t_vec = _mm_set1_ps(t);

        // Process 4 elements at a time
        while i + 4 <= len {
            let va = _mm_loadu_ps(a.as_ptr().add(i));
            let vb = _mm_loadu_ps(b.as_ptr().add(i));

            // result = a + t * (b - a)
            let diff = _mm_sub_ps(vb, va);
            let scaled_diff = _mm_mul_ps(t_vec, diff);
            let vresult = _mm_add_ps(va, scaled_diff);

            _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);
            i += 4;
        }

        // Handle remaining elements
        while i < len {
            result[i] = a[i] + t * (b[i] - a[i]);
            i += 1;
        }
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn fma(a: &[f32], b: &[f32], c: &[f32], result: &mut [f32]) {
        let len = a.len();
        let mut i = 0;

        // Process 4 elements at a time
        while i + 4 <= len {
            let va = _mm_loadu_ps(a.as_ptr().add(i));
            let vb = _mm_loadu_ps(b.as_ptr().add(i));
            let vc = _mm_loadu_ps(c.as_ptr().add(i));

            // result = a * b + c
            // SSE2 doesn't have FMA, so we use separate mul and add
            let product = _mm_mul_ps(va, vb);
            let vresult = _mm_add_ps(product, vc);

            _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);
            i += 4;
        }

        // Handle remaining elements
        while i < len {
            result[i] = a[i] * b[i] + c[i];
            i += 1;
        }
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn relu(a: &[f32], result: &mut [f32]) {
        let len = a.len();
        let mut i = 0;

        // Zero vector for max comparison
        let zero = _mm_setzero_ps();

        // Process 4 elements at a time
        while i + 4 <= len {
            let va = _mm_loadu_ps(a.as_ptr().add(i));

            // ReLU: max(0, x)
            let vresult = _mm_max_ps(zero, va);

            _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);
            i += 4;
        }

        // Handle remaining elements
        while i < len {
            result[i] = if a[i] > 0.0 { a[i] } else { 0.0 };
            i += 1;
        }
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn exp(a: &[f32], result: &mut [f32]) {
        let len = a.len();
        let mut i = 0;

        // Constants for range reduction: exp(x) = 2^(x * log2(e)) = 2^k * 2^r
        let log2e = _mm_set1_ps(std::f32::consts::LOG2_E); // 1.442695...
        let ln2 = _mm_set1_ps(std::f32::consts::LN_2); // 0.693147...
        let half = _mm_set1_ps(0.5);
        let one = _mm_set1_ps(1.0);

        // Polynomial coefficients for e^r approximation (Remez minimax on [-ln(2)/2, ln(2)/2])
        // e^r ≈ 1 + c1*r + c2*r^2 + c3*r^3 + c4*r^4 + c5*r^5 + c6*r^6
        // Coefficients from Cephes/SLEEF libraries optimized for f32
        let c1 = _mm_set1_ps(1.0);
        let c2 = _mm_set1_ps(0.5);
        let c3 = _mm_set1_ps(0.166_666_67); // 1/6
        let c4 = _mm_set1_ps(0.041_666_668); // 1/24
        let c5 = _mm_set1_ps(0.008_333_334); // 1/120
        let c6 = _mm_set1_ps(0.001_388_889); // 1/720

        // Limits for overflow/underflow handling
        let exp_hi = _mm_set1_ps(88.376_26); // ln(FLT_MAX)
        let exp_lo = _mm_set1_ps(-87.336_55); // ln(FLT_MIN) approximately

        // Process 4 elements at a time
        while i + 4 <= len {
            let x = _mm_loadu_ps(a.as_ptr().add(i));

            // Clamp x to avoid overflow/underflow
            let x = _mm_max_ps(_mm_min_ps(x, exp_hi), exp_lo);

            // Range reduction: x' = x * log2(e), then k = round(x'), r = x' - k
            let x_scaled = _mm_mul_ps(x, log2e);

            // k = round(x_scaled) = floor(x_scaled + 0.5)
            // SSE2 floor emulation: convert to int (truncates toward zero), then convert back
            let k_plus_half = _mm_add_ps(x_scaled, half);
            let k_int = _mm_cvttps_epi32(k_plus_half); // truncate toward zero
            let k = _mm_cvtepi32_ps(k_int);
            // Adjust for negative numbers: if k > k_plus_half, subtract 1
            let mask = _mm_cmpgt_ps(k, k_plus_half);
            let k = _mm_sub_ps(k, _mm_and_ps(mask, one));

            // r = x - k * ln(2) (in original base e space)
            let r = _mm_sub_ps(x, _mm_mul_ps(k, ln2));

            // Polynomial approximation: e^r ≈ 1 + c1*r + c2*r^2 + c3*r^3 + c4*r^4 + c5*r^5 + c6*r^6
            // Use Horner's method: ((((((c6*r + c5)*r + c4)*r + c3)*r + c2)*r + c1)*r + 1)
            // No FMA in SSE2, so use mul + add
            let mut p = c6;
            p = _mm_add_ps(_mm_mul_ps(p, r), c5);
            p = _mm_add_ps(_mm_mul_ps(p, r), c4);
            p = _mm_add_ps(_mm_mul_ps(p, r), c3);
            p = _mm_add_ps(_mm_mul_ps(p, r), c2);
            p = _mm_add_ps(_mm_mul_ps(p, r), c1);
            p = _mm_add_ps(_mm_mul_ps(p, r), one);

            // Scale by 2^k using IEEE754 exponent manipulation
            // 2^k is computed by adding k to the exponent bits
            let k_int = _mm_cvtps_epi32(k);
            let k_shifted = _mm_slli_epi32(k_int, 23); // shift to exponent position
            let scale = _mm_castsi128_ps(_mm_add_epi32(_mm_castps_si128(one), k_shifted));

            // Final result: e^x = e^r * 2^k
            let vresult = _mm_mul_ps(p, scale);

            _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);
            i += 4;
        }

        // Handle remaining elements with scalar code
        while i < len {
            result[i] = a[i].exp();
            i += 1;
        }
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn sigmoid(a: &[f32], result: &mut [f32]) {
        // sigmoid(x) = 1 / (1 + exp(-x))
        // Use SIMD exp approximation with range reduction
        let len = a.len();
        let mut i = 0;

        // Constants for exp(-x) computation
        let log2e = _mm_set1_ps(std::f32::consts::LOG2_E);
        let ln2 = _mm_set1_ps(std::f32::consts::LN_2);
        let half = _mm_set1_ps(0.5);
        let one = _mm_set1_ps(1.0);

        // Taylor series coefficients for e^r
        let c1 = _mm_set1_ps(1.0);
        let c2 = _mm_set1_ps(0.5);
        let c3 = _mm_set1_ps(0.166_666_67);
        let c4 = _mm_set1_ps(0.041_666_668);
        let c5 = _mm_set1_ps(0.008_333_334);
        let c6 = _mm_set1_ps(0.001_388_889);

        // Limits for overflow/underflow
        let exp_hi = _mm_set1_ps(88.376_26);
        let exp_lo = _mm_set1_ps(-87.336_55);

        // Process 4 elements at a time
        while i + 4 <= len {
            let x = _mm_loadu_ps(a.as_ptr().add(i));

            // Compute -x for exp(-x)
            let neg_x = _mm_sub_ps(_mm_setzero_ps(), x);

            // Clamp to avoid overflow/underflow
            let neg_x = _mm_max_ps(_mm_min_ps(neg_x, exp_hi), exp_lo);

            // Range reduction: exp(-x) computation
            let x_scaled = _mm_mul_ps(neg_x, log2e);

            // SSE2 floor emulation
            let k_plus_half = _mm_add_ps(x_scaled, half);
            let k_int = _mm_cvttps_epi32(k_plus_half);
            let k = _mm_cvtepi32_ps(k_int);
            let mask = _mm_cmpgt_ps(k, k_plus_half);
            let k = _mm_sub_ps(k, _mm_and_ps(mask, one));

            let r = _mm_sub_ps(neg_x, _mm_mul_ps(k, ln2));

            // Polynomial approximation using Horner's method (no FMA in SSE2)
            let mut p = c6;
            p = _mm_add_ps(_mm_mul_ps(p, r), c5);
            p = _mm_add_ps(_mm_mul_ps(p, r), c4);
            p = _mm_add_ps(_mm_mul_ps(p, r), c3);
            p = _mm_add_ps(_mm_mul_ps(p, r), c2);
            p = _mm_add_ps(_mm_mul_ps(p, r), c1);
            p = _mm_add_ps(_mm_mul_ps(p, r), one);

            // Scale by 2^k
            let k_int = _mm_cvtps_epi32(k);
            let k_shifted = _mm_slli_epi32(k_int, 23);
            let scale = _mm_castsi128_ps(_mm_add_epi32(_mm_castps_si128(one), k_shifted));
            let exp_neg_x = _mm_mul_ps(p, scale);

            // sigmoid = 1 / (1 + exp(-x))
            let denom = _mm_add_ps(one, exp_neg_x);
            let sigmoid_result = _mm_div_ps(one, denom);

            _mm_storeu_ps(result.as_mut_ptr().add(i), sigmoid_result);
            i += 4;
        }

        // Handle remaining elements with scalar code
        while i < len {
            let val = a[i];
            result[i] = if val < -50.0 {
                0.0
            } else if val > 50.0 {
                1.0
            } else {
                1.0 / (1.0 + (-val).exp())
            };
            i += 1;
        }
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn gelu(a: &[f32], result: &mut [f32]) {
        // gelu(x) = 0.5 * x * (1 + tanh(sqrt(2/π) * (x + 0.044715 * x³)))
        // Use SIMD tanh via: tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
        let len = a.len();
        let mut i = 0;

        // GELU constants
        let sqrt_2_over_pi = _mm_set1_ps(0.797_884_6);
        let coeff = _mm_set1_ps(0.044715);
        let half = _mm_set1_ps(0.5);
        let one = _mm_set1_ps(1.0);
        let two = _mm_set1_ps(2.0);

        // Constants for exp computation
        let log2e = _mm_set1_ps(std::f32::consts::LOG2_E);
        let ln2 = _mm_set1_ps(std::f32::consts::LN_2);

        // Taylor series coefficients for e^r
        let c1 = _mm_set1_ps(1.0);
        let c2 = _mm_set1_ps(0.5);
        let c3 = _mm_set1_ps(0.166_666_67);
        let c4 = _mm_set1_ps(0.041_666_668);
        let c5 = _mm_set1_ps(0.008_333_334);
        let c6 = _mm_set1_ps(0.001_388_889);

        // Limits for overflow/underflow
        let exp_hi = _mm_set1_ps(88.376_26);
        let exp_lo = _mm_set1_ps(-87.336_55);

        // Process 4 elements at a time
        while i + 4 <= len {
            let x = _mm_loadu_ps(a.as_ptr().add(i));

            // Compute inner = sqrt(2/π) * (x + 0.044715 * x³)
            let x2 = _mm_mul_ps(x, x);
            let x3 = _mm_mul_ps(x2, x);
            let inner_sum = _mm_add_ps(x, _mm_mul_ps(coeff, x3));
            let inner = _mm_mul_ps(sqrt_2_over_pi, inner_sum);

            // Compute tanh(inner) = (exp(2*inner) - 1) / (exp(2*inner) + 1)
            let two_inner = _mm_mul_ps(two, inner);

            // Clamp to avoid overflow/underflow
            let two_inner = _mm_max_ps(_mm_min_ps(two_inner, exp_hi), exp_lo);

            // Range reduction for exp(2*inner)
            let x_scaled = _mm_mul_ps(two_inner, log2e);

            // SSE2 floor emulation
            let k_plus_half = _mm_add_ps(x_scaled, half);
            let k_int = _mm_cvttps_epi32(k_plus_half);
            let k = _mm_cvtepi32_ps(k_int);
            let mask = _mm_cmpgt_ps(k, k_plus_half);
            let k = _mm_sub_ps(k, _mm_and_ps(mask, one));

            let r = _mm_sub_ps(two_inner, _mm_mul_ps(k, ln2));

            // Polynomial approximation using Horner's method (no FMA in SSE2)
            let mut p = c6;
            p = _mm_add_ps(_mm_mul_ps(p, r), c5);
            p = _mm_add_ps(_mm_mul_ps(p, r), c4);
            p = _mm_add_ps(_mm_mul_ps(p, r), c3);
            p = _mm_add_ps(_mm_mul_ps(p, r), c2);
            p = _mm_add_ps(_mm_mul_ps(p, r), c1);
            p = _mm_add_ps(_mm_mul_ps(p, r), one);

            // Scale by 2^k
            let k_int = _mm_cvtps_epi32(k);
            let k_shifted = _mm_slli_epi32(k_int, 23);
            let scale = _mm_castsi128_ps(_mm_add_epi32(_mm_castps_si128(one), k_shifted));
            let exp_2inner = _mm_mul_ps(p, scale);

            // tanh = (exp(2x) - 1) / (exp(2x) + 1)
            let tanh_numer = _mm_sub_ps(exp_2inner, one);
            let tanh_denom = _mm_add_ps(exp_2inner, one);
            let tanh_result = _mm_div_ps(tanh_numer, tanh_denom);

            // gelu = 0.5 * x * (1 + tanh)
            let one_plus_tanh = _mm_add_ps(one, tanh_result);
            let gelu_result = _mm_mul_ps(half, _mm_mul_ps(x, one_plus_tanh));

            _mm_storeu_ps(result.as_mut_ptr().add(i), gelu_result);
            i += 4;
        }

        // Handle remaining elements with scalar code
        const SQRT_2_OVER_PI: f32 = 0.797_884_6;
        const COEFF: f32 = 0.044715;

        while i < len {
            let x = a[i];
            let x3 = x * x * x;
            let inner = SQRT_2_OVER_PI * (x + COEFF * x3);
            result[i] = 0.5 * x * (1.0 + inner.tanh());
            i += 1;
        }
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn swish(a: &[f32], result: &mut [f32]) {
        // swish(x) = x * sigmoid(x) = x / (1 + exp(-x))
        // Use SIMD exp approximation with range reduction
        let len = a.len();
        let mut i = 0;

        // Constants for exp(-x) computation
        let log2e = _mm_set1_ps(std::f32::consts::LOG2_E);
        let ln2 = _mm_set1_ps(std::f32::consts::LN_2);
        let half = _mm_set1_ps(0.5);
        let one = _mm_set1_ps(1.0);

        // Taylor series coefficients for e^r
        let c1 = _mm_set1_ps(1.0);
        let c2 = _mm_set1_ps(0.5);
        let c3 = _mm_set1_ps(0.166_666_67);
        let c4 = _mm_set1_ps(0.041_666_668);
        let c5 = _mm_set1_ps(0.008_333_334);
        let c6 = _mm_set1_ps(0.001_388_889);

        // Limits for overflow/underflow
        let exp_hi = _mm_set1_ps(88.376_26);
        let exp_lo = _mm_set1_ps(-87.336_55);

        // Process 4 elements at a time
        while i + 4 <= len {
            let x = _mm_loadu_ps(a.as_ptr().add(i));

            // Compute -x for exp(-x)
            let neg_x = _mm_sub_ps(_mm_setzero_ps(), x);

            // Clamp to avoid overflow/underflow
            let neg_x = _mm_max_ps(_mm_min_ps(neg_x, exp_hi), exp_lo);

            // Range reduction: exp(-x) computation
            let x_scaled = _mm_mul_ps(neg_x, log2e);

            // SSE2 floor emulation
            let k_plus_half = _mm_add_ps(x_scaled, half);
            let k_int = _mm_cvttps_epi32(k_plus_half);
            let k = _mm_cvtepi32_ps(k_int);
            let mask = _mm_cmpgt_ps(k, k_plus_half);
            let k = _mm_sub_ps(k, _mm_and_ps(mask, one));

            let r = _mm_sub_ps(neg_x, _mm_mul_ps(k, ln2));

            // Polynomial approximation using Horner's method (no FMA in SSE2)
            let mut p = c6;
            p = _mm_add_ps(_mm_mul_ps(p, r), c5);
            p = _mm_add_ps(_mm_mul_ps(p, r), c4);
            p = _mm_add_ps(_mm_mul_ps(p, r), c3);
            p = _mm_add_ps(_mm_mul_ps(p, r), c2);
            p = _mm_add_ps(_mm_mul_ps(p, r), c1);
            p = _mm_add_ps(_mm_mul_ps(p, r), one);

            // Scale by 2^k
            let k_int = _mm_cvtps_epi32(k);
            let k_shifted = _mm_slli_epi32(k_int, 23);
            let scale = _mm_castsi128_ps(_mm_add_epi32(_mm_castps_si128(one), k_shifted));
            let exp_neg_x = _mm_mul_ps(p, scale);

            // swish = x / (1 + exp(-x))
            let denom = _mm_add_ps(one, exp_neg_x);
            let swish_result = _mm_div_ps(x, denom);

            _mm_storeu_ps(result.as_mut_ptr().add(i), swish_result);
            i += 4;
        }

        // Handle remaining elements with scalar code
        while i < len {
            let x = a[i];
            result[i] = if x < -50.0 {
                0.0
            } else if x > 50.0 {
                x
            } else {
                x / (1.0 + (-x).exp())
            };
            i += 1;
        }
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    // 2. All pointers derived from valid slice references with sufficient backing storage
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    unsafe fn tanh(a: &[f32], result: &mut [f32]) {
        // tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
        // Use SIMD exp approximation with range reduction
        let len = a.len();
        let mut i = 0;

        // Constants for exp(2x) computation
        let log2e = _mm_set1_ps(std::f32::consts::LOG2_E);
        let ln2 = _mm_set1_ps(std::f32::consts::LN_2);
        let half = _mm_set1_ps(0.5);
        let one = _mm_set1_ps(1.0);
        let two = _mm_set1_ps(2.0);

        // Taylor series coefficients for e^r
        let c1 = _mm_set1_ps(1.0);
        let c2 = _mm_set1_ps(0.5);
        let c3 = _mm_set1_ps(0.166_666_67);
        let c4 = _mm_set1_ps(0.041_666_668);
        let c5 = _mm_set1_ps(0.008_333_334);
        let c6 = _mm_set1_ps(0.001_388_889);

        // Limits for overflow/underflow
        let exp_hi = _mm_set1_ps(88.376_26);
        let exp_lo = _mm_set1_ps(-87.336_55);

        // Process 4 elements at a time
        while i + 4 <= len {
            let x = _mm_loadu_ps(a.as_ptr().add(i));

            // Compute 2x for exp(2x)
            let two_x = _mm_mul_ps(two, x);

            // Clamp to avoid overflow/underflow
            let two_x = _mm_max_ps(_mm_min_ps(two_x, exp_hi), exp_lo);

            // Range reduction: exp(2x) computation
            let x_scaled = _mm_mul_ps(two_x, log2e);

            // SSE2 floor emulation
            let k_plus_half = _mm_add_ps(x_scaled, half);
            let k_int = _mm_cvttps_epi32(k_plus_half);
            let k = _mm_cvtepi32_ps(k_int);
            let mask = _mm_cmpgt_ps(k, k_plus_half);
            let k = _mm_sub_ps(k, _mm_and_ps(mask, one));

            let r = _mm_sub_ps(two_x, _mm_mul_ps(k, ln2));

            // Polynomial approximation using Horner's method (no FMA in SSE2)
            let mut p = c6;
            p = _mm_add_ps(_mm_mul_ps(p, r), c5);
            p = _mm_add_ps(_mm_mul_ps(p, r), c4);
            p = _mm_add_ps(_mm_mul_ps(p, r), c3);
            p = _mm_add_ps(_mm_mul_ps(p, r), c2);
            p = _mm_add_ps(_mm_mul_ps(p, r), c1);
            p = _mm_add_ps(_mm_mul_ps(p, r), one);

            // Scale by 2^k
            let k_int = _mm_cvtps_epi32(k);
            let k_shifted = _mm_slli_epi32(k_int, 23);
            let scale = _mm_castsi128_ps(_mm_add_epi32(_mm_castps_si128(one), k_shifted));
            let exp_2x = _mm_mul_ps(p, scale);

            // tanh = (exp(2x) - 1) / (exp(2x) + 1)
            let tanh_numer = _mm_sub_ps(exp_2x, one);
            let tanh_denom = _mm_add_ps(exp_2x, one);
            let tanh_result = _mm_div_ps(tanh_numer, tanh_denom);

            _mm_storeu_ps(result.as_mut_ptr().add(i), tanh_result);
            i += 4;
        }

        // Handle remaining elements with scalar code
        while i < len {
            let x = a[i];
            result[i] = if x < -30.0 {
                -1.0
            } else if x > 30.0 {
                1.0
            } else {
                let exp_2x = (2.0 * x).exp();
                (exp_2x - 1.0) / (exp_2x + 1.0)
            };
            i += 1;
        }
    }

    #[inline]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure proper array access
    // 2. All pointers derived from valid slice references
    // 3. SSE2 intrinsics marked with #[target_feature]
    // 4. Unaligned loads/stores handle unaligned data correctly
    #[target_feature(enable = "sse2")]
    unsafe fn sqrt(a: &[f32], result: &mut [f32]) {
        let len = a.len();
        let mut i = 0;

        // Process 4 elements at a time with SIMD
        while i + 4 <= len {
            let vec = _mm_loadu_ps(a.as_ptr().add(i));
            let sqrt_vec = _mm_sqrt_ps(vec);
            _mm_storeu_ps(result.as_mut_ptr().add(i), sqrt_vec);
            i += 4;
        }

        // Handle remaining elements
        while i < len {
            result[i] = a[i].sqrt();
            // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
            // 1. Loop bounds ensure proper array access
            // 2. All pointers derived from valid slice references
            // 3. SSE2 intrinsics marked with #[target_feature]
            // 4. Unaligned loads/stores handle unaligned data correctly
            i += 1;
        }
    }

    #[inline]
    #[target_feature(enable = "sse2")]
    // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    // 1. Loop bounds ensure `i + 4 <= len` before calling `.add(i)`
    // 2. All pointers derived from valid slice references
    // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    // 4. Unaligned loads/stores handle unaligned data correctly
    unsafe fn recip(a: &[f32], result: &mut [f32]) {
        let len = a.len();
        let mut i = 0;

        // Process 4 elements at a time with SIMD
        // Note: _mm_rcp_ps is an approximation with ~12-bit precision
        // For exact results, we use division
        let one = _mm_set1_ps(1.0);
        while i + 4 <= len {
            let vec = _mm_loadu_ps(a.as_ptr().add(i));
            let recip_vec = _mm_div_ps(one, vec);
            _mm_storeu_ps(result.as_mut_ptr().add(i), recip_vec);
            i += 4;
        }

        // Handle remaining elements
        while i < len {
            result[i] = a[i].recip();
            i += 1;
        }
    }

    unsafe fn ln(a: &[f32], result: &mut [f32]) {
        // Scalar fallback: SIMD transcendental functions require polynomial approximations
        super::scalar::ScalarBackend::ln(a, result);
    }

    unsafe fn log2(a: &[f32], result: &mut [f32]) {
        // Scalar fallback: SIMD transcendental functions require polynomial approximations
        super::scalar::ScalarBackend::log2(a, result);
    }

    unsafe fn log10(a: &[f32], result: &mut [f32]) {
        // Scalar fallback: SIMD transcendental functions require polynomial approximations
        super::scalar::ScalarBackend::log10(a, result);
    }

    unsafe fn sin(a: &[f32], result: &mut [f32]) {
        // Scalar fallback: SIMD transcendental functions require polynomial approximations
        super::scalar::ScalarBackend::sin(a, result);
    }

    unsafe fn cos(a: &[f32], result: &mut [f32]) {
        // Scalar fallback: SIMD transcendental functions require polynomial approximations
        super::scalar::ScalarBackend::cos(a, result);
    }

    unsafe fn tan(a: &[f32], result: &mut [f32]) {
        // Scalar fallback: SIMD transcendental functions require polynomial approximations
        super::scalar::ScalarBackend::tan(a, result);
    }

    unsafe fn floor(a: &[f32], result: &mut [f32]) {
        // Scalar fallback: floor requires SSE4.1 (not available in SSE2)
        super::scalar::ScalarBackend::floor(a, result);
    }

    unsafe fn ceil(a: &[f32], result: &mut [f32]) {
        // Scalar fallback: ceil requires SSE4.1 (not available in SSE2)
        super::scalar::ScalarBackend::ceil(a, result);
    }

    unsafe fn round(a: &[f32], result: &mut [f32]) {
        // Scalar fallback: round requires SSE4.1 (not available in SSE2)
        super::scalar::ScalarBackend::round(a, result);
    }
}

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

    #[test]
    fn test_sse2_add() {
        let a = [1.0, 2.0, 3.0, 4.0, 5.0];
        let b = [5.0, 6.0, 7.0, 8.0, 9.0];
        let mut result = [0.0; 5];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            Sse2Backend::add(&a, &b, &mut result);
        }

        assert_eq!(result, [6.0, 8.0, 10.0, 12.0, 14.0]);
    }

    #[test]
    fn test_sse2_mul() {
        let a = [1.0, 2.0, 3.0, 4.0, 5.0];
        let b = [2.0, 3.0, 4.0, 5.0, 6.0];
        let mut result = [0.0; 5];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            Sse2Backend::mul(&a, &b, &mut result);
        }

        assert_eq!(result, [2.0, 6.0, 12.0, 20.0, 30.0]);
    }

    #[test]
    fn test_sse2_dot() {
        let a = [1.0, 2.0, 3.0, 4.0];
        let b = [4.0, 5.0, 6.0, 7.0];

        // SAFETY: Test code calling backend trait methods marked unsafe
        let result = unsafe { Sse2Backend::dot(&a, &b) };

        assert_eq!(result, 60.0); // 1*4 + 2*5 + 3*6 + 4*7 = 60
    }

    #[test]
    fn test_sse2_sum() {
        let a = [1.0, 2.0, 3.0, 4.0, 5.0];
        // SAFETY: Test code calling backend trait methods marked unsafe
        let result = unsafe { Sse2Backend::sum(&a) };
        assert_eq!(result, 15.0);
    }

    #[test]
    fn test_sse2_max() {
        let a = [1.0, 5.0, 3.0, 2.0, 4.0];
        // SAFETY: Test code calling backend trait methods marked unsafe
        let result = unsafe { Sse2Backend::max(&a) };
        assert_eq!(result, 5.0);
    }

    #[test]
    fn test_sse2_min() {
        let a = [1.0, 5.0, 3.0, 2.0, 4.0];
        // SAFETY: Test code calling backend trait methods marked unsafe
        let result = unsafe { Sse2Backend::min(&a) };
        assert_eq!(result, 1.0);
    }

    #[test]
    fn test_sse2_matches_scalar() {
        // Verify SSE2 produces same results as scalar
        let a = [1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5];
        let b = [8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5];

        let mut scalar_result = [0.0; 7];
        let mut sse2_result = [0.0; 7];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            super::super::scalar::ScalarBackend::add(&a, &b, &mut scalar_result);
            Sse2Backend::add(&a, &b, &mut sse2_result);
        }

        assert_eq!(scalar_result, sse2_result);
    }

    #[test]
    fn test_sse2_relu() {
        let a = [-3.0, -1.0, 0.0, 1.0, 3.0, -2.0, 2.0, -0.5];
        let mut result = [0.0; 8];
        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            Sse2Backend::relu(&a, &mut result);
        }
        assert_eq!(result, [0.0, 0.0, 0.0, 1.0, 3.0, 0.0, 2.0, 0.0]);
    }

    #[test]
    fn test_sse2_relu_matches_scalar() {
        // Verify SSE2 relu produces same results as scalar
        let a = [-5.0, -3.0, -1.0, 0.0, 1.0, 3.0, 5.0];

        let mut scalar_result = [0.0; 7];
        let mut sse2_result = [0.0; 7];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            super::super::scalar::ScalarBackend::relu(&a, &mut scalar_result);
            Sse2Backend::relu(&a, &mut sse2_result);
        }

        assert_eq!(scalar_result, sse2_result);
    }

    #[test]
    fn test_sse2_sigmoid_matches_scalar() {
        // Verify SSE2 sigmoid produces same results as scalar
        let a = [-10.0, -1.0, 0.0, 1.0, 10.0];

        let mut scalar_result = [0.0; 5];
        let mut sse2_result = [0.0; 5];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            super::super::scalar::ScalarBackend::sigmoid(&a, &mut scalar_result);
            Sse2Backend::sigmoid(&a, &mut sse2_result);
        }

        for (s, e) in scalar_result.iter().zip(sse2_result.iter()) {
            assert!(
                (s - e).abs() < 1e-6,
                "sigmoid mismatch: scalar={}, sse2={}",
                s,
                e
            );
        }
    }

    #[cfg(target_arch = "x86_64")]
    #[test]
    fn test_sse2_exp_matches_scalar() {
        if !is_x86_feature_detected!("sse2") {
            eprintln!("Skipping SSE2 test: CPU does not support SSE2");
            return;
        }

        use super::super::scalar::ScalarBackend;

        // Test various ranges: negative, zero, positive, large values
        let test_values = vec![
            -10.0, -5.0, -2.0, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 50.0, -50.0, 87.0,
            -87.0, // near overflow/underflow limits
        ];
        let mut sse2_result = vec![0.0; test_values.len()];
        let mut scalar_result = vec![0.0; test_values.len()];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            Sse2Backend::exp(&test_values, &mut sse2_result);
            ScalarBackend::exp(&test_values, &mut scalar_result);
        }

        for (i, (sse2, scalar)) in sse2_result.iter().zip(scalar_result.iter()).enumerate() {
            let rel_error = if scalar.abs() > 1e-10 {
                (sse2 - scalar).abs() / scalar.abs()
            } else {
                (sse2 - scalar).abs()
            };
            assert!(
                rel_error < 1e-5,
                "exp({}) mismatch: sse2={}, scalar={}, rel_error={}",
                test_values[i],
                sse2,
                scalar,
                rel_error
            );
        }
    }

    #[test]
    fn test_sse2_gelu_matches_scalar() {
        // Verify SSE2 gelu produces same results as scalar
        let a = [-2.0, -1.0, 0.0, 1.0, 2.0];

        let mut scalar_result = [0.0; 5];
        let mut sse2_result = [0.0; 5];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            super::super::scalar::ScalarBackend::gelu(&a, &mut scalar_result);
            Sse2Backend::gelu(&a, &mut sse2_result);
        }

        for (s, e) in scalar_result.iter().zip(sse2_result.iter()) {
            assert!(
                (s - e).abs() < 1e-5,
                "gelu mismatch: scalar={}, sse2={}",
                s,
                e
            );
        }
    }

    #[test]
    fn test_sse2_swish_matches_scalar() {
        // Verify SSE2 swish produces same results as scalar
        let a = [-10.0, -1.0, 0.0, 1.0, 10.0];

        let mut scalar_result = [0.0; 5];
        let mut sse2_result = [0.0; 5];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            super::super::scalar::ScalarBackend::swish(&a, &mut scalar_result);
            Sse2Backend::swish(&a, &mut sse2_result);
        }

        for (s, e) in scalar_result.iter().zip(sse2_result.iter()) {
            assert!(
                (s - e).abs() < 1e-5,
                "swish mismatch: scalar={}, sse2={}",
                s,
                e
            );
        }
    }

    #[test]
    fn test_sse2_sub_matches_scalar() {
        let a = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0];
        let b = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0];

        let mut scalar_result = [0.0; 7];
        let mut sse2_result = [0.0; 7];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            super::super::scalar::ScalarBackend::sub(&a, &b, &mut scalar_result);
            Sse2Backend::sub(&a, &b, &mut sse2_result);
        }

        assert_eq!(scalar_result, sse2_result);
    }

    #[test]
    fn test_sse2_div_matches_scalar() {
        let a = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0];
        let b = [2.0, 4.0, 5.0, 8.0, 10.0, 12.0, 14.0];

        let mut scalar_result = [0.0; 7];
        let mut sse2_result = [0.0; 7];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            super::super::scalar::ScalarBackend::div(&a, &b, &mut scalar_result);
            Sse2Backend::div(&a, &b, &mut sse2_result);
        }

        // Use tolerance-based comparison since rcp+refinement has ~5e-7 relative error
        for (i, (&s, &sse2)) in scalar_result.iter().zip(sse2_result.iter()).enumerate() {
            let rel_error = ((s - sse2) / s).abs();
            assert!(
                rel_error < 1e-5,
                "Div mismatch at index {}: scalar={}, sse2={}, rel_error={}",
                i,
                s,
                sse2,
                rel_error
            );
        }
    }

    #[test]
    fn test_sse2_scale_matches_scalar() {
        let a = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0];
        let scalar = 2.5;

        let mut scalar_result = [0.0; 7];
        let mut sse2_result = [0.0; 7];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            super::super::scalar::ScalarBackend::scale(&a, scalar, &mut scalar_result);
            Sse2Backend::scale(&a, scalar, &mut sse2_result);
        }

        assert_eq!(scalar_result, sse2_result);
    }

    #[test]
    fn test_sse2_clamp_matches_scalar() {
        let a = [1.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0];

        let mut scalar_result = [0.0; 7];
        let mut sse2_result = [0.0; 7];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            super::super::scalar::ScalarBackend::clamp(&a, 5.0, 20.0, &mut scalar_result);
            Sse2Backend::clamp(&a, 5.0, 20.0, &mut sse2_result);
        }

        assert_eq!(scalar_result, sse2_result);
    }

    #[test]
    fn test_sse2_fma_matches_scalar() {
        let a = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0];
        let b = [2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0];
        let c = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0];

        let mut scalar_result = [0.0; 7];
        let mut sse2_result = [0.0; 7];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            super::super::scalar::ScalarBackend::fma(&a, &b, &c, &mut scalar_result);
            Sse2Backend::fma(&a, &b, &c, &mut sse2_result);
        }

        assert_eq!(scalar_result, sse2_result);
    }

    #[test]
    fn test_sse2_lerp_matches_scalar() {
        let a = [0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0];
        let b = [100.0, 110.0, 120.0, 130.0, 140.0, 150.0, 160.0];

        let mut scalar_result = [0.0; 7];
        let mut sse2_result = [0.0; 7];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            super::super::scalar::ScalarBackend::lerp(&a, &b, 0.25, &mut scalar_result);
            Sse2Backend::lerp(&a, &b, 0.25, &mut sse2_result);
        }

        for (s, e) in scalar_result.iter().zip(sse2_result.iter()) {
            assert!(
                (s - e).abs() < 1e-5,
                "lerp mismatch: scalar={}, sse2={}",
                s,
                e
            );
        }
    }

    #[test]
    fn test_sse2_argmax_matches_scalar() {
        let a = [1.0, 5.0, 3.0, 10.0, 2.0, 8.0, 4.0];

        // SAFETY: Test code calling backend trait methods marked unsafe
        let scalar_result = unsafe { super::super::scalar::ScalarBackend::argmax(&a) };
        // SAFETY: CPU feature verified at runtime, slices bounds-checked
        let sse2_result = unsafe { Sse2Backend::argmax(&a) };

        assert_eq!(scalar_result, sse2_result);
    }

    #[test]
    fn test_sse2_argmin_matches_scalar() {
        let a = [5.0, 1.0, 3.0, 10.0, 2.0, 8.0, 4.0];

        // SAFETY: Test code calling backend trait methods marked unsafe
        let scalar_result = unsafe { super::super::scalar::ScalarBackend::argmin(&a) };
        // SAFETY: CPU feature verified at runtime, slices bounds-checked
        let sse2_result = unsafe { Sse2Backend::argmin(&a) };

        assert_eq!(scalar_result, sse2_result);
    }

    #[test]
    fn test_sse2_sum_kahan_matches_scalar() {
        let a = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0];

        // SAFETY: Test code calling backend trait methods marked unsafe
        let scalar_result = unsafe { super::super::scalar::ScalarBackend::sum_kahan(&a) };
        // SAFETY: CPU feature verified at runtime, slices bounds-checked
        let sse2_result = unsafe { Sse2Backend::sum_kahan(&a) };

        assert!((scalar_result - sse2_result).abs() < 1e-5);
    }

    #[test]
    fn test_sse2_norm_l1_matches_scalar() {
        let a = [1.0, -2.0, 3.0, -4.0, 5.0, -6.0, 7.0];

        // SAFETY: Test code calling backend trait methods marked unsafe
        let scalar_result = unsafe { super::super::scalar::ScalarBackend::norm_l1(&a) };
        // SAFETY: CPU feature verified at runtime, slices bounds-checked
        let sse2_result = unsafe { Sse2Backend::norm_l1(&a) };

        assert!((scalar_result - sse2_result).abs() < 1e-5);
    }

    #[test]
    fn test_sse2_norm_l2_matches_scalar() {
        let a = [3.0, 4.0, 0.0, 0.0, 5.0, 12.0, 0.0];

        // SAFETY: Test code calling backend trait methods marked unsafe
        let scalar_result = unsafe { super::super::scalar::ScalarBackend::norm_l2(&a) };
        // SAFETY: CPU feature verified at runtime, slices bounds-checked
        let sse2_result = unsafe { Sse2Backend::norm_l2(&a) };

        assert!((scalar_result - sse2_result).abs() < 1e-5);
    }

    #[test]
    fn test_sse2_dot_matches_scalar() {
        let a = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0];
        let b = [7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0];

        // SAFETY: Test code calling backend trait methods marked unsafe
        let scalar_result = unsafe { super::super::scalar::ScalarBackend::dot(&a, &b) };
        // SAFETY: CPU feature verified at runtime, slices bounds-checked
        let sse2_result = unsafe { Sse2Backend::dot(&a, &b) };

        assert!((scalar_result - sse2_result).abs() < 1e-5);
    }

    #[test]
    fn test_sse2_mul_matches_scalar() {
        let a = [1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5];
        let b = [2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0];

        let mut scalar_result = [0.0; 7];
        let mut sse2_result = [0.0; 7];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            super::super::scalar::ScalarBackend::mul(&a, &b, &mut scalar_result);
            Sse2Backend::mul(&a, &b, &mut sse2_result);
        }

        assert_eq!(scalar_result, sse2_result);
    }

    #[test]
    fn test_sse2_add_matches_scalar() {
        let a = [1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5];
        let b = [8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5];

        let mut scalar_result = [0.0; 7];
        let mut sse2_result = [0.0; 7];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            super::super::scalar::ScalarBackend::add(&a, &b, &mut scalar_result);
            Sse2Backend::add(&a, &b, &mut sse2_result);
        }

        assert_eq!(scalar_result, sse2_result);
    }

    #[test]
    fn test_sse2_sum_matches_scalar() {
        let a = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0];

        // SAFETY: Test code calling backend trait methods marked unsafe
        let scalar_result = unsafe { super::super::scalar::ScalarBackend::sum(&a) };
        // SAFETY: CPU feature verified at runtime, slices bounds-checked
        let sse2_result = unsafe { Sse2Backend::sum(&a) };

        assert!((scalar_result - sse2_result).abs() < 1e-5);
    }

    #[test]
    fn test_sse2_max_matches_scalar() {
        let a = [1.0, 5.0, 3.0, 7.0, 2.0, 8.0, 4.0];

        // SAFETY: Test code calling backend trait methods marked unsafe
        let scalar_result = unsafe { super::super::scalar::ScalarBackend::max(&a) };
        // SAFETY: CPU feature verified at runtime, slices bounds-checked
        let sse2_result = unsafe { Sse2Backend::max(&a) };

        assert_eq!(scalar_result, sse2_result);
    }

    #[test]
    fn test_sse2_min_matches_scalar() {
        let a = [5.0, 1.0, 3.0, 7.0, 2.0, 8.0, 4.0];

        // SAFETY: Test code calling backend trait methods marked unsafe
        let scalar_result = unsafe { super::super::scalar::ScalarBackend::min(&a) };
        // SAFETY: CPU feature verified at runtime, slices bounds-checked
        let sse2_result = unsafe { Sse2Backend::min(&a) };

        assert_eq!(scalar_result, sse2_result);
    }

    #[test]
    fn test_sse2_tanh_matches_scalar() {
        // Verify SSE2 tanh produces same results as scalar
        let a = [-10.0, -1.0, 0.0, 1.0, 10.0];

        let mut scalar_result = [0.0; 5];
        let mut sse2_result = [0.0; 5];

        // SAFETY: Test code calling backend trait methods marked unsafe
        unsafe {
            super::super::scalar::ScalarBackend::tanh(&a, &mut scalar_result);
            Sse2Backend::tanh(&a, &mut sse2_result);
        }

        for (s, e) in scalar_result.iter().zip(sse2_result.iter()) {
            assert!(
                (s - e).abs() < 1e-5,
                "tanh mismatch: scalar={}, sse2={}",
                s,
                e
            );
        }
    }

    #[test]
    fn test_sse2_norm_linf_matches_scalar() {
        // Verify SSE2 norm_linf produces same results as scalar
        let test_cases = vec![
            vec![],                                       // empty
            vec![5.0],                                    // single element
            vec![-3.0, 1.0, -4.0, 1.0, 5.0],              // various values
            vec![-10.0, 5.0, 3.0, 7.0, -2.0, 8.0, 4.0],   // 7 elements (remainder)
            vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0], // 8 elements (aligned)
        ];

        for test_vec in test_cases {
            // SAFETY: Test code calling backend trait methods marked unsafe
            let scalar_result =
                // SAFETY: CPU feature verified at runtime, slices bounds-checked
                unsafe { super::super::scalar::ScalarBackend::norm_linf(&test_vec) };
            // SAFETY: CPU feature verified at runtime, slices bounds-checked
            let sse2_result = unsafe { Sse2Backend::norm_linf(&test_vec) };

            assert!(
                (scalar_result - sse2_result).abs() < 1e-5,
                "norm_linf mismatch for {:?}: scalar={}, sse2={}",
                test_vec,
                scalar_result,
                sse2_result
            );
        }
    }

    #[test]
    fn test_sse2_abs_matches_scalar() {
        // Verify SSE2 abs produces same results as scalar
        let test_cases = vec![
            vec![],                                             // empty
            vec![-5.0],                                         // single negative
            vec![5.0],                                          // single positive
            vec![-3.0, 1.0, -4.0, 1.5],                         // 4 elements (aligned)
            vec![-3.0, 1.0, -4.0, 1.5, -9.0, 2.0, -6.0],        // 7 elements (remainder)
            vec![-1.0, 2.0, -3.0, 4.0, -5.0, 6.0, -7.0, 8.0],   // 8 elements
            vec![0.0, -0.0, f32::INFINITY, f32::NEG_INFINITY],  // special values
        ];

        for test_vec in test_cases {
            let mut scalar_result = vec![0.0f32; test_vec.len()];
            let mut sse2_result = vec![0.0f32; test_vec.len()];

            // SAFETY: CPU feature verified at runtime, slices bounds-checked
            unsafe {
                super::super::scalar::ScalarBackend::abs(&test_vec, &mut scalar_result);
                Sse2Backend::abs(&test_vec, &mut sse2_result);
            }

            for (i, (&s, &e)) in scalar_result.iter().zip(sse2_result.iter()).enumerate() {
                // Handle NaN comparison
                if s.is_nan() && e.is_nan() {
                    continue;
                }
                assert!(
                    (s - e).abs() < 1e-5 || (s.is_infinite() && e.is_infinite() && s.signum() == e.signum()),
                    "abs mismatch at index {} for {:?}: scalar={}, sse2={}",
                    i,
                    test_vec,
                    s,
                    e
                );
            }
        }
    }

    #[test]
    fn test_sse2_tanh_saturation() {
        // Test tanh saturation at extreme values (scalar remainder path)
        let extreme_values = vec![-100.0, -50.0, -31.0, 31.0, 50.0, 100.0];
        let mut result = vec![0.0f32; extreme_values.len()];

        // SAFETY: CPU feature verified at runtime
        unsafe {
            Sse2Backend::tanh(&extreme_values, &mut result);
        }

        // Values < -30 should saturate to -1.0
        assert!((result[0] - (-1.0)).abs() < 1e-5, "tanh(-100) should be -1.0");
        assert!((result[1] - (-1.0)).abs() < 1e-5, "tanh(-50) should be -1.0");
        assert!((result[2] - (-1.0)).abs() < 1e-5, "tanh(-31) should be -1.0");

        // Values > 30 should saturate to 1.0
        assert!((result[3] - 1.0).abs() < 1e-5, "tanh(31) should be 1.0");
        assert!((result[4] - 1.0).abs() < 1e-5, "tanh(50) should be 1.0");
        assert!((result[5] - 1.0).abs() < 1e-5, "tanh(100) should be 1.0");
    }

    #[test]
    fn test_sse2_gelu_edge_cases() {
        // Test GELU with edge case values
        let values = vec![-10.0, -5.0, 0.0, 5.0, 10.0];
        let mut result = vec![0.0f32; values.len()];

        // SAFETY: CPU feature verified at runtime
        unsafe {
            Sse2Backend::gelu(&values, &mut result);
        }

        // GELU(-10) ≈ 0 (very negative values)
        assert!(result[0].abs() < 1e-3, "GELU(-10) should be near 0");

        // GELU(0) = 0
        assert!(result[2].abs() < 1e-5, "GELU(0) should be 0");

        // GELU(10) ≈ 10 (very positive values)
        assert!((result[4] - 10.0).abs() < 0.1, "GELU(10) should be near 10");
    }

    #[test]
    fn test_sse2_sigmoid_edge_cases() {
        // Test sigmoid saturation
        let values = vec![-100.0, -20.0, 0.0, 20.0, 100.0];
        let mut result = vec![0.0f32; values.len()];

        // SAFETY: CPU feature verified at runtime
        unsafe {
            Sse2Backend::sigmoid(&values, &mut result);
        }

        // sigmoid(-100) ≈ 0
        assert!(result[0] < 1e-5, "sigmoid(-100) should be near 0");

        // sigmoid(0) = 0.5
        assert!((result[2] - 0.5).abs() < 1e-5, "sigmoid(0) should be 0.5");

        // sigmoid(100) ≈ 1
        assert!((result[4] - 1.0).abs() < 1e-5, "sigmoid(100) should be near 1");
    }

    #[test]
    fn test_sse2_exp_edge_cases() {
        // Test exp with edge values that exercise saturation paths
        let values = vec![-100.0, -50.0, 0.0, 50.0, 88.0];
        let mut result = vec![0.0f32; values.len()];

        // SAFETY: CPU feature verified at runtime
        unsafe {
            Sse2Backend::exp(&values, &mut result);
        }

        // exp(-100) ≈ 0
        assert!(result[0] < 1e-30, "exp(-100) should be near 0");

        // exp(0) = 1
        assert!((result[2] - 1.0).abs() < 1e-5, "exp(0) should be 1");

        // exp(88) is large but finite
        assert!(result[4].is_finite(), "exp(88) should be finite");
    }
}

Coverage Report

Created: 2026-01-23 22:55