// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// simd_backend.hpp
//
// SIMD-optimized backend implementation
//
// This file provides optimized implementations of linear algebra operations
// using SIMD instruction sets. It accelerates vector and matrix operations
// with instruction sets such as AVX, AVX2, AVX512, SSE2, and SSE4.1.
//
// Main features:
// - Support for different SIMD instruction sets
// - Optimization of basic vector operations
// - SIMD instruction set detection and selection

#ifndef SANGI_SIMD_BACKEND_HPP
#define SANGI_SIMD_BACKEND_HPP

#include "math/core/common.hpp"
#include "math/core/traits.hpp"

#include <cstddef>
#include <type_traits>
#include <algorithm>
#include <cmath>

// Detect the SIMD instruction set
#if defined(__AVX512F__)
#define SANGI_HAS_AVX512 1
#endif

#if defined(__AVX2__)
#define SANGI_HAS_AVX2 1
#endif

#if defined(__AVX__)
#define SANGI_HAS_AVX 1
#endif

#if defined(__SSE4_1__)
#define SANGI_HAS_SSE4_1 1
#endif

#if defined(__SSE2__) || defined(_M_X64) || defined(_M_AMD64) || (defined(_M_IX86_FP) && _M_IX86_FP >= 2)
#define SANGI_HAS_SSE2 1
#endif

// NEON detection for ARM platforms
#if defined(__ARM_NEON) || defined(__ARM_NEON__)
#define SANGI_HAS_NEON 1
#endif

// Includes for SIMD instruction sets
#if defined(SANGI_HAS_AVX512)
#include <immintrin.h>  // Covers AVX512 and below
#elif defined(SANGI_HAS_AVX2)
#include <immintrin.h>  // For AVX2
#elif defined(SANGI_HAS_AVX)
#include <immintrin.h>  // For AVX
#elif defined(SANGI_HAS_SSE4_1)
#include <smmintrin.h>  // For SSE4.1
#elif defined(SANGI_HAS_SSE2)
#include <emmintrin.h>  // For SSE2
#endif

#if defined(SANGI_HAS_NEON)
#include <arm_neon.h>   // For ARM NEON
#endif

namespace sangi {
namespace computation {
namespace simd {

// Enumeration of SIMD instruction sets
enum class SimdInstructionSet {
    None,
    SSE2,
    SSE4_1,
    AVX,
    AVX2,
    AVX512,
    NEON
};

// Detect the available SIMD instruction set
constexpr SimdInstructionSet detect_simd_instruction_set() {
#if defined(SANGI_HAS_AVX512)
    return SimdInstructionSet::AVX512;
#elif defined(SANGI_HAS_AVX2)
    return SimdInstructionSet::AVX2;
#elif defined(SANGI_HAS_AVX)
    return SimdInstructionSet::AVX;
#elif defined(SANGI_HAS_SSE4_1)
    return SimdInstructionSet::SSE4_1;
#elif defined(SANGI_HAS_SSE2)
    return SimdInstructionSet::SSE2;
#elif defined(SANGI_HAS_NEON)
    return SimdInstructionSet::NEON;
#else
    return SimdInstructionSet::None;
#endif
}

// Available SIMD instruction set
constexpr SimdInstructionSet available_simd_level = detect_simd_instruction_set();

//-----------------------------------------------------------------------------
// SIMD-optimization helper functions
//-----------------------------------------------------------------------------

// SIMD-optimized dot product (float)
inline float dot_product_simd_float(const float* a, const float* b, std::size_t size) {
    float result = 0.0f;
    
    // Select the optimal version based on the SIMD instruction set
#if defined(SANGI_HAS_AVX512)
    // AVX512 implementation
    constexpr std::size_t simd_size = 16; // 512 bits = 16 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m512 sum = _mm512_setzero_ps();
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m512 va = _mm512_loadu_ps(a + i);
        __m512 vb = _mm512_loadu_ps(b + i);
        sum = _mm512_fmadd_ps(va, vb, sum);
    }
    
