## Coverage Report

#### Created: 2026-01-23 22:55

next uncovered line (L), next uncovered region (R), next uncovered branch (B)

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

    Line
    Count
    Source
    1
    //! SSE2 backend implementation (x86_64 baseline SIMD)
    2
    //!
    3
    //! This backend uses SSE2 intrinsics for 128-bit SIMD operations.
    4
    //! SSE2 is available on all x86_64 CPUs as a baseline requirement.
    5
    //!
    6
    //! # Performance
    7
    //!
    8
    //! Expected speedup: 4x for operations on aligned f32 vectors (4 elements per register)
    9
    //!
    10
    //! # Safety
    11
    //!
    12
    //! All SSE2 intrinsics are marked `unsafe` by Rust. This module carefully isolates
    13
    //! all unsafe code and verifies correctness through comprehensive testing.
    14
    15
    #[cfg(target_arch = "x86_64")]
    16
    use std::arch::x86_64::*;
    17
    18
    use super::VectorBackend;
    19
    20
    /// SSE2 backend (128-bit SIMD for x86_64)
    21
    pub struct Sse2Backend;
    22
    23
    impl VectorBackend for Sse2Backend {
    24
        #[inline]
    25
        #[target_feature(enable = "sse2")]
    26
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    27
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    28
        // 2. All pointers derived from valid slice references with sufficient backing storage
    29
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    30
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    31
    0
        unsafe fn add(a: &[f32], b: &[f32], result: &mut [f32]) {
    32
    0
            let len = a.len();
    33
    0
            let mut i = 0;
    34
    35
            // Process 4 elements at a time using SSE2 (128-bit = 4 x f32)
    36
    0
            while i + 4 <= len {
    37
    0
                // Load 4 floats from a and b
    38
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    39
    0
                let vb = _mm_loadu_ps(b.as_ptr().add(i));
    40
    0
    41
    0
                // Add them
    42
    0
                let vresult = _mm_add_ps(va, vb);
    43
    0
    44
    0
                // Store result
    45
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);
    46
    0
    47
    0
                i += 4;
    48
    0
            }
    49
    50
            // Handle remaining elements with scalar code
    51
    0
            for j in i..len {
    52
    0
                result[j] = a[j] + b[j];
    53
    0
            }
    54
    0
        }
    55
    56
        #[inline]
    57
        #[target_feature(enable = "sse2")]
    58
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    59
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    60
        // 2. All pointers derived from valid slice references with sufficient backing storage
    61
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    62
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    63
    0
        unsafe fn sub(a: &[f32], b: &[f32], result: &mut [f32]) {
    64
    0
            let len = a.len();
    65
    0
            let mut i = 0;
    66
    67
            // Process 4 elements at a time using SSE2 (128-bit = 4 x f32)
    68
    0
            while i + 4 <= len {
    69
    0
                // Load 4 floats from a and b
    70
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    71
    0
                let vb = _mm_loadu_ps(b.as_ptr().add(i));
    72
    0
    73
    0
                // Subtract them
    74
    0
                let vresult = _mm_sub_ps(va, vb);
    75
    0
    76
    0
                // Store result
    77
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);
    78
    0
    79
    0
                i += 4;
    80
    0
            }
    81
    82
            // Handle remaining elements with scalar code
    83
    0
            for j in i..len {
    84
    0
                result[j] = a[j] - b[j];
    85
    0
            }
    86
    0
        }
    87
    88
        #[inline]
    89
        #[target_feature(enable = "sse2")]
    90
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    91
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    92
        // 2. All pointers derived from valid slice references with sufficient backing storage
    93
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    94
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    95
    0
        unsafe fn mul(a: &[f32], b: &[f32], result: &mut [f32]) {
    96
    0
            let len = a.len();
    97
    0
            let mut i = 0;
    98
    99
            // Process 4 elements at a time
    100
    0
            while i + 4 <= len {
    101
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    102
    0
                let vb = _mm_loadu_ps(b.as_ptr().add(i));
    103
    0
                let vresult = _mm_mul_ps(va, vb);
    104
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);
    105
    0
                i += 4;
    106
    0
            }
    107
    108
            // Handle remaining elements
    109
    0
            for j in i..len {
    110
    0
                result[j] = a[j] * b[j];
    111
    0
            }
    112
    0
        }
    113
    114
        #[inline]
    115
        #[target_feature(enable = "sse2")]
    116
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    117
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    118
        // 2. All pointers derived from valid slice references with sufficient backing storage
    119
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    120
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    121
    0
        unsafe fn div(a: &[f32], b: &[f32], result: &mut [f32]) {
    122
    0
            let len = a.len();
    123
    0
            let mut i = 0;
    124
    125
            // Use direct division instruction
    126
            // Previous reciprocal approximation + Newton-Raphson added too much overhead
    127
            // for small workloads. Direct divps is simpler and performs better overall.
    128
            //
    129
            // Performance (measured 2025-11-21):
    130
            // - 100 elem: reciprocal=90.4ns, direct=~expected 80-85ns
    131
            // - 1000 elem: reciprocal=295ns (1.08x), direct expected similar
    132
            // - Trade-off: Simpler code, better small workload performance
    133
    134
            // Process 4 elements at a time
    135
    0
            while i + 4 <= len {
    136
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    137
    0
                let vb = _mm_loadu_ps(b.as_ptr().add(i));
    138
    0
    139
    0
                // Direct division (13-14 cycle latency, but simpler)
    140
    0
                let vresult = _mm_div_ps(va, vb);
    141
    0
    142
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);
    143
    0
                i += 4;
    144
    0
            }
    145
    146
            // Handle remaining elements
    147
    0
            for j in i..len {
    148
    0
                result[j] = a[j] / b[j];
    149
    0
            }
    150
    0
        }
    151
    152
        #[inline]
    153
        #[target_feature(enable = "sse2")]
    154
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    155
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    156
        // 2. All pointers derived from valid slice references with sufficient backing storage
    157
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    158
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    159
    0
        unsafe fn dot(a: &[f32], b: &[f32]) -> f32 {
    160
    0
            let len = a.len();
    161
    0
            let mut i = 0;
    162
    163
            // Accumulator for SIMD portion
    164
    0
            let mut sum_vec = _mm_setzero_ps();
    165
    166
            // Process 4 elements at a time
    167
    0
            while i + 4 <= len {
    168
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    169
    0
                let vb = _mm_loadu_ps(b.as_ptr().add(i));
    170
    0
                let vmul = _mm_mul_ps(va, vb);
    171
    0
                sum_vec = _mm_add_ps(sum_vec, vmul);
    172
    0
                i += 4;
    173
    0
            }
    174
    175
            // Horizontal sum using faster movehl/shuffle pattern
    176
    0
            let mut sum = {
    177
    0
                let temp = _mm_add_ps(sum_vec, _mm_movehl_ps(sum_vec, sum_vec));
    178
    0
                let temp = _mm_add_ss(temp, _mm_shuffle_ps(temp, temp, 1));
    179
    0
                _mm_cvtss_f32(temp)
    180
            };
    181
    182
            // Handle remaining elements with scalar code
    183
    0
            for j in i..len {
    184
    0
                sum += a[j] * b[j];
    185
    0
            }
    186
    187
    0
            sum
    188
    0
        }
    189
    190
        #[inline]
    191
        #[target_feature(enable = "sse2")]
    192
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    193
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    194
        // 2. All pointers derived from valid slice references with sufficient backing storage
    195
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    196
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    197
    0
        unsafe fn sum(a: &[f32]) -> f32 {
    198
    0
            let len = a.len();
    199
    0
            let mut i = 0;
    200
    0
            let mut sum_vec = _mm_setzero_ps();
    201
    202
            // Process 4 elements at a time
    203
    0
            while i + 4 <= len {
    204
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    205
    0
                sum_vec = _mm_add_ps(sum_vec, va);
    206
    0
                i += 4;
    207
    0
            }
    208
    209
            // Horizontal sum using faster movehl/shuffle pattern
    210
    0
            let mut sum = {
    211
    0
                let temp = _mm_add_ps(sum_vec, _mm_movehl_ps(sum_vec, sum_vec));
    212
    0
                let temp = _mm_add_ss(temp, _mm_shuffle_ps(temp, temp, 1));
    213
    0
                _mm_cvtss_f32(temp)
    214
            };
    215
    216
            // Handle remaining elements
    217
    0
            sum += a[i..len].iter().sum::<f32>();
    218
    219
    0
            sum
    220
    0
        }
    221
    222
        #[inline]
    223
        #[target_feature(enable = "sse2")]
    224
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    225
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    226
        // 2. All pointers derived from valid slice references with sufficient backing storage
    227
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    228
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    229
    0
        unsafe fn max(a: &[f32]) -> f32 {
    230
    0
            let len = a.len();
    231
    0
            let mut i = 0;
    232
    233
            // Initialize with first element broadcast to all lanes
    234
    0
            let mut max_vec = _mm_set1_ps(a[0]);
    235
    236
            // Process 4 elements at a time
    237
    0
            while i + 4 <= len {
    238
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    239
    0
                max_vec = _mm_max_ps(max_vec, va);
    240
    0
                i += 4;
    241
    0
            }
    242
    243
            // Horizontal max using faster movehl/shuffle pattern
    244
    0
            let mut maximum = {
    245
    0
                let temp = _mm_max_ps(max_vec, _mm_movehl_ps(max_vec, max_vec));
    246
    0
                let temp = _mm_max_ss(temp, _mm_shuffle_ps(temp, temp, 1));
    247
    0
                _mm_cvtss_f32(temp)
    248
            };
    249
    250
            // Handle remaining elements
    251
    0
            for &val in &a[i..len] {
    252
    0
                if val > maximum {
    253
    0
                    maximum = val;
    254
    0
                }
    255
            }
    256
    257
    0
            maximum
    258
    0
        }
    259
    260
        #[inline]
    261
        #[target_feature(enable = "sse2")]
    262
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    263
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    264
        // 2. All pointers derived from valid slice references with sufficient backing storage
    265
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    266
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    267
    0
        unsafe fn min(a: &[f32]) -> f32 {
    268
    0
            let len = a.len();
    269
    0
            let mut i = 0;
    270
    271
            // Initialize with first element broadcast to all lanes
    272
    0
            let mut min_vec = _mm_set1_ps(a[0]);
    273
    274
            // Process 4 elements at a time
    275
    0
            while i + 4 <= len {
    276
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    277
    0
                min_vec = _mm_min_ps(min_vec, va);
    278
    0
                i += 4;
    279
    0
            }
    280
    281
            // Horizontal min using faster movehl/shuffle pattern
    282
    0
            let mut minimum = {
    283
    0
                let temp = _mm_min_ps(min_vec, _mm_movehl_ps(min_vec, min_vec));
    284
    0
                let temp = _mm_min_ss(temp, _mm_shuffle_ps(temp, temp, 1));
    285
    0
                _mm_cvtss_f32(temp)
    286
            };
    287
    288
            // Handle remaining elements
    289
    0
            for &val in &a[i..len] {
    290
    0
                if val < minimum {
    291
    0
                    minimum = val;
    292
    0
                }
    293
            }
    294
    295
    0
            minimum
    296
    0
        }
    297
    298
        #[inline]
    299
        #[target_feature(enable = "sse2")]
    300
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    301
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    302
        // 2. All pointers derived from valid slice references with sufficient backing storage
    303
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    304
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    305
    0
        unsafe fn argmax(a: &[f32]) -> usize {
    306
    0
            let len = a.len();
    307
    0
            let mut i = 0;
    308
    309
            // Initialize SIMD vectors with first element value and index 0
    310
    0
            let mut vmax = _mm_set1_ps(a[0]);
    311
    0
            let mut vmax_idx = _mm_set1_ps(0.0); // Track indices as floats
    312
    313
            // Initialize index vector [0, 1, 2, 3] and increment constant
    314
    0
            let mut vidx_current = _mm_set_ps(3.0, 2.0, 1.0, 0.0);
    315
    0
            let vinc = _mm_set1_ps(4.0);
    316
    317
            // Process 4 elements at a time with index tracking
    318
    0
            while i + 4 <= len {
    319
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    320
    0
    321
    0
                // Compare: va > vmax (strict greater-than to preserve first occurrence)
    322
    0
                let mask = _mm_cmpgt_ps(va, vmax);
    323
    0
    324
    0
                // Conditionally update max values and indices using SSE2 blend emulation
    325
    0
                // blend = (mask & new) | (~mask & old)
    326
    0
                vmax = _mm_or_ps(_mm_and_ps(mask, va), _mm_andnot_ps(mask, vmax));
    327
    0
                vmax_idx = _mm_or_ps(
    328
    0
                    _mm_and_ps(mask, vidx_current),
    329
    0
                    _mm_andnot_ps(mask, vmax_idx),
    330
    0
                );
    331
    0
    332
    0
                // Increment index vector for next iteration
    333
    0
                vidx_current = _mm_add_ps(vidx_current, vinc);
    334
    0
                i += 4;
    335
    0
            }
    336
    337
            // Horizontal reduction: find max and its index across 4 lanes
    338
    0
            let mut max_array = [0.0f32; 4];
    339
    0
            let mut idx_array = [0.0f32; 4];
    340
    0
            _mm_storeu_ps(max_array.as_mut_ptr(), vmax);
    341
    0
            _mm_storeu_ps(idx_array.as_mut_ptr(), vmax_idx);
    342
    343
    0
            let mut max_value = max_array[0];
    344
    0
            let mut max_index = idx_array[0] as usize;
    345
    0
            for j in 1..4 {
    346
    0
                if max_array[j] > max_value {
    347
    0
                    max_value = max_array[j];
    348
    0
                    max_index = idx_array[j] as usize;
    349
    0
                }
    350
            }
    351
    352
            // Handle remaining elements
    353
    0
            for (idx, &val) in a[i..].iter().enumerate() {
    354
    0
                if val > max_value {
    355
    0
                    max_value = val;
    356
    0
                    max_index = i + idx;
    357
    0
                }
    358
            }
    359
    360
    0
            max_index
    361
    0
        }
    362
    363
        #[inline]
    364
        #[target_feature(enable = "sse2")]
    365
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    366
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    367
        // 2. All pointers derived from valid slice references with sufficient backing storage
    368
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    369
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    370
    0
        unsafe fn argmin(a: &[f32]) -> usize {
    371
    0
            let len = a.len();
    372
    0
            let mut i = 0;
    373
    374
            // Initialize SIMD vectors with first element value and index 0
    375
    0
            let mut vmin = _mm_set1_ps(a[0]);
    376
    0
            let mut vmin_idx = _mm_set1_ps(0.0); // Track indices as floats
    377
    378
            // Initialize index vector [0, 1, 2, 3] and increment constant
    379
    0
            let mut vidx_current = _mm_set_ps(3.0, 2.0, 1.0, 0.0);
    380
    0
            let vinc = _mm_set1_ps(4.0);
    381
    382
            // Process 4 elements at a time with index tracking
    383
    0
            while i + 4 <= len {
    384
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    385
    0
    386
    0
                // Compare: va < vmin (strict less-than to preserve first occurrence)
    387
    0
                let mask = _mm_cmplt_ps(va, vmin);
    388
    0
    389
    0
                // Conditionally update min values and indices using SSE2 blend emulation
    390
    0
                // blend = (mask & new) | (~mask & old)
    391
    0
                vmin = _mm_or_ps(_mm_and_ps(mask, va), _mm_andnot_ps(mask, vmin));
    392
    0
                vmin_idx = _mm_or_ps(
    393
    0
                    _mm_and_ps(mask, vidx_current),
    394
    0
                    _mm_andnot_ps(mask, vmin_idx),
    395
    0
                );
    396
    0
    397
    0
                // Increment index vector for next iteration
    398
    0
                vidx_current = _mm_add_ps(vidx_current, vinc);
    399
    0
                i += 4;
    400
    0
            }
    401
    402
            // Horizontal reduction: find min and its index across 4 lanes
    403
    0
            let mut min_array = [0.0f32; 4];
    404
    0
            let mut idx_array = [0.0f32; 4];
    405
    0
            _mm_storeu_ps(min_array.as_mut_ptr(), vmin);
    406
    0
            _mm_storeu_ps(idx_array.as_mut_ptr(), vmin_idx);
    407
    408
    0
            let mut min_value = min_array[0];
    409
    0
            let mut min_index = idx_array[0] as usize;
    410
    0
            for j in 1..4 {
    411
    0
                if min_array[j] < min_value {
    412
    0
                    min_value = min_array[j];
    413
    0
                    min_index = idx_array[j] as usize;
    414
    0
                }
    415
            }
    416
    417
            // Handle remaining elements
    418
    0
            for (idx, &val) in a[i..].iter().enumerate() {
    419
    0
                if val < min_value {
    420
    0
                    min_value = val;
    421
    0
                    min_index = i + idx;
    422
    0
                }
    423
            }
    424
    425
    0
            min_index
    426
    0
        }
    427
    428
        #[inline]
    429
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    430
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    431
        // 2. All pointers derived from valid slice references with sufficient backing storage
    432
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    433
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    434
    0
        unsafe fn sum_kahan(a: &[f32]) -> f32 {
    435
            // Kahan summation is inherently sequential, use scalar implementation
    436
    0
            super::scalar::ScalarBackend::sum_kahan(a)
    437
    0
        }
    438
    439
        #[inline]
    440
        #[target_feature(enable = "sse2")]
    441
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    442
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    443
        // 2. All pointers derived from valid slice references with sufficient backing storage
    444
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    445
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    446
    0
        unsafe fn norm_l2(a: &[f32]) -> f32 {
    447
    0
            if a.is_empty() {
    448
    0
                return 0.0;
    449
    0
            }
    450
    451
            // L2 norm is sqrt(dot(a, a))
    452
    0
            let sum_of_squares = Self::dot(a, a);
    453
    0
            sum_of_squares.sqrt()
    454
    0
        }
    455
    456
        #[inline]
    457
        #[target_feature(enable = "sse2")]
    458
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    459
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    460
        // 2. All pointers derived from valid slice references with sufficient backing storage
    461
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    462
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    463
    0
        unsafe fn norm_l1(a: &[f32]) -> f32 {
    464
    0
            if a.is_empty() {
    465
    0
                return 0.0;
    466
    0
            }
    467
    468
    0
            let len = a.len();
    469
    0
            let mut i = 0;
    470
    471
            // Accumulator for 4-way parallel accumulation
    472
    0
            let mut acc = _mm_setzero_ps();
    473
    474
            // SSE2 doesn't have abs for floats, use bitwise AND to clear sign bit
    475
    0
            let sign_mask = _mm_set1_ps(f32::from_bits(0x7FFF_FFFF));
    476
    477
            // Process 4 elements at a time
    478
    0
            while i + 4 <= len {
    479
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    480
    0
    481
    0
                // Compute absolute value by clearing sign bit
    482
    0
                let abs_va = _mm_and_ps(va, sign_mask);
    483
    0
    484
    0
                // Accumulate
    485
    0
                acc = _mm_add_ps(acc, abs_va);
    486
    0
    487
    0
                i += 4;
    488
    0
            }
    489
    490
            // Horizontal sum: extract all 4 lanes and sum them
    491
    0
            let mut result = {
    492
    0
                let temp = _mm_add_ps(acc, _mm_movehl_ps(acc, acc));
    493
    0
                let temp = _mm_add_ss(temp, _mm_shuffle_ps(temp, temp, 1));
    494
    0
                _mm_cvtss_f32(temp)
    495
            };
    496
    497
            // Handle remaining elements with scalar code
    498
    0
            for &val in &a[i..] {
    499
    0
                result += val.abs();
    500
    0
            }
    501
    502
    0
            result
    503
    0
        }
    504
    505
        #[inline]
    506
        #[target_feature(enable = "sse2")]
    507
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    508
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    509
        // 2. All pointers derived from valid slice references with sufficient backing storage
    510
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    511
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    512
    0
        unsafe fn norm_linf(a: &[f32]) -> f32 {
    513
    0
            if a.is_empty() {
    514
    0
                return 0.0;
    515
    0
            }
    516
    517
    0
            let len = a.len();
    518
    0
            let mut i = 0;
    519
    520
            // Accumulator for maximum value
    521
    0
            let mut max_vec = _mm_setzero_ps();
    522
    523
            // SSE2 doesn't have abs for floats, use bitwise AND to clear sign bit
    524
    0
            let sign_mask = _mm_set1_ps(f32::from_bits(0x7FFF_FFFF));
    525
    526
            // Process 4 elements at a time
    527
    0
            while i + 4 <= len {
    528
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    529
    0
                // Compute absolute value
    530
    0
                let abs_va = _mm_and_ps(va, sign_mask);
    531
    0
                // Update maximum
    532
    0
                max_vec = _mm_max_ps(max_vec, abs_va);
    533
    0
                i += 4;
    534
    0
            }
    535
    536
            // Horizontal max: extract all 4 lanes and take maximum
    537
    0
            let mut result = {
    538
    0
                let temp = _mm_max_ps(max_vec, _mm_movehl_ps(max_vec, max_vec));
    539
    0
                let temp = _mm_max_ss(temp, _mm_shuffle_ps(temp, temp, 1));
    540
    0
                _mm_cvtss_f32(temp)
    541
            };
    542
    543
            // Handle remaining elements with scalar code
    544
    0
            for &val in &a[i..] {
    545
    0
                let abs_val = val.abs();
    546
    0
                if abs_val > result {
    547
    0
                    result = abs_val;
    548
    0
                }
    549
            }
    550
    551
    0
            result
    552
    0
        }
    553
    554
        #[inline]
    555
        #[target_feature(enable = "sse2")]
    556
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    557
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    558
        // 2. All pointers derived from valid slice references with sufficient backing storage
    559
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    560
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    561
    0
        unsafe fn scale(a: &[f32], scalar: f32, result: &mut [f32]) {
    562
    0
            let len = a.len();
    563
    0
            let mut i = 0;
    564
    565
            // Broadcast scalar to all 4 lanes
    566
    0
            let scalar_vec = _mm_set1_ps(scalar);
    567
    568
            // Process 4 elements at a time
    569
    0
            while i + 4 <= len {
    570
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    571
    0
                let vresult = _mm_mul_ps(va, scalar_vec);
    572
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);
    573
    0
                i += 4;
    574
    0
            }
    575
    576
            // Handle remaining elements
    577
    0
            while i < len {
    578
    0
                result[i] = a[i] * scalar;
    579
    0
                i += 1;
    580
    0
            }
    581
    0
        }
    582
    583
        #[inline]
    584
        #[target_feature(enable = "sse2")]
    585
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    586
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    587
        // 2. All pointers derived from valid slice references with sufficient backing storage
    588
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    589
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    590
    0
        unsafe fn abs(a: &[f32], result: &mut [f32]) {
    591
    0
            let len = a.len();
    592
    0
            let mut i = 0;
    593
    594
            // Create mask to clear sign bit (0x7FFFFFFF for all elements)
    595
    0
            let sign_mask = _mm_set1_ps(f32::from_bits(0x7FFF_FFFF));
    596
    597
            // Process 4 elements at a time using SSE2 (128-bit = 4 x f32)
    598
    0
            while i + 4 <= len {
    599
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    600
    0
    601
    0
                // Compute absolute value by clearing sign bit
    602
    0
                let abs_va = _mm_and_ps(va, sign_mask);
    603
    0
    604
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), abs_va);
    605
    0
                i += 4;
    606
    0
            }
    607
    608
            // Handle remaining elements with scalar code
    609
    0
            while i < len {
    610
    0
                result[i] = a[i].abs();
    611
    0
                i += 1;
    612
    0
            }
    613
    0
        }
    614
    615
        #[inline]
    616
        #[target_feature(enable = "sse2")]
    617
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    618
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    619
        // 2. All pointers derived from valid slice references with sufficient backing storage
    620
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    621
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    622
    0
        unsafe fn clamp(a: &[f32], min_val: f32, max_val: f32, result: &mut [f32]) {
    623
    0
            let len = a.len();
    624
    0
            let mut i = 0;
    625
    626
            // Broadcast min and max to all 4 lanes
    627
    0
            let min_vec = _mm_set1_ps(min_val);
    628
    0
            let max_vec = _mm_set1_ps(max_val);
    629
    630
            // Process 4 elements at a time
    631
    0
            while i + 4 <= len {
    632
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    633
    0
                let clamped = _mm_min_ps(_mm_max_ps(va, min_vec), max_vec);
    634
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), clamped);
    635
    0
                i += 4;
    636
    0
            }
    637
    638
            // Handle remaining elements
    639
    0
            while i < len {
    640
    0
                result[i] = a[i].max(min_val).min(max_val);
    641
    0
                i += 1;
    642
    0
            }
    643
    0
        }
    644
    645
        #[inline]
    646
        #[target_feature(enable = "sse2")]
    647
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    648
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    649
        // 2. All pointers derived from valid slice references with sufficient backing storage
    650
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    651
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    652
    0
        unsafe fn lerp(a: &[f32], b: &[f32], t: f32, result: &mut [f32]) {
    653
    0
            let len = a.len();
    654
    0
            let mut i = 0;
    655
    656
            // Broadcast t to all 4 lanes
    657
    0
            let t_vec = _mm_set1_ps(t);
    658
    659
            // Process 4 elements at a time
    660
    0
            while i + 4 <= len {
    661
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    662
    0
                let vb = _mm_loadu_ps(b.as_ptr().add(i));
    663
    0
    664
    0
                // result = a + t * (b - a)
    665
    0
                let diff = _mm_sub_ps(vb, va);
    666
    0
                let scaled_diff = _mm_mul_ps(t_vec, diff);
    667
    0
                let vresult = _mm_add_ps(va, scaled_diff);
    668
    0
    669
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);
    670
    0
                i += 4;
    671
    0
            }
    672
    673
            // Handle remaining elements
    674
    0
            while i < len {
    675
    0
                result[i] = a[i] + t * (b[i] - a[i]);
    676
    0
                i += 1;
    677
    0
            }
    678
    0
        }
    679
    680
        #[inline]
    681
        #[target_feature(enable = "sse2")]
    682
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    683
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    684
        // 2. All pointers derived from valid slice references with sufficient backing storage
    685
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    686
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    687
    0
        unsafe fn fma(a: &[f32], b: &[f32], c: &[f32], result: &mut [f32]) {
    688
    0
            let len = a.len();
    689
    0
            let mut i = 0;
    690
    691
            // Process 4 elements at a time
    692
    0
            while i + 4 <= len {
    693
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    694
    0
                let vb = _mm_loadu_ps(b.as_ptr().add(i));
    695
    0
                let vc = _mm_loadu_ps(c.as_ptr().add(i));
    696
    0
    697
    0
                // result = a * b + c
    698
    0
                // SSE2 doesn't have FMA, so we use separate mul and add
    699
    0
                let product = _mm_mul_ps(va, vb);
    700
    0
                let vresult = _mm_add_ps(product, vc);
    701
    0
    702
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);
    703
    0
                i += 4;
    704
    0
            }
    705
    706
            // Handle remaining elements
    707
    0
            while i < len {
    708
    0
                result[i] = a[i] * b[i] + c[i];
    709
    0
                i += 1;
    710
    0
            }
    711
    0
        }
    712
    713
        #[inline]
    714
        #[target_feature(enable = "sse2")]
    715
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    716
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    717
        // 2. All pointers derived from valid slice references with sufficient backing storage
    718
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    719
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    720
    0
        unsafe fn relu(a: &[f32], result: &mut [f32]) {
    721
    0
            let len = a.len();
    722
    0
            let mut i = 0;
    723
    724
            // Zero vector for max comparison
    725
    0
            let zero = _mm_setzero_ps();
    726
    727
            // Process 4 elements at a time
    728
    0
            while i + 4 <= len {
    729
    0
                let va = _mm_loadu_ps(a.as_ptr().add(i));
    730
    0
    731
    0
                // ReLU: max(0, x)
    732
    0
                let vresult = _mm_max_ps(zero, va);
    733
    0
    734
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);
    735
    0
                i += 4;
    736
    0
            }
    737
    738
            // Handle remaining elements
    739
    0
            while i < len {
    740
    0
                result[i] = if a[i] > 0.0 { a[i] } else { 0.0 };
    741
    0
                i += 1;
    742
            }
    743
    0
        }
    744
    745
        #[inline]
    746
        #[target_feature(enable = "sse2")]
    747
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    748
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    749
        // 2. All pointers derived from valid slice references with sufficient backing storage
    750
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    751
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    752
    0
        unsafe fn exp(a: &[f32], result: &mut [f32]) {
    753
    0
            let len = a.len();
    754
    0
            let mut i = 0;
    755
    756
            // Constants for range reduction: exp(x) = 2^(x * log2(e)) = 2^k * 2^r
    757
    0
            let log2e = _mm_set1_ps(std::f32::consts::LOG2_E); // 1.442695...
    758
    0
            let ln2 = _mm_set1_ps(std::f32::consts::LN_2); // 0.693147...
    759
    0
            let half = _mm_set1_ps(0.5);
    760
    0
            let one = _mm_set1_ps(1.0);
    761
    762
            // Polynomial coefficients for e^r approximation (Remez minimax on [-ln(2)/2, ln(2)/2])
    763
            // e^r ≈ 1 + c1*r + c2*r^2 + c3*r^3 + c4*r^4 + c5*r^5 + c6*r^6
    764
            // Coefficients from Cephes/SLEEF libraries optimized for f32
    765
    0
            let c1 = _mm_set1_ps(1.0);
    766
    0
            let c2 = _mm_set1_ps(0.5);
    767
    0
            let c3 = _mm_set1_ps(0.166_666_67); // 1/6
    768
    0
            let c4 = _mm_set1_ps(0.041_666_668); // 1/24
    769
    0
            let c5 = _mm_set1_ps(0.008_333_334); // 1/120
    770
    0
            let c6 = _mm_set1_ps(0.001_388_889); // 1/720
    771
    772
            // Limits for overflow/underflow handling
    773
    0
            let exp_hi = _mm_set1_ps(88.376_26); // ln(FLT_MAX)
    774
    0
            let exp_lo = _mm_set1_ps(-87.336_55); // ln(FLT_MIN) approximately
    775
    776
            // Process 4 elements at a time
    777
    0
            while i + 4 <= len {
    778
    0
                let x = _mm_loadu_ps(a.as_ptr().add(i));
    779
    0
    780
    0
                // Clamp x to avoid overflow/underflow
    781
    0
                let x = _mm_max_ps(_mm_min_ps(x, exp_hi), exp_lo);
    782
    0
    783
    0
                // Range reduction: x' = x * log2(e), then k = round(x'), r = x' - k
    784
    0
                let x_scaled = _mm_mul_ps(x, log2e);
    785
    0
    786
    0
                // k = round(x_scaled) = floor(x_scaled + 0.5)
    787
    0
                // SSE2 floor emulation: convert to int (truncates toward zero), then convert back
    788
    0
                let k_plus_half = _mm_add_ps(x_scaled, half);
    789
    0
                let k_int = _mm_cvttps_epi32(k_plus_half); // truncate toward zero
    790
    0
                let k = _mm_cvtepi32_ps(k_int);
    791
    0
                // Adjust for negative numbers: if k > k_plus_half, subtract 1
    792
    0
                let mask = _mm_cmpgt_ps(k, k_plus_half);
    793
    0
                let k = _mm_sub_ps(k, _mm_and_ps(mask, one));
    794
    0
    795
    0
                // r = x - k * ln(2) (in original base e space)
    796
    0
                let r = _mm_sub_ps(x, _mm_mul_ps(k, ln2));
    797
    0
    798
    0
                // Polynomial approximation: e^r ≈ 1 + c1*r + c2*r^2 + c3*r^3 + c4*r^4 + c5*r^5 + c6*r^6
    799
    0
                // Use Horner's method: ((((((c6*r + c5)*r + c4)*r + c3)*r + c2)*r + c1)*r + 1)
    800
    0
                // No FMA in SSE2, so use mul + add
    801
    0
                let mut p = c6;
    802
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c5);
    803
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c4);
    804
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c3);
    805
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c2);
    806
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c1);
    807
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), one);
    808
    0
    809
    0
                // Scale by 2^k using IEEE754 exponent manipulation
    810
    0
                // 2^k is computed by adding k to the exponent bits
    811
    0
                let k_int = _mm_cvtps_epi32(k);
    812
    0
                let k_shifted = _mm_slli_epi32(k_int, 23); // shift to exponent position
    813
    0
                let scale = _mm_castsi128_ps(_mm_add_epi32(_mm_castps_si128(one), k_shifted));
    814
    0
    815
    0
                // Final result: e^x = e^r * 2^k
    816
    0
                let vresult = _mm_mul_ps(p, scale);
    817
    0
    818
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), vresult);
    819
    0
                i += 4;
    820
    0
            }
    821
    822
            // Handle remaining elements with scalar code
    823
    0
            while i < len {
    824
    0
                result[i] = a[i].exp();
    825
    0
                i += 1;
    826
    0
            }
    827
    0
        }
    828
    829
        #[inline]
    830
        #[target_feature(enable = "sse2")]
    831
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    832
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    833
        // 2. All pointers derived from valid slice references with sufficient backing storage
    834
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    835
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    836
    0
        unsafe fn sigmoid(a: &[f32], result: &mut [f32]) {
    837
            // sigmoid(x) = 1 / (1 + exp(-x))
    838
            // Use SIMD exp approximation with range reduction
    839
    0
            let len = a.len();
    840
    0
            let mut i = 0;
    841
    842
            // Constants for exp(-x) computation
    843
    0
            let log2e = _mm_set1_ps(std::f32::consts::LOG2_E);
    844
    0
            let ln2 = _mm_set1_ps(std::f32::consts::LN_2);
    845
    0
            let half = _mm_set1_ps(0.5);
    846
    0
            let one = _mm_set1_ps(1.0);
    847
    848
            // Taylor series coefficients for e^r
    849
    0
            let c1 = _mm_set1_ps(1.0);
    850
    0
            let c2 = _mm_set1_ps(0.5);
    851
    0
            let c3 = _mm_set1_ps(0.166_666_67);
    852
    0
            let c4 = _mm_set1_ps(0.041_666_668);
    853
    0
            let c5 = _mm_set1_ps(0.008_333_334);
    854
    0
            let c6 = _mm_set1_ps(0.001_388_889);
    855
    856
            // Limits for overflow/underflow
    857
    0
            let exp_hi = _mm_set1_ps(88.376_26);
    858
    0
            let exp_lo = _mm_set1_ps(-87.336_55);
    859
    860
            // Process 4 elements at a time
    861
    0
            while i + 4 <= len {
    862
    0
                let x = _mm_loadu_ps(a.as_ptr().add(i));
    863
    0
    864
    0
                // Compute -x for exp(-x)
    865
    0
                let neg_x = _mm_sub_ps(_mm_setzero_ps(), x);
    866
    0
    867
    0
                // Clamp to avoid overflow/underflow
    868
    0
                let neg_x = _mm_max_ps(_mm_min_ps(neg_x, exp_hi), exp_lo);
    869
    0
    870
    0
                // Range reduction: exp(-x) computation
    871
    0
                let x_scaled = _mm_mul_ps(neg_x, log2e);
    872
    0
    873
    0
                // SSE2 floor emulation
    874
    0
                let k_plus_half = _mm_add_ps(x_scaled, half);
    875
    0
                let k_int = _mm_cvttps_epi32(k_plus_half);
    876
    0
                let k = _mm_cvtepi32_ps(k_int);
    877
    0
                let mask = _mm_cmpgt_ps(k, k_plus_half);
    878
    0
                let k = _mm_sub_ps(k, _mm_and_ps(mask, one));
    879
    0
    880
    0
                let r = _mm_sub_ps(neg_x, _mm_mul_ps(k, ln2));
    881
    0
    882
    0
                // Polynomial approximation using Horner's method (no FMA in SSE2)
    883
    0
                let mut p = c6;
    884
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c5);
    885
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c4);
    886
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c3);
    887
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c2);
    888
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c1);
    889
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), one);
    890
    0
    891
    0
                // Scale by 2^k
    892
    0
                let k_int = _mm_cvtps_epi32(k);
    893
    0
                let k_shifted = _mm_slli_epi32(k_int, 23);
    894
    0
                let scale = _mm_castsi128_ps(_mm_add_epi32(_mm_castps_si128(one), k_shifted));
    895
    0
                let exp_neg_x = _mm_mul_ps(p, scale);
    896
    0
    897
    0
                // sigmoid = 1 / (1 + exp(-x))
    898
    0
                let denom = _mm_add_ps(one, exp_neg_x);
    899
    0
                let sigmoid_result = _mm_div_ps(one, denom);
    900
    0
    901
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), sigmoid_result);
    902
    0
                i += 4;
    903
    0
            }
    904
    905
            // Handle remaining elements with scalar code
    906
    0
            while i < len {
    907
    0
                let val = a[i];
    908
    0
                result[i] = if val < -50.0 {
    909
    0
                    0.0
    910
    0
                } else if val > 50.0 {
    911
    0
                    1.0
    912
                } else {
    913
    0
                    1.0 / (1.0 + (-val).exp())
    914
                };
    915
    0
                i += 1;
    916
            }
    917
    0
        }
    918
    919
        #[inline]
    920
        #[target_feature(enable = "sse2")]
    921
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    922
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    923
        // 2. All pointers derived from valid slice references with sufficient backing storage
    924
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    925
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    926
    0
        unsafe fn gelu(a: &[f32], result: &mut [f32]) {
    927
            // gelu(x) = 0.5 * x * (1 + tanh(sqrt(2/π) * (x + 0.044715 * x³)))
    928
            // Use SIMD tanh via: tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
    929
    0
            let len = a.len();
    930
    0
            let mut i = 0;
    931
    932
            // GELU constants
    933
    0
            let sqrt_2_over_pi = _mm_set1_ps(0.797_884_6);
    934
    0
            let coeff = _mm_set1_ps(0.044715);
    935
    0
            let half = _mm_set1_ps(0.5);
    936
    0
            let one = _mm_set1_ps(1.0);
    937
    0
            let two = _mm_set1_ps(2.0);
    938
    939
            // Constants for exp computation
    940
    0
            let log2e = _mm_set1_ps(std::f32::consts::LOG2_E);
    941
    0
            let ln2 = _mm_set1_ps(std::f32::consts::LN_2);
    942
    943
            // Taylor series coefficients for e^r
    944
    0
            let c1 = _mm_set1_ps(1.0);
    945
    0
            let c2 = _mm_set1_ps(0.5);
    946
    0
            let c3 = _mm_set1_ps(0.166_666_67);
    947
    0
            let c4 = _mm_set1_ps(0.041_666_668);
    948
    0
            let c5 = _mm_set1_ps(0.008_333_334);
    949
    0
            let c6 = _mm_set1_ps(0.001_388_889);
    950
    951
            // Limits for overflow/underflow
    952
    0
            let exp_hi = _mm_set1_ps(88.376_26);
    953
    0
            let exp_lo = _mm_set1_ps(-87.336_55);
    954
    955
            // Process 4 elements at a time
    956
    0
            while i + 4 <= len {
    957
    0
                let x = _mm_loadu_ps(a.as_ptr().add(i));
    958
    0
    959
    0
                // Compute inner = sqrt(2/π) * (x + 0.044715 * x³)
    960
    0
                let x2 = _mm_mul_ps(x, x);
    961
    0
                let x3 = _mm_mul_ps(x2, x);
    962
    0
                let inner_sum = _mm_add_ps(x, _mm_mul_ps(coeff, x3));
    963
    0
                let inner = _mm_mul_ps(sqrt_2_over_pi, inner_sum);
    964
    0
    965
    0
                // Compute tanh(inner) = (exp(2*inner) - 1) / (exp(2*inner) + 1)
    966
    0
                let two_inner = _mm_mul_ps(two, inner);
    967
    0
    968
    0
                // Clamp to avoid overflow/underflow
    969
    0
                let two_inner = _mm_max_ps(_mm_min_ps(two_inner, exp_hi), exp_lo);
    970
    0
    971
    0
                // Range reduction for exp(2*inner)
    972
    0
                let x_scaled = _mm_mul_ps(two_inner, log2e);
    973
    0
    974
    0
                // SSE2 floor emulation
    975
    0
                let k_plus_half = _mm_add_ps(x_scaled, half);
    976
    0
                let k_int = _mm_cvttps_epi32(k_plus_half);
    977
    0
                let k = _mm_cvtepi32_ps(k_int);
    978
    0
                let mask = _mm_cmpgt_ps(k, k_plus_half);
    979
    0
                let k = _mm_sub_ps(k, _mm_and_ps(mask, one));
    980
    0
    981
    0
                let r = _mm_sub_ps(two_inner, _mm_mul_ps(k, ln2));
    982
    0
    983
    0
                // Polynomial approximation using Horner's method (no FMA in SSE2)
    984
    0
                let mut p = c6;
    985
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c5);
    986
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c4);
    987
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c3);
    988
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c2);
    989
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c1);
    990
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), one);
    991
    0
    992
    0
                // Scale by 2^k
    993
    0
                let k_int = _mm_cvtps_epi32(k);
    994
    0
                let k_shifted = _mm_slli_epi32(k_int, 23);
    995
    0
                let scale = _mm_castsi128_ps(_mm_add_epi32(_mm_castps_si128(one), k_shifted));
    996
    0
                let exp_2inner = _mm_mul_ps(p, scale);
    997
    0
    998
    0
                // tanh = (exp(2x) - 1) / (exp(2x) + 1)
    999
    0
                let tanh_numer = _mm_sub_ps(exp_2inner, one);
    1000
    0
                let tanh_denom = _mm_add_ps(exp_2inner, one);
    1001
    0
                let tanh_result = _mm_div_ps(tanh_numer, tanh_denom);
    1002
    0
    1003
    0
                // gelu = 0.5 * x * (1 + tanh)
    1004
    0
                let one_plus_tanh = _mm_add_ps(one, tanh_result);
    1005
    0
                let gelu_result = _mm_mul_ps(half, _mm_mul_ps(x, one_plus_tanh));
    1006
    0
    1007
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), gelu_result);
    1008
    0
                i += 4;
    1009
    0
            }
    1010
    1011
            // Handle remaining elements with scalar code
    1012
            const SQRT_2_OVER_PI: f32 = 0.797_884_6;
    1013
            const COEFF: f32 = 0.044715;
    1014
    1015
    0
            while i < len {
    1016
    0
                let x = a[i];
    1017
    0
                let x3 = x * x * x;
    1018
    0
                let inner = SQRT_2_OVER_PI * (x + COEFF * x3);
    1019
    0
                result[i] = 0.5 * x * (1.0 + inner.tanh());
    1020
    0
                i += 1;
    1021
    0
            }
    1022
    0
        }
    1023
    1024
        #[inline]
    1025
        #[target_feature(enable = "sse2")]
    1026
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    1027
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    1028
        // 2. All pointers derived from valid slice references with sufficient backing storage
    1029
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    1030
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    1031
    0
        unsafe fn swish(a: &[f32], result: &mut [f32]) {
    1032
            // swish(x) = x * sigmoid(x) = x / (1 + exp(-x))
    1033
            // Use SIMD exp approximation with range reduction
    1034
    0
            let len = a.len();
    1035
    0
            let mut i = 0;
    1036
    1037
            // Constants for exp(-x) computation
    1038
    0
            let log2e = _mm_set1_ps(std::f32::consts::LOG2_E);
    1039
    0
            let ln2 = _mm_set1_ps(std::f32::consts::LN_2);
    1040
    0
            let half = _mm_set1_ps(0.5);
    1041
    0
            let one = _mm_set1_ps(1.0);
    1042
    1043
            // Taylor series coefficients for e^r
    1044
    0
            let c1 = _mm_set1_ps(1.0);
    1045
    0
            let c2 = _mm_set1_ps(0.5);
    1046
    0
            let c3 = _mm_set1_ps(0.166_666_67);
    1047
    0
            let c4 = _mm_set1_ps(0.041_666_668);
    1048
    0
            let c5 = _mm_set1_ps(0.008_333_334);
    1049
    0
            let c6 = _mm_set1_ps(0.001_388_889);
    1050
    1051
            // Limits for overflow/underflow
    1052
    0
            let exp_hi = _mm_set1_ps(88.376_26);
    1053
    0
            let exp_lo = _mm_set1_ps(-87.336_55);
    1054
    1055
            // Process 4 elements at a time
    1056
    0
            while i + 4 <= len {
    1057
    0
                let x = _mm_loadu_ps(a.as_ptr().add(i));
    1058
    0
    1059
    0
                // Compute -x for exp(-x)
    1060
    0
                let neg_x = _mm_sub_ps(_mm_setzero_ps(), x);
    1061
    0
    1062
    0
                // Clamp to avoid overflow/underflow
    1063
    0
                let neg_x = _mm_max_ps(_mm_min_ps(neg_x, exp_hi), exp_lo);
    1064
    0
    1065
    0
                // Range reduction: exp(-x) computation
    1066
    0
                let x_scaled = _mm_mul_ps(neg_x, log2e);
    1067
    0
    1068
    0
                // SSE2 floor emulation
    1069
    0
                let k_plus_half = _mm_add_ps(x_scaled, half);
    1070
    0
                let k_int = _mm_cvttps_epi32(k_plus_half);
    1071
    0
                let k = _mm_cvtepi32_ps(k_int);
    1072
    0
                let mask = _mm_cmpgt_ps(k, k_plus_half);
    1073
    0
                let k = _mm_sub_ps(k, _mm_and_ps(mask, one));
    1074
    0
    1075
    0
                let r = _mm_sub_ps(neg_x, _mm_mul_ps(k, ln2));
    1076
    0
    1077
    0
                // Polynomial approximation using Horner's method (no FMA in SSE2)
    1078
    0
                let mut p = c6;
    1079
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c5);
    1080
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c4);
    1081
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c3);
    1082
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c2);
    1083
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c1);
    1084
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), one);
    1085
    0
    1086
    0
                // Scale by 2^k
    1087
    0
                let k_int = _mm_cvtps_epi32(k);
    1088
    0
                let k_shifted = _mm_slli_epi32(k_int, 23);
    1089
    0
                let scale = _mm_castsi128_ps(_mm_add_epi32(_mm_castps_si128(one), k_shifted));
    1090
    0
                let exp_neg_x = _mm_mul_ps(p, scale);
    1091
    0
    1092
    0
                // swish = x / (1 + exp(-x))
    1093
    0
                let denom = _mm_add_ps(one, exp_neg_x);
    1094
    0
                let swish_result = _mm_div_ps(x, denom);
    1095
    0
    1096
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), swish_result);
    1097
    0
                i += 4;
    1098
    0
            }
    1099
    1100
            // Handle remaining elements with scalar code
    1101
    0
            while i < len {
    1102
    0
                let x = a[i];
    1103
    0
                result[i] = if x < -50.0 {
    1104
    0
                    0.0
    1105
    0
                } else if x > 50.0 {
    1106
    0
                    x
    1107
                } else {
    1108
    0
                    x / (1.0 + (-x).exp())
    1109
                };
    1110
    0
                i += 1;
    1111
            }
    1112
    0
        }
    1113
    1114
        #[inline]
    1115
        #[target_feature(enable = "sse2")]
    1116
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    1117
        // 1. Loop bounds ensure `i + N <= len` before calling `.add(i)` (N=4 for SSE2)
    1118
        // 2. All pointers derived from valid slice references with sufficient backing storage
    1119
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    1120
        // 4. Unaligned loads/stores used (_mm_loadu_ps/_mm_storeu_ps) - no alignment requirement
    1121
    0
        unsafe fn tanh(a: &[f32], result: &mut [f32]) {
    1122
            // tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
    1123
            // Use SIMD exp approximation with range reduction
    1124
    0
            let len = a.len();
    1125
    0
            let mut i = 0;
    1126
    1127
            // Constants for exp(2x) computation
    1128
    0
            let log2e = _mm_set1_ps(std::f32::consts::LOG2_E);
    1129
    0
            let ln2 = _mm_set1_ps(std::f32::consts::LN_2);
    1130
    0
            let half = _mm_set1_ps(0.5);
    1131
    0
            let one = _mm_set1_ps(1.0);
    1132
    0
            let two = _mm_set1_ps(2.0);
    1133
    1134
            // Taylor series coefficients for e^r
    1135
    0
            let c1 = _mm_set1_ps(1.0);
    1136
    0
            let c2 = _mm_set1_ps(0.5);
    1137
    0
            let c3 = _mm_set1_ps(0.166_666_67);
    1138
    0
            let c4 = _mm_set1_ps(0.041_666_668);
    1139
    0
            let c5 = _mm_set1_ps(0.008_333_334);
    1140
    0
            let c6 = _mm_set1_ps(0.001_388_889);
    1141
    1142
            // Limits for overflow/underflow
    1143
    0
            let exp_hi = _mm_set1_ps(88.376_26);
    1144
    0
            let exp_lo = _mm_set1_ps(-87.336_55);
    1145
    1146
            // Process 4 elements at a time
    1147
    0
            while i + 4 <= len {
    1148
    0
                let x = _mm_loadu_ps(a.as_ptr().add(i));
    1149
    0
    1150
    0
                // Compute 2x for exp(2x)
    1151
    0
                let two_x = _mm_mul_ps(two, x);
    1152
    0
    1153
    0
                // Clamp to avoid overflow/underflow
    1154
    0
                let two_x = _mm_max_ps(_mm_min_ps(two_x, exp_hi), exp_lo);
    1155
    0
    1156
    0
                // Range reduction: exp(2x) computation
    1157
    0
                let x_scaled = _mm_mul_ps(two_x, log2e);
    1158
    0
    1159
    0
                // SSE2 floor emulation
    1160
    0
                let k_plus_half = _mm_add_ps(x_scaled, half);
    1161
    0
                let k_int = _mm_cvttps_epi32(k_plus_half);
    1162
    0
                let k = _mm_cvtepi32_ps(k_int);
    1163
    0
                let mask = _mm_cmpgt_ps(k, k_plus_half);
    1164
    0
                let k = _mm_sub_ps(k, _mm_and_ps(mask, one));
    1165
    0
    1166
    0
                let r = _mm_sub_ps(two_x, _mm_mul_ps(k, ln2));
    1167
    0
    1168
    0
                // Polynomial approximation using Horner's method (no FMA in SSE2)
    1169
    0
                let mut p = c6;
    1170
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c5);
    1171
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c4);
    1172
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c3);
    1173
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c2);
    1174
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), c1);
    1175
    0
                p = _mm_add_ps(_mm_mul_ps(p, r), one);
    1176
    0
    1177
    0
                // Scale by 2^k
    1178
    0
                let k_int = _mm_cvtps_epi32(k);
    1179
    0
                let k_shifted = _mm_slli_epi32(k_int, 23);
    1180
    0
                let scale = _mm_castsi128_ps(_mm_add_epi32(_mm_castps_si128(one), k_shifted));
    1181
    0
                let exp_2x = _mm_mul_ps(p, scale);
    1182
    0
    1183
    0
                // tanh = (exp(2x) - 1) / (exp(2x) + 1)
    1184
    0
                let tanh_numer = _mm_sub_ps(exp_2x, one);
    1185
    0
                let tanh_denom = _mm_add_ps(exp_2x, one);
    1186
    0
                let tanh_result = _mm_div_ps(tanh_numer, tanh_denom);
    1187
    0
    1188
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), tanh_result);
    1189
    0
                i += 4;
    1190
    0
            }
    1191
    1192
            // Handle remaining elements with scalar code
    1193
    0
            while i < len {
    1194
    0
                let x = a[i];
    1195
    0
                result[i] = if x < -30.0 {
    1196
    0
                    -1.0
    1197
    0
                } else if x > 30.0 {
    1198
    0
                    1.0
    1199
                } else {
    1200
    0
                    let exp_2x = (2.0 * x).exp();
    1201
    0
                    (exp_2x - 1.0) / (exp_2x + 1.0)
    1202
                };
    1203
    0
                i += 1;
    1204
            }
    1205
    0
        }
    1206
    1207
        #[inline]
    1208
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    1209
        // 1. Loop bounds ensure proper array access
    1210
        // 2. All pointers derived from valid slice references
    1211
        // 3. SSE2 intrinsics marked with #[target_feature]
    1212
        // 4. Unaligned loads/stores handle unaligned data correctly
    1213
        #[target_feature(enable = "sse2")]
    1214
    0
        unsafe fn sqrt(a: &[f32], result: &mut [f32]) {
    1215
    0
            let len = a.len();
    1216
    0
            let mut i = 0;
    1217
    1218
            // Process 4 elements at a time with SIMD
    1219
    0
            while i + 4 <= len {
    1220
    0
                let vec = _mm_loadu_ps(a.as_ptr().add(i));
    1221
    0
                let sqrt_vec = _mm_sqrt_ps(vec);
    1222
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), sqrt_vec);
    1223
    0
                i += 4;
    1224
    0
            }
    1225
    1226
            // Handle remaining elements
    1227
    0
            while i < len {
    1228
    0
                result[i] = a[i].sqrt();
    1229
    0
                // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    1230
    0
                // 1. Loop bounds ensure proper array access
    1231
    0
                // 2. All pointers derived from valid slice references
    1232
    0
                // 3. SSE2 intrinsics marked with #[target_feature]
    1233
    0
                // 4. Unaligned loads/stores handle unaligned data correctly
    1234
    0
                i += 1;
    1235
    0
            }
    1236
    0
        }
    1237
    1238
        #[inline]
    1239
        #[target_feature(enable = "sse2")]
    1240
        // SAFETY: Pointer arithmetic and SIMD intrinsics are safe because:
    1241
        // 1. Loop bounds ensure `i + 4 <= len` before calling `.add(i)`
    1242
        // 2. All pointers derived from valid slice references
    1243
        // 3. SSE2 intrinsics marked with #[target_feature(enable = "sse2")]
    1244
        // 4. Unaligned loads/stores handle unaligned data correctly
    1245
    0
        unsafe fn recip(a: &[f32], result: &mut [f32]) {
    1246
    0
            let len = a.len();
    1247
    0
            let mut i = 0;
    1248
    1249
            // Process 4 elements at a time with SIMD
    1250
            // Note: _mm_rcp_ps is an approximation with ~12-bit precision
    1251
            // For exact results, we use division
    1252
    0
            let one = _mm_set1_ps(1.0);
    1253
    0
            while i + 4 <= len {
    1254
    0
                let vec = _mm_loadu_ps(a.as_ptr().add(i));
    1255
    0
                let recip_vec = _mm_div_ps(one, vec);
    1256
    0
                _mm_storeu_ps(result.as_mut_ptr().add(i), recip_vec);
    1257
    0
                i += 4;
    1258
    0
            }
    1259
    1260
            // Handle remaining elements
    1261
    0
            while i < len {
    1262
    0
                result[i] = a[i].recip();
    1263
    0
                i += 1;
    1264
    0
            }
    1265
    0
        }
    1266
    1267
    0
        unsafe fn ln(a: &[f32], result: &mut [f32]) {
    1268
            // Scalar fallback: SIMD transcendental functions require polynomial approximations
    1269
    0
            super::scalar::ScalarBackend::ln(a, result);
    1270
    0
        }
    1271
    1272
    0
        unsafe fn log2(a: &[f32], result: &mut [f32]) {
    1273
            // Scalar fallback: SIMD transcendental functions require polynomial approximations
    1274
    0
            super::scalar::ScalarBackend::log2(a, result);
    1275
    0
        }
    1276
    1277
    0
        unsafe fn log10(a: &[f32], result: &mut [f32]) {
    1278
            // Scalar fallback: SIMD transcendental functions require polynomial approximations
    1279
    0
            super::scalar::ScalarBackend::log10(a, result);
    1280
    0
        }
    1281
    1282
    0
        unsafe fn sin(a: &[f32], result: &mut [f32]) {
    1283
            // Scalar fallback: SIMD transcendental functions require polynomial approximations
    1284
    0
            super::scalar::ScalarBackend::sin(a, result);
    1285
    0
        }
    1286
    1287
    0
        unsafe fn cos(a: &[f32], result: &mut [f32]) {
    1288
            // Scalar fallback: SIMD transcendental functions require polynomial approximations
    1289
    0
            super::scalar::ScalarBackend::cos(a, result);
    1290
    0
        }
    1291
    1292
    0
        unsafe fn tan(a: &[f32], result: &mut [f32]) {
    1293
            // Scalar fallback: SIMD transcendental functions require polynomial approximations
    1294
    0
            super::scalar::ScalarBackend::tan(a, result);
    1295
    0
        }
    1296
    1297
    0
        unsafe fn floor(a: &[f32], result: &mut [f32]) {
    1298
            // Scalar fallback: floor requires SSE4.1 (not available in SSE2)
    1299
    0
            super::scalar::ScalarBackend::floor(a, result);
    1300
    0
        }
    1301
    1302
    0
        unsafe fn ceil(a: &[f32], result: &mut [f32]) {
    1303
            // Scalar fallback: ceil requires SSE4.1 (not available in SSE2)
    1304
    0
            super::scalar::ScalarBackend::ceil(a, result);
    1305
    0
        }
    1306
    1307
    0
        unsafe fn round(a: &[f32], result: &mut [f32]) {
    1308
            // Scalar fallback: round requires SSE4.1 (not available in SSE2)
    1309
    0
            super::scalar::ScalarBackend::round(a, result);
    1310
    0
        }
    1311
    }
    1312
    1313
    #[cfg(test)]
    1314
    mod tests {
    1315
        use super::*;
    1316
    1317
        #[test]
    1318
        fn test_sse2_add() {
    1319
            let a = [1.0, 2.0, 3.0, 4.0, 5.0];
    1320
            let b = [5.0, 6.0, 7.0, 8.0, 9.0];
    1321
            let mut result = [0.0; 5];
    1322
    1323
            // SAFETY: Test code calling backend trait methods marked unsafe
    1324
            unsafe {
    1325
                Sse2Backend::add(&a, &b, &mut result);
    1326
            }
    1327
    1328
            assert_eq!(result, [6.0, 8.0, 10.0, 12.0, 14.0]);
    1329
        }
    1330
    1331
        #[test]
    1332
        fn test_sse2_mul() {
    1333
            let a = [1.0, 2.0, 3.0, 4.0, 5.0];
    1334
            let b = [2.0, 3.0, 4.0, 5.0, 6.0];
    1335
            let mut result = [0.0; 5];
    1336
    1337
            // SAFETY: Test code calling backend trait methods marked unsafe
    1338
            unsafe {
    1339
                Sse2Backend::mul(&a, &b, &mut result);
    1340
            }
    1341
    1342
            assert_eq!(result, [2.0, 6.0, 12.0, 20.0, 30.0]);
    1343
        }
    1344
    1345
        #[test]
    1346
        fn test_sse2_dot() {
    1347
            let a = [1.0, 2.0, 3.0, 4.0];
    1348
            let b = [4.0, 5.0, 6.0, 7.0];
    1349
    1350
            // SAFETY: Test code calling backend trait methods marked unsafe
    1351
            let result = unsafe { Sse2Backend::dot(&a, &b) };
    1352
    1353
            assert_eq!(result, 60.0); // 1*4 + 2*5 + 3*6 + 4*7 = 60
    1354
        }
    1355
    1356
        #[test]
    1357
        fn test_sse2_sum() {
    1358
            let a = [1.0, 2.0, 3.0, 4.0, 5.0];
    1359
            // SAFETY: Test code calling backend trait methods marked unsafe
    1360
            let result = unsafe { Sse2Backend::sum(&a) };
    1361
            assert_eq!(result, 15.0);
    1362
        }
    1363
    1364
        #[test]
    1365
        fn test_sse2_max() {
    1366
            let a = [1.0, 5.0, 3.0, 2.0, 4.0];
    1367
            // SAFETY: Test code calling backend trait methods marked unsafe
    1368
            let result = unsafe { Sse2Backend::max(&a) };
    1369
            assert_eq!(result, 5.0);
    1370
        }
    1371
    1372
        #[test]
    1373
        fn test_sse2_min() {
    1374
            let a = [1.0, 5.0, 3.0, 2.0, 4.0];
    1375
            // SAFETY: Test code calling backend trait methods marked unsafe
    1376
            let result = unsafe { Sse2Backend::min(&a) };
    1377
            assert_eq!(result, 1.0);
    1378
        }
    1379
    1380
        #[test]
    1381
        fn test_sse2_matches_scalar() {
    1382
            // Verify SSE2 produces same results as scalar
    1383
            let a = [1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5];
    1384
            let b = [8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5];
    1385
    1386
            let mut scalar_result = [0.0; 7];
    1387
            let mut sse2_result = [0.0; 7];
    1388
    1389
            // SAFETY: Test code calling backend trait methods marked unsafe
    1390
            unsafe {
    1391
                super::super::scalar::ScalarBackend::add(&a, &b, &mut scalar_result);
    1392
                Sse2Backend::add(&a, &b, &mut sse2_result);
    1393
            }
    1394
    1395
            assert_eq!(scalar_result, sse2_result);
    1396
        }
    1397
    1398
        #[test]
    1399
        fn test_sse2_relu() {
    1400
            let a = [-3.0, -1.0, 0.0, 1.0, 3.0, -2.0, 2.0, -0.5];
    1401
            let mut result = [0.0; 8];
    1402
            // SAFETY: Test code calling backend trait methods marked unsafe
    1403
            unsafe {
    1404
                Sse2Backend::relu(&a, &mut result);
    1405
            }
    1406
            assert_eq!(result, [0.0, 0.0, 0.0, 1.0, 3.0, 0.0, 2.0, 0.0]);
    1407
        }
    1408
    1409
        #[test]
    1410
        fn test_sse2_relu_matches_scalar() {
    1411
            // Verify SSE2 relu produces same results as scalar
    1412
            let a = [-5.0, -3.0, -1.0, 0.0, 1.0, 3.0, 5.0];
    1413
    1414
            let mut scalar_result = [0.0; 7];
    1415
            let mut sse2_result = [0.0; 7];
    1416
    1417
            // SAFETY: Test code calling backend trait methods marked unsafe
    1418
            unsafe {
    1419
                super::super::scalar::ScalarBackend::relu(&a, &mut scalar_result);
    1420
                Sse2Backend::relu(&a, &mut sse2_result);
    1421
            }
    1422
    1423
            assert_eq!(scalar_result, sse2_result);
    1424
        }
    1425
    1426
        #[test]
    1427
        fn test_sse2_sigmoid_matches_scalar() {
    1428
            // Verify SSE2 sigmoid produces same results as scalar
    1429
            let a = [-10.0, -1.0, 0.0, 1.0, 10.0];
    1430
    1431
            let mut scalar_result = [0.0; 5];
    1432
            let mut sse2_result = [0.0; 5];
    1433
    1434
            // SAFETY: Test code calling backend trait methods marked unsafe
    1435
            unsafe {
    1436
                super::super::scalar::ScalarBackend::sigmoid(&a, &mut scalar_result);
    1437
                Sse2Backend::sigmoid(&a, &mut sse2_result);
    1438
            }
    1439
    1440
            for (s, e) in scalar_result.iter().zip(sse2_result.iter()) {
    1441
                assert!(
    1442
                    (s - e).abs() < 1e-6,
    1443
                    "sigmoid mismatch: scalar={}, sse2={}",
    1444
                    s,
    1445
                    e
    1446
                );
    1447
            }
    1448
        }
    1449
    1450
        #[cfg(target_arch = "x86_64")]
    1451
        #[test]
    1452
        fn test_sse2_exp_matches_scalar() {
    1453
            if !is_x86_feature_detected!("sse2") {
    1454
                eprintln!("Skipping SSE2 test: CPU does not support SSE2");
    1455
                return;
    1456
            }
    1457
    1458
            use super::super::scalar::ScalarBackend;
    1459
    1460
            // Test various ranges: negative, zero, positive, large values
    1461
            let test_values = vec![
    1462
                -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,
    1463
                -87.0, // near overflow/underflow limits
    1464
            ];
    1465
            let mut sse2_result = vec![0.0; test_values.len()];
    1466
            let mut scalar_result = vec![0.0; test_values.len()];
    1467
    1468
            // SAFETY: Test code calling backend trait methods marked unsafe
    1469
            unsafe {
    1470
                Sse2Backend::exp(&test_values, &mut sse2_result);
    1471
                ScalarBackend::exp(&test_values, &mut scalar_result);
    1472
            }
    1473
    1474
            for (i, (sse2, scalar)) in sse2_result.iter().zip(scalar_result.iter()).enumerate() {
    1475
                let rel_error = if scalar.abs() > 1e-10 {
    1476
                    (sse2 - scalar).abs() / scalar.abs()
    1477
                } else {
    1478
                    (sse2 - scalar).abs()
    1479
                };
    1480
                assert!(
    1481
                    rel_error < 1e-5,
    1482
                    "exp({}) mismatch: sse2={}, scalar={}, rel_error={}",
    1483
                    test_values[i],
    1484
                    sse2,
    1485
                    scalar,
    1486
                    rel_error
    1487
                );
    1488
            }
    1489
        }
    1490
    1491
        #[test]
    1492
        fn test_sse2_gelu_matches_scalar() {
    1493
            // Verify SSE2 gelu produces same results as scalar
    1494
            let a = [-2.0, -1.0, 0.0, 1.0, 2.0];
    1495
    1496
            let mut scalar_result = [0.0; 5];
    1497
            let mut sse2_result = [0.0; 5];
    1498
    1499
            // SAFETY: Test code calling backend trait methods marked unsafe
    1500
            unsafe {
    1501
                super::super::scalar::ScalarBackend::gelu(&a, &mut scalar_result);
    1502
                Sse2Backend::gelu(&a, &mut sse2_result);
    1503
            }
    1504
    1505
            for (s, e) in scalar_result.iter().zip(sse2_result.iter()) {
    1506
                assert!(
    1507
                    (s - e).abs() < 1e-5,
    1508
                    "gelu mismatch: scalar={}, sse2={}",
    1509
                    s,
    1510
                    e
    1511
                );
    1512
            }
    1513
        }
    1514
    1515
        #[test]
    1516
        fn test_sse2_swish_matches_scalar() {
    1517
            // Verify SSE2 swish produces same results as scalar
    1518
            let a = [-10.0, -1.0, 0.0, 1.0, 10.0];
    1519
    1520
            let mut scalar_result = [0.0; 5];
    1521
            let mut sse2_result = [0.0; 5];
    1522
    1523
            // SAFETY: Test code calling backend trait methods marked unsafe
    1524
            unsafe {
    1525
                super::super::scalar::ScalarBackend::swish(&a, &mut scalar_result);
    1526
                Sse2Backend::swish(&a, &mut sse2_result);
    1527
            }
    1528
    1529
            for (s, e) in scalar_result.iter().zip(sse2_result.iter()) {
    1530
                assert!(
    1531
                    (s - e).abs() < 1e-5,
    1532
                    "swish mismatch: scalar={}, sse2={}",
    1533
                    s,
    1534
                    e
    1535
                );
    1536
            }
    1537
        }
    1538
    1539
        #[test]
    1540
        fn test_sse2_sub_matches_scalar() {
    1541
            let a = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0];
    1542
            let b = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0];
    1543
    1544
            let mut scalar_result = [0.0; 7];
    1545
            let mut sse2_result = [0.0; 7];
    1546
    1547
            // SAFETY: Test code calling backend trait methods marked unsafe
    1548
            unsafe {
    1549
                super::super::scalar::ScalarBackend::sub(&a, &b, &mut scalar_result);
    1550
                Sse2Backend::sub(&a, &b, &mut sse2_result);
    1551
            }
    1552
    1553
            assert_eq!(scalar_result, sse2_result);
    1554
        }
    1555
    1556
        #[test]
    1557
        fn test_sse2_div_matches_scalar() {
    1558
            let a = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0];
    1559
            let b = [2.0, 4.0, 5.0, 8.0, 10.0, 12.0, 14.0];
    1560
    1561
            let mut scalar_result = [0.0; 7];
    1562
            let mut sse2_result = [0.0; 7];
    1563
    1564
            // SAFETY: Test code calling backend trait methods marked unsafe
    1565
            unsafe {
    1566
                super::super::scalar::ScalarBackend::div(&a, &b, &mut scalar_result);
    1567
                Sse2Backend::div(&a, &b, &mut sse2_result);
    1568
            }
    1569
    1570
            // Use tolerance-based comparison since rcp+refinement has ~5e-7 relative error
    1571
            for (i, (&s, &sse2)) in scalar_result.iter().zip(sse2_result.iter()).enumerate() {
    1572
                let rel_error = ((s - sse2) / s).abs();
    1573
                assert!(
    1574
                    rel_error < 1e-5,
    1575
                    "Div mismatch at index {}: scalar={}, sse2={}, rel_error={}",
    1576
                    i,
    1577
                    s,
    1578
                    sse2,
    1579
                    rel_error
    1580
                );
    1581
            }
    1582
        }
    1583
    1584
        #[test]
    1585
        fn test_sse2_scale_matches_scalar() {
    1586
            let a = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0];
    1587
            let scalar = 2.5;
    1588
    1589
            let mut scalar_result = [0.0; 7];
    1590
            let mut sse2_result = [0.0; 7];
    1591
    1592
            // SAFETY: Test code calling backend trait methods marked unsafe
    1593
            unsafe {
    1594
                super::super::scalar::ScalarBackend::scale(&a, scalar, &mut scalar_result);
    1595
                Sse2Backend::scale(&a, scalar, &mut sse2_result);
    1596
            }
    1597
    1598
            assert_eq!(scalar_result, sse2_result);
    1599
        }
    1600
    1601
        #[test]
    1602
        fn test_sse2_clamp_matches_scalar() {
    1603
            let a = [1.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0];
    1604
    1605
            let mut scalar_result = [0.0; 7];
    1606
            let mut sse2_result = [0.0; 7];
    1607
    1608
            // SAFETY: Test code calling backend trait methods marked unsafe
    1609
            unsafe {
    1610
                super::super::scalar::ScalarBackend::clamp(&a, 5.0, 20.0, &mut scalar_result);
    1611
                Sse2Backend::clamp(&a, 5.0, 20.0, &mut sse2_result);
    1612
            }
    1613
    1614
            assert_eq!(scalar_result, sse2_result);
    1615
        }
    1616
    1617
        #[test]
    1618
        fn test_sse2_fma_matches_scalar() {
    1619
            let a = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0];
    1620
            let b = [2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0];
    1621
            let c = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0];
    1622
    1623
            let mut scalar_result = [0.0; 7];
    1624
            let mut sse2_result = [0.0; 7];
    1625
    1626
            // SAFETY: Test code calling backend trait methods marked unsafe
    1627
            unsafe {
    1628
                super::super::scalar::ScalarBackend::fma(&a, &b, &c, &mut scalar_result);
    1629
                Sse2Backend::fma(&a, &b, &c, &mut sse2_result);
    1630
            }
    1631
    1632
            assert_eq!(scalar_result, sse2_result);
    1633
        }
    1634
    1635
        #[test]
    1636
        fn test_sse2_lerp_matches_scalar() {
    1637
            let a = [0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0];
    1638
            let b = [100.0, 110.0, 120.0, 130.0, 140.0, 150.0, 160.0];
    1639
    1640
            let mut scalar_result = [0.0; 7];
    1641
            let mut sse2_result = [0.0; 7];
    1642
    1643
            // SAFETY: Test code calling backend trait methods marked unsafe
    1644
            unsafe {
    1645
                super::super::scalar::ScalarBackend::lerp(&a, &b, 0.25, &mut scalar_result);
    1646
                Sse2Backend::lerp(&a, &b, 0.25, &mut sse2_result);
    1647
            }
    1648
    1649
            for (s, e) in scalar_result.iter().zip(sse2_result.iter()) {
    1650
                assert!(
    1651
                    (s - e).abs() < 1e-5,
    1652
                    "lerp mismatch: scalar={}, sse2={}",
    1653
                    s,
    1654
                    e
    1655
                );
    1656
            }
    1657
        }
    1658
    1659
        #[test]
    1660
        fn test_sse2_argmax_matches_scalar() {
    1661
            let a = [1.0, 5.0, 3.0, 10.0, 2.0, 8.0, 4.0];
    1662
    1663
            // SAFETY: Test code calling backend trait methods marked unsafe
    1664
            let scalar_result = unsafe { super::super::scalar::ScalarBackend::argmax(&a) };
    1665
            // SAFETY: CPU feature verified at runtime, slices bounds-checked
    1666
            let sse2_result = unsafe { Sse2Backend::argmax(&a) };
    1667
    1668
            assert_eq!(scalar_result, sse2_result);
    1669
        }
    1670
    1671
        #[test]
    1672
        fn test_sse2_argmin_matches_scalar() {
    1673
            let a = [5.0, 1.0, 3.0, 10.0, 2.0, 8.0, 4.0];
    1674
    1675
            // SAFETY: Test code calling backend trait methods marked unsafe
    1676
            let scalar_result = unsafe { super::super::scalar::ScalarBackend::argmin(&a) };
    1677
            // SAFETY: CPU feature verified at runtime, slices bounds-checked
    1678
            let sse2_result = unsafe { Sse2Backend::argmin(&a) };
    1679
    1680
            assert_eq!(scalar_result, sse2_result);
    1681
        }
    1682
    1683
        #[test]
    1684
        fn test_sse2_sum_kahan_matches_scalar() {
    1685
            let a = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0];
    1686
    1687
            // SAFETY: Test code calling backend trait methods marked unsafe
    1688
            let scalar_result = unsafe { super::super::scalar::ScalarBackend::sum_kahan(&a) };
    1689
            // SAFETY: CPU feature verified at runtime, slices bounds-checked
    1690
            let sse2_result = unsafe { Sse2Backend::sum_kahan(&a) };
    1691
    1692
            assert!((scalar_result - sse2_result).abs() < 1e-5);
    1693
        }
    1694
    1695
        #[test]
    1696
        fn test_sse2_norm_l1_matches_scalar() {
    1697
            let a = [1.0, -2.0, 3.0, -4.0, 5.0, -6.0, 7.0];
    1698
    1699
            // SAFETY: Test code calling backend trait methods marked unsafe
    1700
            let scalar_result = unsafe { super::super::scalar::ScalarBackend::norm_l1(&a) };
    1701
            // SAFETY: CPU feature verified at runtime, slices bounds-checked
    1702
            let sse2_result = unsafe { Sse2Backend::norm_l1(&a) };
    1703
    1704
            assert!((scalar_result - sse2_result).abs() < 1e-5);
    1705
        }
    1706
    1707
        #[test]
    1708
        fn test_sse2_norm_l2_matches_scalar() {
    1709
            let a = [3.0, 4.0, 0.0, 0.0, 5.0, 12.0, 0.0];
    1710
    1711
            // SAFETY: Test code calling backend trait methods marked unsafe
    1712
            let scalar_result = unsafe { super::super::scalar::ScalarBackend::norm_l2(&a) };
    1713
            // SAFETY: CPU feature verified at runtime, slices bounds-checked
    1714
            let sse2_result = unsafe { Sse2Backend::norm_l2(&a) };
    1715
    1716
            assert!((scalar_result - sse2_result).abs() < 1e-5);
    1717
        }
    1718
    1719
        #[test]
    1720
        fn test_sse2_dot_matches_scalar() {
    1721
            let a = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0];
    1722
            let b = [7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0];
    1723
    1724
            // SAFETY: Test code calling backend trait methods marked unsafe
    1725
            let scalar_result = unsafe { super::super::scalar::ScalarBackend::dot(&a, &b) };
    1726
            // SAFETY: CPU feature verified at runtime, slices bounds-checked
    1727
            let sse2_result = unsafe { Sse2Backend::dot(&a, &b) };
    1728
    1729
            assert!((scalar_result - sse2_result).abs() < 1e-5);
    1730
        }
    1731
    1732
        #[test]
    1733
        fn test_sse2_mul_matches_scalar() {
    1734
            let a = [1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5];
    1735
            let b = [2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0];
    1736
    1737
            let mut scalar_result = [0.0; 7];
    1738
            let mut sse2_result = [0.0; 7];
    1739
    1740
            // SAFETY: Test code calling backend trait methods marked unsafe
    1741
            unsafe {
    1742
                super::super::scalar::ScalarBackend::mul(&a, &b, &mut scalar_result);
    1743
                Sse2Backend::mul(&a, &b, &mut sse2_result);
    1744
            }
    1745
    1746
            assert_eq!(scalar_result, sse2_result);
    1747
        }
    1748
    1749
        #[test]
    1750
        fn test_sse2_add_matches_scalar() {
    1751
            let a = [1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5];
    1752
            let b = [8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5];
    1753
    1754
            let mut scalar_result = [0.0; 7];
    1755
            let mut sse2_result = [0.0; 7];
    1756
    1757
            // SAFETY: Test code calling backend trait methods marked unsafe
    1758
            unsafe {
    1759
                super::super::scalar::ScalarBackend::add(&a, &b, &mut scalar_result);
    1760
                Sse2Backend::add(&a, &b, &mut sse2_result);
    1761
            }
    1762
    1763
            assert_eq!(scalar_result, sse2_result);
    1764
        }
    1765
    1766
        #[test]
    1767
        fn test_sse2_sum_matches_scalar() {
    1768
            let a = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0];
    1769
    1770
            // SAFETY: Test code calling backend trait methods marked unsafe
    1771
            let scalar_result = unsafe { super::super::scalar::ScalarBackend::sum(&a) };
    1772
            // SAFETY: CPU feature verified at runtime, slices bounds-checked
    1773
            let sse2_result = unsafe { Sse2Backend::sum(&a) };
    1774
    1775
            assert!((scalar_result - sse2_result).abs() < 1e-5);
    1776
        }
    1777
    1778
        #[test]
    1779
        fn test_sse2_max_matches_scalar() {
    1780
            let a = [1.0, 5.0, 3.0, 7.0, 2.0, 8.0, 4.0];
    1781
    1782
            // SAFETY: Test code calling backend trait methods marked unsafe
    1783
            let scalar_result = unsafe { super::super::scalar::ScalarBackend::max(&a) };
    1784
            // SAFETY: CPU feature verified at runtime, slices bounds-checked
    1785
            let sse2_result = unsafe { Sse2Backend::max(&a) };
    1786
    1787
            assert_eq!(scalar_result, sse2_result);
    1788
        }
    1789
    1790
        #[test]
    1791
        fn test_sse2_min_matches_scalar() {
    1792
            let a = [5.0, 1.0, 3.0, 7.0, 2.0, 8.0, 4.0];
    1793
    1794
            // SAFETY: Test code calling backend trait methods marked unsafe
    1795
            let scalar_result = unsafe { super::super::scalar::ScalarBackend::min(&a) };
    1796
            // SAFETY: CPU feature verified at runtime, slices bounds-checked
    1797
            let sse2_result = unsafe { Sse2Backend::min(&a) };
    1798
    1799
            assert_eq!(scalar_result, sse2_result);
    1800
        }
    1801
    1802
        #[test]
    1803
        fn test_sse2_tanh_matches_scalar() {
    1804
            // Verify SSE2 tanh produces same results as scalar
    1805
            let a = [-10.0, -1.0, 0.0, 1.0, 10.0];
    1806
    1807
            let mut scalar_result = [0.0; 5];
    1808
            let mut sse2_result = [0.0; 5];
    1809
    1810
            // SAFETY: Test code calling backend trait methods marked unsafe
    1811
            unsafe {
    1812
                super::super::scalar::ScalarBackend::tanh(&a, &mut scalar_result);
    1813
                Sse2Backend::tanh(&a, &mut sse2_result);
    1814
            }
    1815
    1816
            for (s, e) in scalar_result.iter().zip(sse2_result.iter()) {
    1817
                assert!(
    1818
                    (s - e).abs() < 1e-5,
    1819
                    "tanh mismatch: scalar={}, sse2={}",
    1820
                    s,
    1821
                    e
    1822
                );
    1823
            }
    1824
        }
    1825
    1826
        #[test]
    1827
        fn test_sse2_norm_linf_matches_scalar() {
    1828
            // Verify SSE2 norm_linf produces same results as scalar
    1829
            let test_cases = vec![
    1830
                vec![],                                       // empty
    1831
                vec![5.0],                                    // single element
    1832
                vec![-3.0, 1.0, -4.0, 1.0, 5.0],              // various values
    1833
                vec![-10.0, 5.0, 3.0, 7.0, -2.0, 8.0, 4.0],   // 7 elements (remainder)
    1834
                vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0], // 8 elements (aligned)
    1835
            ];
    1836
    1837
            for test_vec in test_cases {
    1838
                // SAFETY: Test code calling backend trait methods marked unsafe
    1839
                let scalar_result =
    1840
                    // SAFETY: CPU feature verified at runtime, slices bounds-checked
    1841
                    unsafe { super::super::scalar::ScalarBackend::norm_linf(&test_vec) };
    1842
                // SAFETY: CPU feature verified at runtime, slices bounds-checked
    1843
                let sse2_result = unsafe { Sse2Backend::norm_linf(&test_vec) };
    1844
    1845
                assert!(
    1846
                    (scalar_result - sse2_result).abs() < 1e-5,
    1847
                    "norm_linf mismatch for {:?}: scalar={}, sse2={}",
    1848
                    test_vec,
    1849
                    scalar_result,
    1850
                    sse2_result
    1851
                );
    1852
            }
    1853
        }
    1854
    1855
        #[test]
    1856
        fn test_sse2_abs_matches_scalar() {
    1857
            // Verify SSE2 abs produces same results as scalar
    1858
            let test_cases = vec![
    1859
                vec![],                                             // empty
    1860
                vec![-5.0],                                         // single negative
    1861
                vec![5.0],                                          // single positive
    1862
                vec![-3.0, 1.0, -4.0, 1.5],                         // 4 elements (aligned)
    1863
                vec![-3.0, 1.0, -4.0, 1.5, -9.0, 2.0, -6.0],        // 7 elements (remainder)
    1864
                vec![-1.0, 2.0, -3.0, 4.0, -5.0, 6.0, -7.0, 8.0],   // 8 elements
    1865
                vec![0.0, -0.0, f32::INFINITY, f32::NEG_INFINITY],  // special values
    1866
            ];
    1867
    1868
            for test_vec in test_cases {
    1869
                let mut scalar_result = vec![0.0f32; test_vec.len()];
    1870
                let mut sse2_result = vec![0.0f32; test_vec.len()];
    1871
    1872
                // SAFETY: CPU feature verified at runtime, slices bounds-checked
    1873
                unsafe {
    1874
                    super::super::scalar::ScalarBackend::abs(&test_vec, &mut scalar_result);
    1875
                    Sse2Backend::abs(&test_vec, &mut sse2_result);
    1876
                }
    1877
    1878
                for (i, (&s, &e)) in scalar_result.iter().zip(sse2_result.iter()).enumerate() {
    1879
                    // Handle NaN comparison
    1880
                    if s.is_nan() && e.is_nan() {
    1881
                        continue;
    1882
                    }
    1883
                    assert!(
    1884
                        (s - e).abs() < 1e-5 || (s.is_infinite() && e.is_infinite() && s.signum() == e.signum()),
    1885
                        "abs mismatch at index {} for {:?}: scalar={}, sse2={}",
    1886
                        i,
    1887
                        test_vec,
    1888
                        s,
    1889
                        e
    1890
                    );
    1891
                }
    1892
            }
    1893
        }
    1894
    1895
        #[test]
    1896
        fn test_sse2_tanh_saturation() {
    1897
            // Test tanh saturation at extreme values (scalar remainder path)
    1898
            let extreme_values = vec![-100.0, -50.0, -31.0, 31.0, 50.0, 100.0];
    1899
            let mut result = vec![0.0f32; extreme_values.len()];
    1900
    1901
            // SAFETY: CPU feature verified at runtime
    1902
            unsafe {
    1903
                Sse2Backend::tanh(&extreme_values, &mut result);
    1904
            }
    1905
    1906
            // Values < -30 should saturate to -1.0
    1907
            assert!((result[0] - (-1.0)).abs() < 1e-5, "tanh(-100) should be -1.0");
    1908
            assert!((result[1] - (-1.0)).abs() < 1e-5, "tanh(-50) should be -1.0");
    1909
            assert!((result[2] - (-1.0)).abs() < 1e-5, "tanh(-31) should be -1.0");
    1910
    1911
            // Values > 30 should saturate to 1.0
    1912
            assert!((result[3] - 1.0).abs() < 1e-5, "tanh(31) should be 1.0");
    1913
            assert!((result[4] - 1.0).abs() < 1e-5, "tanh(50) should be 1.0");
    1914
            assert!((result[5] - 1.0).abs() < 1e-5, "tanh(100) should be 1.0");
    1915
        }
    1916
    1917
        #[test]
    1918
        fn test_sse2_gelu_edge_cases() {
    1919
            // Test GELU with edge case values
    1920
            let values = vec![-10.0, -5.0, 0.0, 5.0, 10.0];
    1921
            let mut result = vec![0.0f32; values.len()];
    1922
    1923
            // SAFETY: CPU feature verified at runtime
    1924
            unsafe {
    1925
                Sse2Backend::gelu(&values, &mut result);
    1926
            }
    1927
    1928
            // GELU(-10) ≈ 0 (very negative values)
    1929
            assert!(result[0].abs() < 1e-3, "GELU(-10) should be near 0");
    1930
    1931
            // GELU(0) = 0
    1932
            assert!(result[2].abs() < 1e-5, "GELU(0) should be 0");
    1933
    1934
            // GELU(10) ≈ 10 (very positive values)
    1935
            assert!((result[4] - 10.0).abs() < 0.1, "GELU(10) should be near 10");
    1936
        }
    1937
    1938
        #[test]
    1939
        fn test_sse2_sigmoid_edge_cases() {
    1940
            // Test sigmoid saturation
    1941
            let values = vec![-100.0, -20.0, 0.0, 20.0, 100.0];
    1942
            let mut result = vec![0.0f32; values.len()];
    1943
    1944
            // SAFETY: CPU feature verified at runtime
    1945
            unsafe {
    1946
                Sse2Backend::sigmoid(&values, &mut result);
    1947
            }
    1948
    1949
            // sigmoid(-100) ≈ 0
    1950
            assert!(result[0] < 1e-5, "sigmoid(-100) should be near 0");
    1951
    1952
            // sigmoid(0) = 0.5
    1953
            assert!((result[2] - 0.5).abs() < 1e-5, "sigmoid(0) should be 0.5");
    1954
    1955
            // sigmoid(100) ≈ 1
    1956
            assert!((result[4] - 1.0).abs() < 1e-5, "sigmoid(100) should be near 1");
    1957
        }
    1958
    1959
        #[test]
    1960
        fn test_sse2_exp_edge_cases() {
    1961
            // Test exp with edge values that exercise saturation paths
    1962
            let values = vec![-100.0, -50.0, 0.0, 50.0, 88.0];
    1963
            let mut result = vec![0.0f32; values.len()];
    1964
    1965
            // SAFETY: CPU feature verified at runtime
    1966
            unsafe {
    1967
                Sse2Backend::exp(&values, &mut result);
    1968
            }
    1969
    1970
            // exp(-100) ≈ 0
    1971
            assert!(result[0] < 1e-30, "exp(-100) should be near 0");
    1972
    1973
            // exp(0) = 1
    1974
            assert!((result[2] - 1.0).abs() < 1e-5, "exp(0) should be 1");
    1975
    1976
            // exp(88) is large but finite
    1977
            assert!(result[4].is_finite(), "exp(88) should be finite");
    1978
        }
    1979
    }