## Coverage Report

#### Created: 2026-01-25 15:05

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

    /home/noah/src/realizar/src/quantize/simd.rs

    Line
    Count
    Source
    1
    //! Quantization SIMD Helpers (PMAT-802)
    2
    //!
    3
    //! Extracted from quantize/mod.rs - Shared SIMD utility functions.
    4
    //!
    5
    //! ## Contents
    6
    //! - f16 conversion: `f16_to_f32`, `read_f16`
    7
    //! - Scale extraction: `extract_scale_min`, `extract_scale_min_from_slice`
    8
    //! - Horizontal sum helpers: `hsum_epi32_128`, `hsum_epi32_256`, etc.
    9
    //! - SIMD activations: `softmax_simd`, `fused_swiglu_simd`, `apply_rope_rotation_simd`
    10
    11
    // ============================================================================
    12
    // f16 Conversion (Manual Implementation)
    13
    // ============================================================================
    14
    15
    /// Convert IEEE 754 half-precision (f16) to single-precision (f32)
    16
    ///
    17
    /// Handles normal values, subnormals, infinities, and NaN.
    18
    #[inline]
    19
    111
    pub fn f16_to_f32(h: u16) -> f32 {
    20
    111
        let sign = (h >> 15) & 1;
    21
    111
        let exp = (h >> 10) & 0x1F;
    22
    111
        let mantissa = h & 0x3FF;
    23
    24
    111
        if exp == 0 {
    25
            // Subnormal or zero
    26
    5
            if mantissa == 0 {
    27
                // Zero (preserve sign)
    28
    2
                if sign == 1 {
    29
    1
                    -0.0
    30
                } else {
    31
    1
                    0.0
    32
                }
    33
            } else {
    34
                // Subnormal: (mantissa / 1024) * 2^-14
    35
    3
                let value = (mantissa as f32 / 1024.0) * (2.0_f32).powi(-14);
    36
    3
                if sign == 1 {
    37
    0
                    -value
    38
                } else {
    39
    3
                    value
    40
                }
    41
            }
    42
    106
        } else if exp == 31 {
    43
            // Infinity or NaN
    44
    3
            if mantissa == 0 {
    45
    2
                if sign == 1 {
    46
    1
                    f32::NEG_INFINITY
    47
                } else {
    48
    1
                    f32::INFINITY
    49
                }
    50
            } else {
    51
    1
                f32::NAN
    52
            }
    53
        } else {
    54
            // Normal value: (1 + mantissa/1024) * 2^(exp-15)
    55
    103
            let value = (1.0 + mantissa as f32 / 1024.0) * (2.0_f32).powi(exp as i32 - 15);
    56
    103
            if sign == 1 {
    57
    1
                -value
    58
            } else {
    59
    102
                value
    60
            }
    61
        }
    62
    111
    }
    63
    64
    /// Helper: Read f16 from bytes and convert to f32
    65
    #[inline]
    66
    6
    pub fn read_f16(bytes: &[u8]) -> f32 {
    67
    6
        let bits = u16::from_le_bytes([bytes[0], bytes[1]]);
    68
    6
        half::f16::from_bits(bits).to_f32()
    69
    6
    }
    70
    71
    // ============================================================================
    72
    // Scale Extraction for K-Quantization
    73
    // ============================================================================
    74
    75
    /// Extract 6-bit scale and min values from packed scales array
    76
    ///
    77
    /// PAR-001 FIX: Matches llama.cpp's get_scale_min_k4 packing scheme:
    78
    /// - Blocks 0-3: scale = q[j] & 63, min = q[j+4] & 63
    79
    /// - Blocks 4-7: scale = (q[j+4] & 0xF) | ((q[j-4] >> 6) << 4)
    80
    ///   min = (q[j+4] >> 4) | ((q[j] >> 6) << 4)
    81
    #[inline]
    82
    10.0M
    pub fn extract_scale_min(scales: &[u8; 12], block_idx: usize) -> (f32, f32) {
    83
    10.0M
        let j = block_idx;
    84
    10.0M
        let (scale_bits, min_bits) = if j < 4 {
    85
            // First 4 blocks: simple layout
    86
    5.04M
            let d = scales[j] & 63;
    87
    5.04M
            let m = scales[j + 4] & 63;
    88
    5.04M
            (d, m)
    89
        } else {
    90
            // Last 4 blocks: packed layout using high bits from first 4 bytes
    91
    5.04M
            let d = (scales[j + 4] & 0x0F) | ((scales[j - 4] >> 6) << 4);
    92
    5.04M
            let m = (scales[j + 4] >> 4) | ((scales[j] >> 6) << 4);
    93
    5.04M
            (d, m)
    94
        };
    95
    96
        // Return raw 6-bit values as floats
    97
        // The GGUF header's d/dmin values already include the /63 normalization
    98
    10.0M
        let scale = f32::from(scale_bits);
    99
    10.0M
        let min = f32::from(min_bits);
    100
    101
    10.0M
        (scale, min)
    102
    10.0M
    }
    103
    104
    /// Extract scale and min from packed 6-bit scales (helper for InterleavedQ4K)
    105
    1
    pub fn extract_scale_min_from_slice(scales: &[u8], idx: usize) -> (f32, f32) {
    106
        // Same logic as extract_scale_min but works with slice
    107
    1
        let scale_idx = idx / 2;
    108
    1
        let min_idx = idx / 2 + 4;
    109
    110
    1
        let (scale_raw, min_raw) = if idx.is_multiple_of(2) {
    111
    1
            (scales[scale_idx] & 0x3F, scales[min_idx] & 0x3F)
    112
        } else {
    113
    0
            (
    114
    0
                (scales[scale_idx] >> 6) | ((scales[scale_idx + 2] & 0x0F) << 2),
    115
    0
                (scales[min_idx] >> 6) | ((scales[min_idx + 2] & 0x0F) << 2),
    116
    0
            )
    117
        };
    118
    119
    1
        (scale_raw as f32, min_raw as f32)
    120
    1
    }
    121
    122
    // ============================================================================
    123
    // x86_64 SIMD Horizontal Sum Helpers
    124
    // ============================================================================
    125
    126
    /// Fast horizontal sum of 4 i32 in __m128i
    127
    ///
    128
    /// # Safety
    129
    /// Requires AVX2 support. Caller must verify CPU feature availability.
    130
    #[cfg(target_arch = "x86_64")]
    131
    #[target_feature(enable = "avx2")]
    132
    #[inline]
    133
    0
    pub unsafe fn hsum_epi32_128(v: std::arch::x86_64::__m128i) -> i32 {
    134
        use std::arch::x86_64::{_mm_cvtsi128_si32, _mm_hadd_epi32};
    135
    0
        let sum64 = _mm_hadd_epi32(v, v);
    136
    0
        let sum32 = _mm_hadd_epi32(sum64, sum64);
    137
    0
        _mm_cvtsi128_si32(sum32)
    138
    0
    }
    139
    140
    /// Fast horizontal sum of 8 i32 in __m256i
    141
    ///
    142
    /// # Safety
    143
    /// Requires AVX2 support. Caller must verify CPU feature availability.
    144
    #[cfg(target_arch = "x86_64")]
    145
    #[target_feature(enable = "avx2")]
    146
    #[inline]
    147
    0
    pub unsafe fn hsum_epi32_256(v: std::arch::x86_64::__m256i) -> i32 {
    148
        use std::arch::x86_64::{_mm256_castsi256_si128, _mm256_extracti128_si256, _mm_add_epi32};
    149
        // SAFETY: Unsafe operation with validated invariants
    150
        unsafe {
    151
    0
            let lo = _mm256_castsi256_si128(v);
    152
    0
            let hi = _mm256_extracti128_si256(v, 1);
    153
    0
            hsum_epi32_128(_mm_add_epi32(lo, hi))
    154
        }
    155
    0
    }
    156
    157
    /// Helper: horizontal sum of 8 int32 values in a 256-bit register
    158
    ///
    159
    /// # Safety
    160
    /// Requires AVX2 support. Caller must verify CPU feature availability.
    161
    #[cfg(target_arch = "x86_64")]
    162
    #[target_feature(enable = "avx2")]
    163
    #[inline]
    164
    0
    pub unsafe fn horizontal_sum_epi32_256(v: std::arch::x86_64::__m256i) -> i32 {
    165
        use std::arch::x86_64::{
    166
            _mm256_castsi256_si128, _mm256_extracti128_si256, _mm_add_epi32, _mm_cvtsi128_si32,
    167
            _mm_hadd_epi32,
    168
        };
    169
    170
        // Add high 128 bits to low 128 bits
    171
    0
        let hi = _mm256_extracti128_si256(v, 1);
    172
    0
        let lo = _mm256_castsi256_si128(v);
    173
    0
        let sum128 = _mm_add_epi32(lo, hi);
    174
    175
        // Horizontal add within 128 bits
    176
    0
        let sum64 = _mm_hadd_epi32(sum128, sum128);
    177
    0
        let sum32 = _mm_hadd_epi32(sum64, sum64);
    178
    179
    0
        _mm_cvtsi128_si32(sum32)
    180
    0
    }
    181
    182
    /// Helper: horizontal sum of 16 int16 values in a 256-bit register
    183
    ///
    184
    /// # Safety
    185
    /// Requires AVX2 support. Caller must verify CPU feature availability.
    186
    #[cfg(target_arch = "x86_64")]
    187
    #[target_feature(enable = "avx2")]
    188
    #[inline]
    189
    #[allow(unsafe_op_in_unsafe_fn)]
    190
    0
    pub unsafe fn horizontal_sum_epi16_256(v: std::arch::x86_64::__m256i) -> i32 {
    191
        use std::arch::x86_64::{_mm256_madd_epi16, _mm256_set1_epi16};
    192
    193
        // Use madd to sum pairs of i16 to i32
    194
    0
        let ones = _mm256_set1_epi16(1);
    195
    0
        let sum_i32 = _mm256_madd_epi16(v, ones);
    196
    197
        // Now sum the 8 i32 values
    198
    0
        horizontal_sum_epi32_256(sum_i32)
    199
    0
    }
    200
    201
    /// Helper: Horizontal sum of 8 i32 values to single i32
    202
    ///
    203
    /// # Safety
    204
    /// Requires AVX2 support. Caller must verify CPU feature availability.
    205
    #[cfg(target_arch = "x86_64")]
    206
    #[target_feature(enable = "avx2")]
    207
    #[inline]
    208
    0
    pub unsafe fn hsum_epi32(v: std::arch::x86_64::__m256i) -> i32 {
    209
        #[allow(clippy::wildcard_imports)]
    210
        use std::arch::x86_64::*;
    211
    212
        // All intrinsics are unsafe and we're in an unsafe fn with target_feature
    213
    0
        let sum128 = _mm_add_epi32(_mm256_castsi256_si128(v), _mm256_extracti128_si256(v, 1));
    214
    0
        let sum64 = _mm_add_epi32(sum128, _mm_shuffle_epi32(sum128, 0b10_11_00_01));
    215
    0
        let sum32 = _mm_add_epi32(sum64, _mm_shuffle_epi32(sum64, 0b00_00_10_10));
    216
    0
        _mm_cvtsi128_si32(sum32)
    217
    0
    }
    218
    219
    // ============================================================================
    220
    // SIMD Activation Functions
    221
    // ============================================================================
    222
    223
    /// SIMD-accelerated in-place softmax
    224
    ///
    225
    /// Uses AVX2 for parallel max-finding and normalization.
    226
    /// Falls back to scalar on non-x86_64 or when AVX2 is unavailable.
    227
    ///
    228
    /// # Arguments
    229
    /// * `x` - Input/output slice, modified in-place to contain softmax values
    230
    216
    pub fn softmax_simd(x: &mut [f32]) {
    231
    216
        if x.is_empty() {
    232
    1
            return;
    233
    215
        }
    234
    235
        #[cfg(target_arch = "x86_64")]
    236
        {
    237
    215
            if is_x86_feature_detected!("avx2") && x.len() >= 8 {
    238
                // SAFETY: AVX2 is available and slice is large enough
    239
    192
                unsafe {
    240
    192
                    softmax_simd_avx2(x);
    241
    192
                }
    242
    192
                return;
    243
    23
            }
    244
        }
    245
    246
        // Scalar fallback
    247
    23
        softmax_scalar(x);
    248
    216
    }
    249
    250
    /// Scalar softmax implementation
    251
    23
    fn softmax_scalar(x: &mut [f32]) {
    252
        // Find max for numerical stability
    253
    23
        let max_val = x.iter().copied().fold(f32::NEG_INFINITY, f32::max);
    254
    255
        // Compute exp(x - max) and sum
    256
    23
        let mut sum = 0.0f32;
    257
    98
        for val in x23.iter_mut23() {
    258
    98
            *val = (*val - max_val).exp();
    259
    98
            sum += *val;
    260
    98
        }
    261
    262
        // Normalize
    263
    23
        let inv_sum = 1.0 / sum;
    264
    98
        for val in x23.iter_mut23() {
    265
    98
            *val *= inv_sum;
    266
    98
        }
    267
    23
    }
    268
    269
    /// AVX2 softmax implementation
    270
    #[cfg(target_arch = "x86_64")]
    271
    #[target_feature(enable = "avx2")]
    272
    #[allow(unsafe_op_in_unsafe_fn)]
    273
    192
    unsafe fn softmax_simd_avx2(x: &mut [f32]) {
    274
        #[allow(clippy::wildcard_imports)]
    275
        use std::arch::x86_64::*;
    276
    277
    192
        let len = x.len();
    278
    279
        // Phase 1: Find max using SIMD
    280
    192
        let mut max_vec = _mm256_set1_ps(f32::NEG_INFINITY);
    281
    192
        let chunks = len / 8;
    282
    283
    762
        for i in 0..chunks192 {
    284
    762
            let v = _mm256_loadu_ps(x.as_ptr().add(i * 8));
    285
    762
            max_vec = _mm256_max_ps(max_vec, v);
    286
    762
        }
    287
    288
        // Horizontal max
    289
    192
        let max_128 = _mm_max_ps(_mm256_castps256_ps128(max_vec), _mm256_extractf128_ps(max_vec, 1));
    290
    192
        let max_64 = _mm_max_ps(max_128, _mm_movehl_ps(max_128, max_128));
    291
    192
        let max_32 = _mm_max_ss(max_64, _mm_shuffle_ps(max_64, max_64, 1));
    292
    192
        let mut max_val = _mm_cvtss_f32(max_32);
    293
    294
        // Check remaining elements
    295
    628
        for i in (chunks * 8)192..len192 {
    296
    628
            max_val = max_val.max(x[i]);
    297
    628
        }
    298
    299
        // Phase 2: Compute exp(x - max) and sum using SIMD
    300
    192
        let max_vec = _mm256_set1_ps(max_val);
    301
    192
        let mut sum_vec = _mm256_setzero_ps();
    302
    303
    762
        for i in 0..chunks192 {
    304
    762
            let v = _mm256_loadu_ps(x.as_ptr().add(i * 8));
    305
    762
            let diff = _mm256_sub_ps(v, max_vec);
    306
    307
            // Use polynomial approximation for exp (faster but less accurate)
    308
            // For production, could use vectorized exp from a math library
    309
            // Here we fall back to scalar exp for accuracy
    310
    762
            let mut exp_vals = [0.0f32; 8];
    311
    762
            let diff_arr: [f32; 8] = std::mem::transmute(diff);
    312
    6.09k
            for (j, &d) in diff_arr762.iter762().enumerate762() {
    313
    6.09k
                exp_vals[j] = d.exp();
    314
    6.09k
            }
    315
    762
            let exp_vec = _mm256_loadu_ps(exp_vals.as_ptr());
    316
    317
    762
            _mm256_storeu_ps(x.as_mut_ptr().add(i * 8), exp_vec);
    318
    762
            sum_vec = _mm256_add_ps(sum_vec, exp_vec);
    319
        }
    320
    321
        // Handle remaining elements
    322
    192
        let mut sum_scalar = 0.0f32;
    323
    628
        for i in (chunks * 8)192..len192 {
    324
    628
            x[i] = (x[i] - max_val).exp();
    325
    628
            sum_scalar += x[i];
    326
    628
        }
    327
    328
        // Horizontal sum
    329
    192
        let sum_128 = _mm_add_ps(_mm256_castps256_ps128(sum_vec), _mm256_extractf128_ps(sum_vec, 1));
    330
    192
        let sum_64 = _mm_add_ps(sum_128, _mm_movehl_ps(sum_128, sum_128));
    331
    192
        let sum_32 = _mm_add_ss(sum_64, _mm_shuffle_ps(sum_64, sum_64, 1));
    332
    192
        let sum = _mm_cvtss_f32(sum_32) + sum_scalar;
    333
    334
        // Phase 3: Normalize
    335
    192
        let inv_sum = 1.0 / sum;
    336
    192
        let inv_sum_vec = _mm256_set1_ps(inv_sum);
    337
    338
    762
        for i in 0..chunks192 {
    339
    762
            let v = _mm256_loadu_ps(x.as_ptr().add(i * 8));
    340
    762
            let normalized = _mm256_mul_ps(v, inv_sum_vec);
    341
    762
            _mm256_storeu_ps(x.as_mut_ptr().add(i * 8), normalized);
    342
    762
        }
    343
    344
    628
        for i in (chunks * 8)192..len192 {
    345
    628
            x[i] *= inv_sum;
    346
    628
        }
    347
    192
    }
    348
    349
    /// Fused SiLU (Swish) with gating: gate = gate * sigmoid(gate) * up
    350
    ///
    351
    /// Modifies `gate` in-place: gate[i] = silu(gate[i]) * up[i]
    352
    ///
    353
    /// # Arguments
    354
    /// * `gate` - Gate values, modified in-place
    355
    /// * `up` - Up-projection values (must be same length as gate)
    356
    ///
    357
    /// # Panics
    358
    /// Panics if `gate.len() != up.len()`
    359
    33
    pub fn fused_swiglu_simd(gate: &mut [f32], up: &[f32]) {
    360
    33
        assert_eq!(
    361
    33
            gate.len(),
    362
    33
            up.len(),
    363
    1
            "fused_swiglu_simd: gate and up must have same length"
    364
        );
    365
    366
        #[cfg(target_arch = "x86_64")]
    367
        {
    368
    32
            if is_x86_feature_detected!("avx2") && is_x86_feature_detected!("fma") && gate.len() >= 8 {
    369
                // SAFETY: AVX2+FMA available, slices are same length
    370
    28
                unsafe {
    371
    28
                    fused_swiglu_simd_avx2(gate, up);
    372
    28
                }
    373
    28
                return;
    374
    4
            }
    375
        }
    376
    377
        // Scalar fallback
    378
    4
        fused_swiglu_scalar(gate, up);
    379
    32
    }
    380
    381
    /// Scalar SwiGLU implementation
    382
    4
    fn fused_swiglu_scalar(gate: &mut [f32], up: &[f32]) {
    383
    16
        for (g, &u) in gate4.iter_mut4().zip4(up4.iter4()) {
    384
    16
            // SiLU: x * sigmoid(x) = x / (1 + exp(-x))
    385
    16
            let silu = *g / (1.0 + (-*g).exp());
    386
    16
            *g = silu * u;
    387
    16
        }
    388
    4
    }
    389
    390
    /// AVX2 SwiGLU implementation
    391
    #[cfg(target_arch = "x86_64")]
    392
    #[target_feature(enable = "avx2", enable = "fma")]
    393
    #[allow(unsafe_op_in_unsafe_fn)]
    394
    28
    unsafe fn fused_swiglu_simd_avx2(gate: &mut [f32], up: &[f32]) {
    395
        #[allow(clippy::wildcard_imports)]
    396
        use std::arch::x86_64::*;
    397
    398
    28
        let len = gate.len();
    399
    28
        let chunks = len / 8;
    400
    28
        let _one = _mm256_set1_ps(1.0);
    401
    402
    116
        for i in 0..chunks28 {
    403
    116
            let g = _mm256_loadu_ps(gate.as_ptr().add(i * 8));
    404
    116
            let u = _mm256_loadu_ps(up.as_ptr().add(i * 8));
    405
    406
            // Compute sigmoid using scalar (for accuracy)
    407
            // A fully vectorized version would use polynomial approximation
    408
    116
            let mut sigmoid_vals = [0.0f32; 8];
    409
    116
            let g_arr: [f32; 8] = std::mem::transmute(g);
    410
    928
            for (j, &gv) in g_arr116.iter116().enumerate116() {
    411
    928
                sigmoid_vals[j] = 1.0 / (1.0 + (-gv).exp());
    412
    928
            }
    413
    116
            let sigmoid = _mm256_loadu_ps(sigmoid_vals.as_ptr());
    414
    415
            // SiLU = g * sigmoid(g)
    416
    116
            let silu = _mm256_mul_ps(g, sigmoid);
    417
    418
            // Result = silu * up
    419
    116
            let result = _mm256_mul_ps(silu, u);
    420
    421
    116
            _mm256_storeu_ps(gate.as_mut_ptr().add(i * 8), result);
    422
        }
    423
    424
        // Handle remaining elements
    425
    81
        for i in (chunks * 8)28..len28 {
    426
    81
            let g = gate[i];
    427
    81
            let sigmoid = 1.0 / (1.0 + (-g).exp());
    428
    81
            gate[i] = g * sigmoid * up[i];
    429
    81
        }
    430
    28
    }
    431
    432
    /// Apply RoPE (Rotary Position Embeddings) rotation with SIMD
    433
    ///
    434
    /// Applies rotary position encoding to query/key vectors.
    435
    ///
    436
    /// # Arguments
    437
    /// * `x` - Input tensor, modified in-place (shape: [seq_len, num_heads, head_dim])
    438
    /// * `freqs_cis` - Complex exponentials (cos, sin) for each position
    439
    /// * `head_dim` - Dimension per head (must be even)
    440
    32
    pub fn apply_rope_rotation_simd(x: &mut [f32], freqs_cos: &[f32], freqs_sin: &[f32], head_dim: usize) {
    441
    32
        debug_assert!(head_dim.is_multiple_of(2), "head_dim must be even"0);
    442
    32
        debug_assert_eq!(freqs_cos.len(), freqs_sin.len());
    443
    444
    32
        let half_dim = head_dim / 2;
    445
    446
        #[cfg(target_arch = "x86_64")]
    447
        {
    448
    32
            if is_x86_feature_detected!("avx2") && is_x86_feature_detected!("fma") && half_dim >= 8 {
    449
                // SAFETY: AVX2+FMA available
    450
    1
                unsafe {
    451
    1
                    apply_rope_rotation_avx2(x, freqs_cos, freqs_sin, head_dim);
    452
    1
                }
    453
    1
                return;
    454
    31
            }
    455
        }
    456
    457
        // Scalar fallback
    458
    31
        apply_rope_rotation_scalar(x, freqs_cos, freqs_sin, head_dim);
    459
    32
    }
    460
    461
    /// Scalar RoPE implementation
    462
    31
    fn apply_rope_rotation_scalar(x: &mut [f32], freqs_cos: &[f32], freqs_sin: &[f32], head_dim: usize) {
    463
    31
        let half_dim = head_dim / 2;
    464
    465
        // Process pairs of values
    466
    119
        for i in 0..half_dim31 {
    467
    119
            if i >= freqs_cos.len() {
    468
    1
                break;
    469
    118
            }
    470
    471
    118
            let cos = freqs_cos[i];
    472
    118
            let sin = freqs_sin[i];
    473
    474
    118
            let x0 = x[i];
    475
    118
            let x1 = x[i + half_dim];
    476
    477
            // Apply rotation: [cos -sin; sin cos] @ [x0; x1]
    478
    118
            x[i] = x0 * cos - x1 * sin;
    479
    118
            x[i + half_dim] = x0 * sin + x1 * cos;
    480
        }
    481
    31
    }
    482
    483
    /// AVX2 RoPE implementation
    484
    #[cfg(target_arch = "x86_64")]
    485
    #[target_feature(enable = "avx2", enable = "fma")]
    486
    #[allow(unsafe_op_in_unsafe_fn)]
    487
    1
    unsafe fn apply_rope_rotation_avx2(
    488
    1
        x: &mut [f32],
    489
    1
        freqs_cos: &[f32],
    490
    1
        freqs_sin: &[f32],
    491
    1
        head_dim: usize,
    492
    1
    ) {
    493
        #[allow(clippy::wildcard_imports)]
    494
        use std::arch::x86_64::*;
    495
    496
    1
        let half_dim = head_dim / 2;
    497
    1
        let chunks = half_dim / 8;
    498
    499
    1
        for i in 0..chunks {
    500
    1
            let offset = i * 8;
    501
    1
    502
    1
            // Load x0 and x1 (first and second half of dimension)
    503
    1
            let x0 = _mm256_loadu_ps(x.as_ptr().add(offset));
    504
    1
            let x1 = _mm256_loadu_ps(x.as_ptr().add(offset + half_dim));
    505
    1
    506
    1
            // Load cos and sin values
    507
    1
            let cos = _mm256_loadu_ps(freqs_cos.as_ptr().add(offset));
    508
    1
            let sin = _mm256_loadu_ps(freqs_sin.as_ptr().add(offset));
    509
    1
    510
    1
            // Compute: x0' = x0 * cos - x1 * sin
    511
    1
            let x0_new = _mm256_fmsub_ps(x0, cos, _mm256_mul_ps(x1, sin));
    512
    1
    513
    1
            // Compute: x1' = x0 * sin + x1 * cos
    514
    1
            let x1_new = _mm256_fmadd_ps(x0, sin, _mm256_mul_ps(x1, cos));
    515
    1
    516
    1
            // Store results
    517
    1
            _mm256_storeu_ps(x.as_mut_ptr().add(offset), x0_new);
    518
    1
            _mm256_storeu_ps(x.as_mut_ptr().add(offset + half_dim), x1_new);
    519
    1
        }
    520
    521
        // Handle remaining elements
    522
    1
        let remaining_start = chunks * 8;
    523
    1
        for i0 in remaining_start..half_dim {
    524
    0
            if i >= freqs_cos.len() {
    525
    0
                break;
    526
    0
            }
    527
    528
    0
            let cos = freqs_cos[i];
    529
    0
            let sin = freqs_sin[i];
    530
    531
    0
            let x0 = x[i];
    532
    0
            let x1 = x[i + half_dim];
    533
    534
    0
            x[i] = x0 * cos - x1 * sin;
    535
    0
            x[i + half_dim] = x0 * sin + x1 * cos;
    536
        }
    537
    1
    }