    // Compute the sum via horizontal addition
    result = _mm512_reduce_add_ps(sum);
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#elif defined(SANGI_HAS_AVX2)
    // AVX2 implementation
    constexpr std::size_t simd_size = 8; // 256 bits = 8 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m256 sum = _mm256_setzero_ps();
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m256 va = _mm256_loadu_ps(a + i);
        __m256 vb = _mm256_loadu_ps(b + i);
        // Use the FMA instruction: sum += va * vb
        sum = _mm256_fmadd_ps(va, vb, sum);
    }
    
    // Compute the sum via horizontal addition
    __m128 sum_hi = _mm256_extractf128_ps(sum, 1);
    __m128 sum_lo = _mm256_castps256_ps128(sum);
    __m128 sum_4 = _mm_add_ps(sum_hi, sum_lo);
    __m128 sum_2 = _mm_add_ps(sum_4, _mm_movehl_ps(sum_4, sum_4));
    __m128 sum_1 = _mm_add_ss(sum_2, _mm_shuffle_ps(sum_2, sum_2, 1));
    
    result = _mm_cvtss_f32(sum_1);
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#elif defined(SANGI_HAS_AVX)
    // AVX implementation
    constexpr std::size_t simd_size = 8; // 256 bits = 8 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m256 sum = _mm256_setzero_ps();
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m256 va = _mm256_loadu_ps(a + i);
        __m256 vb = _mm256_loadu_ps(b + i);
        __m256 mul = _mm256_mul_ps(va, vb);
        sum = _mm256_add_ps(sum, mul);
    }
    
    // Compute the sum via horizontal addition
    __m128 sum_hi = _mm256_extractf128_ps(sum, 1);
    __m128 sum_lo = _mm256_castps256_ps128(sum);
    __m128 sum_4 = _mm_add_ps(sum_hi, sum_lo);
    __m128 sum_2 = _mm_add_ps(sum_4, _mm_movehl_ps(sum_4, sum_4));
    __m128 sum_1 = _mm_add_ss(sum_2, _mm_shuffle_ps(sum_2, sum_2, 1));
    
    result = _mm_cvtss_f32(sum_1);
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#elif defined(SANGI_HAS_SSE4_1) || defined(SANGI_HAS_SSE2)
    // SSE2/SSE4.1 implementation
    constexpr std::size_t simd_size = 4; // 128 bits = 4 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m128 sum = _mm_setzero_ps();
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m128 va = _mm_loadu_ps(a + i);
        __m128 vb = _mm_loadu_ps(b + i);
        __m128 mul = _mm_mul_ps(va, vb);
        sum = _mm_add_ps(sum, mul);
    }
    
    // Compute the sum via horizontal addition
    __m128 sum_2 = _mm_add_ps(sum, _mm_movehl_ps(sum, sum));
    __m128 sum_1 = _mm_add_ss(sum_2, _mm_shuffle_ps(sum_2, sum_2, 1));
    
    result = _mm_cvtss_f32(sum_1);
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#elif defined(SANGI_HAS_NEON)
    // ARM NEON implementation
    constexpr std::size_t simd_size = 4; // 128 bits = 4 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    float32x4_t sum = vdupq_n_f32(0.0f);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        float32x4_t va = vld1q_f32(a + i);
        float32x4_t vb = vld1q_f32(b + i);
        // Multiply-accumulate: sum += va * vb
        sum = vmlaq_f32(sum, va, vb);
    }
    
    // Compute the sum via horizontal addition
    float32x2_t sum_2 = vadd_f32(vget_low_f32(sum), vget_high_f32(sum));
    float32x2_t sum_1 = vpadd_f32(sum_2, sum_2);
    
    result = vget_lane_f32(sum_1, 0);
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#else
    // Standard implementation without SIMD
    for (std::size_t i = 0; i < size; ++i) {
        result += a[i] * b[i];
    }
#endif
    
    return result;
}

// SIMD-optimized dot product (double)
inline double dot_product_simd_double(const double* a, const double* b, std::size_t size) {
    double result = 0.0;
    
    // Select the optimal version based on the SIMD instruction set
#if defined(SANGI_HAS_AVX512)
    // AVX512 implementation
    constexpr std::size_t simd_size = 8; // 512 bits = 8 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m512d sum = _mm512_setzero_pd();
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m512d va = _mm512_loadu_pd(a + i);
        __m512d vb = _mm512_loadu_pd(b + i);
        sum = _mm512_fmadd_pd(va, vb, sum);
    }
    
    // Compute the sum via horizontal addition
    result = _mm512_reduce_add_pd(sum);
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#elif defined(SANGI_HAS_AVX2)
    // AVX2 implementation
    constexpr std::size_t simd_size = 4; // 256 bits = 4 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m256d sum = _mm256_setzero_pd();
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m256d va = _mm256_loadu_pd(a + i);
        __m256d vb = _mm256_loadu_pd(b + i);
        // Use the FMA instruction: sum += va * vb
        sum = _mm256_fmadd_pd(va, vb, sum);
    }
    
    // Compute the sum via horizontal addition
    __m128d sum_hi = _mm256_extractf128_pd(sum, 1);
    __m128d sum_lo = _mm256_castpd256_pd128(sum);
    __m128d sum_2 = _mm_add_pd(sum_hi, sum_lo);
    __m128d sum_1 = _mm_add_sd(sum_2, _mm_unpackhi_pd(sum_2, sum_2));
    
    result = _mm_cvtsd_f64(sum_1);
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#elif defined(SANGI_HAS_AVX)
    // AVX implementation
    constexpr std::size_t simd_size = 4; // 256 bits = 4 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m256d sum = _mm256_setzero_pd();
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m256d va = _mm256_loadu_pd(a + i);
        __m256d vb = _mm256_loadu_pd(b + i);
        __m256d mul = _mm256_mul_pd(va, vb);
        sum = _mm256_add_pd(sum, mul);
    }
    
    // Compute the sum via horizontal addition
    __m128d sum_hi = _mm256_extractf128_pd(sum, 1);
    __m128d sum_lo = _mm256_castpd256_pd128(sum);
    __m128d sum_2 = _mm_add_pd(sum_hi, sum_lo);
    __m128d sum_1 = _mm_add_sd(sum_2, _mm_unpackhi_pd(sum_2, sum_2));
    
    result = _mm_cvtsd_f64(sum_1);
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#elif defined(SANGI_HAS_SSE4_1) || defined(SANGI_HAS_SSE2)
    // SSE2/SSE4.1 implementation
    constexpr std::size_t simd_size = 2; // 128 bits = 2 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m128d sum = _mm_setzero_pd();
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m128d va = _mm_loadu_pd(a + i);
        __m128d vb = _mm_loadu_pd(b + i);
        __m128d mul = _mm_mul_pd(va, vb);
        sum = _mm_add_pd(sum, mul);
    }
    
    // Compute the sum via horizontal addition
    __m128d sum_1 = _mm_add_sd(sum, _mm_unpackhi_pd(sum, sum));
    
    result = _mm_cvtsd_f64(sum_1);
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#elif defined(SANGI_HAS_NEON)
    // ARM NEON implementation (note: double-precision NEON instructions are supported on ARMv8 and later)
    constexpr std::size_t simd_size = 2; // 128 bits = 2 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    float64x2_t sum = vdupq_n_f64(0.0);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        float64x2_t va = vld1q_f64(a + i);
        float64x2_t vb = vld1q_f64(b + i);
        // Multiply-accumulate: sum += va * vb
        sum = vmlaq_f64(sum, va, vb);
    }
    
    // Compute the sum via horizontal addition
    result = vgetq_lane_f64(sum, 0) + vgetq_lane_f64(sum, 1);
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#else
    // Standard implementation without SIMD
    for (std::size_t i = 0; i < size; ++i) {
        result += a[i] * b[i];
    }
#endif
    
    return result;
}

// SIMD-optimized dot product (selects the appropriate implementation by type)
template<typename T>
inline T dot_product_simd(const T* a, const T* b, std::size_t size) {
    if constexpr (std::is_same_v<T, float>) {
        return dot_product_simd_float(a, b, size);
    } 
    else if constexpr (std::is_same_v<T, double>) {
        return dot_product_simd_double(a, b, size);
    }
    else {
        // Other types are computed without SIMD optimization
        T result = numeric_traits<T>::zero();
        for (std::size_t i = 0; i < size; ++i) {
            result += a[i] * b[i];
        }
        return result;
    }
}

// SIMD-optimized axpy (y = alpha*x + y) for float
inline void axpy_simd_float(float* y, float alpha, const float* x, std::size_t size) {
    // Select the optimal version based on the SIMD instruction set
#if defined(SANGI_HAS_AVX512)
    // AVX512 implementation
    constexpr std::size_t simd_size = 16; // 512 bits = 16 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m512 valpha = _mm512_set1_ps(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m512 vx = _mm512_loadu_ps(x + i);
        __m512 vy = _mm512_loadu_ps(y + i);
        // Use the FMA instruction: vy = alpha * vx + vy
        vy = _mm512_fmadd_ps(valpha, vx, vy);
        _mm512_storeu_ps(y + i, vy);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#elif defined(SANGI_HAS_AVX2)
    // AVX2 implementation
    constexpr std::size_t simd_size = 8; // 256 bits = 8 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m256 valpha = _mm256_set1_ps(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m256 vx = _mm256_loadu_ps(x + i);
        __m256 vy = _mm256_loadu_ps(y + i);
        // Use the FMA instruction: vy = alpha * vx + vy
        vy = _mm256_fmadd_ps(valpha, vx, vy);
        _mm256_storeu_ps(y + i, vy);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#elif defined(SANGI_HAS_AVX)
    // AVX implementation
    constexpr std::size_t simd_size = 8; // 256 bits = 8 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m256 valpha = _mm256_set1_ps(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m256 vx = _mm256_loadu_ps(x + i);
        __m256 vy = _mm256_loadu_ps(y + i);
        // vy = alpha * vx + vy
        __m256 vax = _mm256_mul_ps(valpha, vx);
        vy = _mm256_add_ps(vy, vax);
        _mm256_storeu_ps(y + i, vy);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#elif defined(SANGI_HAS_SSE4_1) || defined(SANGI_HAS_SSE2)
    // SSE2/SSE4.1 implementation
    constexpr std::size_t simd_size = 4; // 128 bits = 4 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m128 valpha = _mm_set1_ps(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m128 vx = _mm_loadu_ps(x + i);
        __m128 vy = _mm_loadu_ps(y + i);
        // vy = alpha * vx + vy
        __m128 vax = _mm_mul_ps(valpha, vx);
        vy = _mm_add_ps(vy, vax);
        _mm_storeu_ps(y + i, vy);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#elif defined(SANGI_HAS_NEON)
    // ARM NEON implementation
    constexpr std::size_t simd_size = 4; // 128 bits = 4 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    float32x4_t valpha = vdupq_n_f32(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        float32x4_t vx = vld1q_f32(x + i);
        float32x4_t vy = vld1q_f32(y + i);
        // vy = alpha * vx + vy
        vy = vmlaq_f32(vy, valpha, vx);
        vst1q_f32(y + i, vy);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#else
    // Standard implementation without SIMD
    for (std::size_t i = 0; i < size; ++i) {
        y[i] += alpha * x[i];
    }
#endif
}

// SIMD-optimized axpy (y = alpha*x + y) for double
inline void axpy_simd_double(double* y, double alpha, const double* x, std::size_t size) {
    // Select the optimal version based on the SIMD instruction set
#if defined(SANGI_HAS_AVX512)
    // AVX512 implementation
    constexpr std::size_t simd_size = 8; // 512 bits = 8 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m512d valpha = _mm512_set1_pd(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m512d vx = _mm512_loadu_pd(x + i);
        __m512d vy = _mm512_loadu_pd(y + i);
        // Use the FMA instruction: vy = alpha * vx + vy
        vy = _mm512_fmadd_pd(valpha, vx, vy);
        _mm512_storeu_pd(y + i, vy);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#elif defined(SANGI_HAS_AVX2)
    // AVX2 implementation
    constexpr std::size_t simd_size = 4; // 256 bits = 4 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m256d valpha = _mm256_set1_pd(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m256d vx = _mm256_loadu_pd(x + i);
        __m256d vy = _mm256_loadu_pd(y + i);
        // Use the FMA instruction: vy = alpha * vx + vy
        vy = _mm256_fmadd_pd(valpha, vx, vy);
        _mm256_storeu_pd(y + i, vy);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#elif defined(SANGI_HAS_AVX)
    // AVX implementation
    constexpr std::size_t simd_size = 4; // 256 bits = 4 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m256d valpha = _mm256_set1_pd(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m256d vx = _mm256_loadu_pd(x + i);
        __m256d vy = _mm256_loadu_pd(y + i);
        // vy = alpha * vx + vy
        __m256d vax = _mm256_mul_pd(valpha, vx);
        vy = _mm256_add_pd(vy, vax);
        _mm256_storeu_pd(y + i, vy);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#elif defined(SANGI_HAS_SSE4_1) || defined(SANGI_HAS_SSE2)
    // SSE2/SSE4.1 implementation
    constexpr std::size_t simd_size = 2; // 128 bits = 2 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m128d valpha = _mm_set1_pd(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m128d vx = _mm_loadu_pd(x + i);
        __m128d vy = _mm_loadu_pd(y + i);
        // vy = alpha * vx + vy
        __m128d vax = _mm_mul_pd(valpha, vx);
        vy = _mm_add_pd(vy, vax);
        _mm_storeu_pd(y + i, vy);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#elif defined(SANGI_HAS_NEON)
    // ARM NEON implementation
    constexpr std::size_t simd_size = 2; // 128 bits = 2 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    float64x2_t valpha = vdupq_n_f64(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        float64x2_t vx = vld1q_f64(x + i);
        float64x2_t vy = vld1q_f64(y + i);
        // vy = alpha * vx + vy
        vy = vmlaq_f64(vy, valpha, vx);
        vst1q_f64(y + i, vy);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#else
    // Standard implementation without SIMD
    for (std::size_t i = 0; i < size; ++i) {
        y[i] += alpha * x[i];
    }
#endif
}

// SIMD-optimized axpy (selects the appropriate implementation by type)
template<typename T, typename Alpha>
inline void axpy_simd(T* y, Alpha alpha, const T* x, std::size_t size) {
    if constexpr (std::is_same_v<T, float> && std::is_convertible_v<Alpha, float>) {
        axpy_simd_float(y, static_cast<float>(alpha), x, size);
    } 
    else if constexpr (std::is_same_v<T, double> && std::is_convertible_v<Alpha, double>) {
        axpy_simd_double(y, static_cast<double>(alpha), x, size);
    }
    else {
        // Other types are computed without SIMD optimization
        for (std::size_t i = 0; i < size; ++i) {
            y[i] += static_cast<T>(alpha) * x[i];
        }
    }
}

// SIMD-optimized scaling (x = alpha*x) for float
inline void scale_simd_float(float* x, float alpha, std::size_t size) {
    // Select the optimal version based on the SIMD instruction set
#if defined(SANGI_HAS_AVX512)
    // AVX512 implementation
    constexpr std::size_t simd_size = 16; // 512 bits = 16 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m512 valpha = _mm512_set1_ps(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m512 vx = _mm512_loadu_ps(x + i);
        vx = _mm512_mul_ps(valpha, vx);
        _mm512_storeu_ps(x + i, vx);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        x[i] *= alpha;
    }
    
#elif defined(SANGI_HAS_AVX2) || defined(SANGI_HAS_AVX)
    // AVX / AVX2 implementation
    constexpr std::size_t simd_size = 8; // 256 bits = 8 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m256 valpha = _mm256_set1_ps(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m256 vx = _mm256_loadu_ps(x + i);
        vx = _mm256_mul_ps(valpha, vx);
        _mm256_storeu_ps(x + i, vx);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        x[i] *= alpha;
    }
    
#elif defined(SANGI_HAS_SSE4_1) || defined(SANGI_HAS_SSE2)
    // SSE2/SSE4.1 implementation
    constexpr std::size_t simd_size = 4; // 128 bits = 4 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m128 valpha = _mm_set1_ps(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m128 vx = _mm_loadu_ps(x + i);
        vx = _mm_mul_ps(valpha, vx);
        _mm_storeu_ps(x + i, vx);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        x[i] *= alpha;
    }
    
#elif defined(SANGI_HAS_NEON)
    // ARM NEON implementation
    constexpr std::size_t simd_size = 4; // 128 bits = 4 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    float32x4_t valpha = vdupq_n_f32(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        float32x4_t vx = vld1q_f32(x + i);
        vx = vmulq_f32(valpha, vx);
        vst1q_f32(x + i, vx);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        x[i] *= alpha;
    }
    
#else
    // Standard implementation without SIMD
    for (std::size_t i = 0; i < size; ++i) {
        x[i] *= alpha;
    }
#endif
}

// SIMD-optimized scaling (x = alpha*x) for double
inline void scale_simd_double(double* x, double alpha, std::size_t size) {
    // Select the optimal version based on the SIMD instruction set
#if defined(SANGI_HAS_AVX512)
    // AVX512 implementation
    constexpr std::size_t simd_size = 8; // 512 bits = 8 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m512d valpha = _mm512_set1_pd(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m512d vx = _mm512_loadu_pd(x + i);
        vx = _mm512_mul_pd(valpha, vx);
        _mm512_storeu_pd(x + i, vx);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        x[i] *= alpha;
    }
    
#elif defined(SANGI_HAS_AVX2) || defined(SANGI_HAS_AVX)
    // AVX / AVX2 implementation
    constexpr std::size_t simd_size = 4; // 256 bits = 4 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m256d valpha = _mm256_set1_pd(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m256d vx = _mm256_loadu_pd(x + i);
        vx = _mm256_mul_pd(valpha, vx);
        _mm256_storeu_pd(x + i, vx);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        x[i] *= alpha;
    }
    
#elif defined(SANGI_HAS_SSE4_1) || defined(SANGI_HAS_SSE2)
    // SSE2/SSE4.1 implementation
    constexpr std::size_t simd_size = 2; // 128 bits = 2 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    __m128d valpha = _mm_set1_pd(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m128d vx = _mm_loadu_pd(x + i);
        vx = _mm_mul_pd(valpha, vx);
        _mm_storeu_pd(x + i, vx);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        x[i] *= alpha;
    }
    
#elif defined(SANGI_HAS_NEON)
    // ARM NEON implementation
    constexpr std::size_t simd_size = 2; // 128 bits = 2 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    float64x2_t valpha = vdupq_n_f64(alpha);
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        float64x2_t vx = vld1q_f64(x + i);
        vx = vmulq_f64(valpha, vx);
        vst1q_f64(x + i, vx);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        x[i] *= alpha;
    }
    
#else
    // Standard implementation without SIMD
    for (std::size_t i = 0; i < size; ++i) {
        x[i] *= alpha;
    }
#endif
}

// SIMD-optimized scaling (selects the appropriate implementation by type)
template<typename T, typename Alpha>
inline void scale_simd(T* x, Alpha alpha, std::size_t size) {
    if constexpr (std::is_same_v<T, float> && std::is_convertible_v<Alpha, float>) {
        scale_simd_float(x, static_cast<float>(alpha), size);
    } 
    else if constexpr (std::is_same_v<T, double> && std::is_convertible_v<Alpha, double>) {
        scale_simd_double(x, static_cast<double>(alpha), size);
    }
    else {
        // Other types are computed without SIMD optimization
        for (std::size_t i = 0; i < size; ++i) {
            x[i] *= static_cast<T>(alpha);
        }
    }
}

// SIMD-optimized vector addition (z = x + y) for float
inline void add_simd_float(const float* x, const float* y, float* z, std::size_t size) {
    // Select the optimal version based on the SIMD instruction set
#if defined(SANGI_HAS_AVX512)
    // AVX512 implementation
    constexpr std::size_t simd_size = 16; // 512 bits = 16 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m512 vx = _mm512_loadu_ps(x + i);
        __m512 vy = _mm512_loadu_ps(y + i);
        __m512 vz = _mm512_add_ps(vx, vy);
        _mm512_storeu_ps(z + i, vz);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
    
#elif defined(SANGI_HAS_AVX2) || defined(SANGI_HAS_AVX)
    // AVX / AVX2 implementation
    constexpr std::size_t simd_size = 8; // 256 bits = 8 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m256 vx = _mm256_loadu_ps(x + i);
        __m256 vy = _mm256_loadu_ps(y + i);
        __m256 vz = _mm256_add_ps(vx, vy);
        _mm256_storeu_ps(z + i, vz);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
    
#elif defined(SANGI_HAS_SSE4_1) || defined(SANGI_HAS_SSE2)
    // SSE2/SSE4.1 implementation
    constexpr std::size_t simd_size = 4; // 128 bits = 4 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m128 vx = _mm_loadu_ps(x + i);
        __m128 vy = _mm_loadu_ps(y + i);
        __m128 vz = _mm_add_ps(vx, vy);
        _mm_storeu_ps(z + i, vz);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
    
#elif defined(SANGI_HAS_NEON)
    // ARM NEON implementation
    constexpr std::size_t simd_size = 4; // 128 bits = 4 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        float32x4_t vx = vld1q_f32(x + i);
        float32x4_t vy = vld1q_f32(y + i);
        float32x4_t vz = vaddq_f32(vx, vy);
        vst1q_f32(z + i, vz);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
    
#else
    // Standard implementation without SIMD
    for (std::size_t i = 0; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
#endif
}

// SIMD-optimized vector addition (z = x + y) for double
inline void add_simd_double(const double* x, const double* y, double* z, std::size_t size) {
    // Select the optimal version based on the SIMD instruction set
#if defined(SANGI_HAS_AVX512)
    // AVX512 implementation
    constexpr std::size_t simd_size = 8; // 512 bits = 8 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m512d vx = _mm512_loadu_pd(x + i);
        __m512d vy = _mm512_loadu_pd(y + i);
        __m512d vz = _mm512_add_pd(vx, vy);
        _mm512_storeu_pd(z + i, vz);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
    
#elif defined(SANGI_HAS_AVX2) || defined(SANGI_HAS_AVX)
    // AVX / AVX2 implementation
    constexpr std::size_t simd_size = 4; // 256 bits = 4 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m256d vx = _mm256_loadu_pd(x + i);
        __m256d vy = _mm256_loadu_pd(y + i);
        __m256d vz = _mm256_add_pd(vx, vy);
        _mm256_storeu_pd(z + i, vz);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
    
#elif defined(SANGI_HAS_SSE4_1) || defined(SANGI_HAS_SSE2)
    // SSE2/SSE4.1 implementation
    constexpr std::size_t simd_size = 2; // 128 bits = 2 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m128d vx = _mm_loadu_pd(x + i);
        __m128d vy = _mm_loadu_pd(y + i);
        __m128d vz = _mm_add_pd(vx, vy);
        _mm_storeu_pd(z + i, vz);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
    
#elif defined(SANGI_HAS_NEON)
    // ARM NEON implementation
    constexpr std::size_t simd_size = 2; // 128 bits = 2 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        float64x2_t vx = vld1q_f64(x + i);
        float64x2_t vy = vld1q_f64(y + i);
        float64x2_t vz = vaddq_f64(vx, vy);
        vst1q_f64(z + i, vz);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
    
#else
    // Standard implementation without SIMD
    for (std::size_t i = 0; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
#endif
}

// SIMD-optimized vector addition (selects the appropriate implementation by type)
template<typename T>
inline void add_simd(const T* x, const T* y, T* z, std::size_t size) {
    if constexpr (std::is_same_v<T, float>) {
        add_simd_float(x, y, z, size);
    } 
    else if constexpr (std::is_same_v<T, double>) {
        add_simd_double(x, y, z, size);
    }
    else {
        // Other types are computed without SIMD optimization
        for (std::size_t i = 0; i < size; ++i) {
            z[i] = x[i] + y[i];
        }
    }
}

// SIMD-optimized vector subtraction (z = x - y) for float
inline void subtract_simd_float(const float* x, const float* y, float* z, std::size_t size) {
    // Select the optimal version based on the SIMD instruction set
#if defined(SANGI_HAS_AVX512)
    // AVX512 implementation
    constexpr std::size_t simd_size = 16; // 512 bits = 16 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m512 vx = _mm512_loadu_ps(x + i);
        __m512 vy = _mm512_loadu_ps(y + i);
        __m512 vz = _mm512_sub_ps(vx, vy);
        _mm512_storeu_ps(z + i, vz);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
    
#elif defined(SANGI_HAS_AVX2) || defined(SANGI_HAS_AVX)
    // AVX / AVX2 implementation
    constexpr std::size_t simd_size = 8; // 256 bits = 8 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m256 vx = _mm256_loadu_ps(x + i);
        __m256 vy = _mm256_loadu_ps(y + i);
        __m256 vz = _mm256_sub_ps(vx, vy);
        _mm256_storeu_ps(z + i, vz);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
    
#elif defined(SANGI_HAS_SSE4_1) || defined(SANGI_HAS_SSE2)
    // SSE2/SSE4.1 implementation
    constexpr std::size_t simd_size = 4; // 128 bits = 4 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m128 vx = _mm_loadu_ps(x + i);
        __m128 vy = _mm_loadu_ps(y + i);
        __m128 vz = _mm_sub_ps(vx, vy);
        _mm_storeu_ps(z + i, vz);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
    
#elif defined(SANGI_HAS_NEON)
    // ARM NEON implementation
    constexpr std::size_t simd_size = 4; // 128 bits = 4 floats
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        float32x4_t vx = vld1q_f32(x + i);
        float32x4_t vy = vld1q_f32(y + i);
        float32x4_t vz = vsubq_f32(vx, vy);
        vst1q_f32(z + i, vz);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
    
#else
    // Standard implementation without SIMD
    for (std::size_t i = 0; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
#endif
}

// SIMD-optimized vector subtraction (z = x - y) for double
inline void subtract_simd_double(const double* x, const double* y, double* z, std::size_t size) {
    // Select the optimal version based on the SIMD instruction set
#if defined(SANGI_HAS_AVX512)
    // AVX512 implementation
    constexpr std::size_t simd_size = 8; // 512 bits = 8 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m512d vx = _mm512_loadu_pd(x + i);
        __m512d vy = _mm512_loadu_pd(y + i);
        __m512d vz = _mm512_sub_pd(vx, vy);
        _mm512_storeu_pd(z + i, vz);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
    
#elif defined(SANGI_HAS_AVX2) || defined(SANGI_HAS_AVX)
    // AVX / AVX2 implementation
    constexpr std::size_t simd_size = 4; // 256 bits = 4 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m256d vx = _mm256_loadu_pd(x + i);
        __m256d vy = _mm256_loadu_pd(y + i);
        __m256d vz = _mm256_sub_pd(vx, vy);
        _mm256_storeu_pd(z + i, vz);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
    
#elif defined(SANGI_HAS_SSE4_1) || defined(SANGI_HAS_SSE2)
    // SSE2/SSE4.1 implementation
    constexpr std::size_t simd_size = 2; // 128 bits = 2 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        __m128d vx = _mm_loadu_pd(x + i);
        __m128d vy = _mm_loadu_pd(y + i);
        __m128d vz = _mm_sub_pd(vx, vy);
        _mm_storeu_pd(z + i, vz);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
    
#elif defined(SANGI_HAS_NEON)
    // ARM NEON implementation
    constexpr std::size_t simd_size = 2; // 128 bits = 2 doubles
    
    const std::size_t simd_loop_size = (size / simd_size) * simd_size;
    
    for (std::size_t i = 0; i < simd_loop_size; i += simd_size) {
        float64x2_t vx = vld1q_f64(x + i);
        float64x2_t vy = vld1q_f64(y + i);
        float64x2_t vz = vsubq_f64(vx, vy);
        vst1q_f64(z + i, vz);
    }
    
    // Process the remaining elements
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
    
#else
    // Standard implementation without SIMD
    for (std::size_t i = 0; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
#endif
}

// SIMD-optimized vector subtraction (selects the appropriate implementation by type)
template<typename T>
inline void subtract_simd(const T* x, const T* y, T* z, std::size_t size) {
    if constexpr (std::is_same_v<T, float>) {
        subtract_simd_float(x, y, z, size);
    } 
    else if constexpr (std::is_same_v<T, double>) {
        subtract_simd_double(x, y, z, size);
    }
    else {
        // Other types are computed without SIMD optimization
        for (std::size_t i = 0; i < size; ++i) {
            z[i] = x[i] - y[i];
        }
    }
}

// Other SIMD-optimized functions can also be added

} // namespace simd
} // namespace computation
} // namespace sangi

#endif // SANGI_SIMD_BACKEND_HPP