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

// simd_backend.hpp
//
// SIMD最適化バックエンド実装
//
// このファイルでは、SIMD命令セットを使用した線形代数演算の最適化実装を提供します。
// AVX, AVX2, AVX512, SSE2, SSE4.1などの命令セットを使って、
// ベクトルと行列の演算を高速化します。
//
// 主な機能:
// - 異なるSIMD命令セットのサポート
// - 基本的なベクトル演算の最適化
// - SIMD命令セットの検出と選択

#ifndef CALX_SIMD_BACKEND_HPP
#define CALX_SIMD_BACKEND_HPP

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

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

// SIMD命令セットの検出
#if defined(__AVX512F__)
#define CALX_HAS_AVX512 1
#endif

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

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

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

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

// ARMプラットフォーム用のNEON検出
#if defined(__ARM_NEON) || defined(__ARM_NEON__)
#define CALX_HAS_NEON 1
#endif

// SIMD命令セットのインクルード
#if defined(CALX_HAS_AVX512)
#include <immintrin.h>  // AVX512を含むすべてのインクルード
#elif defined(CALX_HAS_AVX2)
#include <immintrin.h>  // AVX2用
#elif defined(CALX_HAS_AVX)
#include <immintrin.h>  // AVX用
#elif defined(CALX_HAS_SSE4_1)
#include <smmintrin.h>  // SSE4.1用
#elif defined(CALX_HAS_SSE2)
#include <emmintrin.h>  // SSE2用
#endif

#if defined(CALX_HAS_NEON)
#include <arm_neon.h>   // ARM NEON用
#endif

namespace calx {
namespace computation {
namespace simd {

// SIMD命令セットの列挙型
enum class SimdInstructionSet {
    None,
    SSE2,
    SSE4_1,
    AVX,
    AVX2,
    AVX512,
    NEON
};

// 利用可能なSIMD命令セットの検出
constexpr SimdInstructionSet detect_simd_instruction_set() {
#if defined(CALX_HAS_AVX512)
    return SimdInstructionSet::AVX512;
#elif defined(CALX_HAS_AVX2)
    return SimdInstructionSet::AVX2;
#elif defined(CALX_HAS_AVX)
    return SimdInstructionSet::AVX;
#elif defined(CALX_HAS_SSE4_1)
    return SimdInstructionSet::SSE4_1;
#elif defined(CALX_HAS_SSE2)
    return SimdInstructionSet::SSE2;
#elif defined(CALX_HAS_NEON)
    return SimdInstructionSet::NEON;
#else
    return SimdInstructionSet::None;
#endif
}

// 利用可能なSIMD命令セット
constexpr SimdInstructionSet available_simd_level = detect_simd_instruction_set();

//-----------------------------------------------------------------------------
// SIMD最適化ヘルパー関数
//-----------------------------------------------------------------------------

// SIMD最適化された内積計算 (float)
inline float dot_product_simd_float(const float* a, const float* b, std::size_t size) {
    float result = 0.0f;
    
    // SIMD命令セットに基づいて最適なバージョンを選択
#if defined(CALX_HAS_AVX512)
    // AVX512実装
    constexpr std::size_t simd_size = 16; // 512ビット = 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);
    }
    
    // 水平加算で合計を計算
    result = _mm512_reduce_add_ps(sum);
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#elif defined(CALX_HAS_AVX2)
    // AVX2実装
    constexpr std::size_t simd_size = 8; // 256ビット = 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);
        // FMA命令を使用: sum += va * vb
        sum = _mm256_fmadd_ps(va, vb, sum);
    }
    
    // 水平加算で合計を計算
    __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);
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#elif defined(CALX_HAS_AVX)
    // AVX実装
    constexpr std::size_t simd_size = 8; // 256ビット = 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);
    }
    
    // 水平加算で合計を計算
    __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);
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#elif defined(CALX_HAS_SSE4_1) || defined(CALX_HAS_SSE2)
    // SSE2/SSE4.1実装
    constexpr std::size_t simd_size = 4; // 128ビット = 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);
    }
    
    // 水平加算で合計を計算
    __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);
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#elif defined(CALX_HAS_NEON)
    // ARM NEON実装
    constexpr std::size_t simd_size = 4; // 128ビット = 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);
        // 積和演算: sum += va * vb
        sum = vmlaq_f32(sum, va, vb);
    }
    
    // 水平加算で合計を計算
    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);
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#else
    // SIMDなしの標準実装
    for (std::size_t i = 0; i < size; ++i) {
        result += a[i] * b[i];
    }
#endif
    
    return result;
}

// SIMD最適化された内積計算 (double)
inline double dot_product_simd_double(const double* a, const double* b, std::size_t size) {
    double result = 0.0;
    
    // SIMD命令セットに基づいて最適なバージョンを選択
#if defined(CALX_HAS_AVX512)
    // AVX512実装
    constexpr std::size_t simd_size = 8; // 512ビット = 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);
    }
    
    // 水平加算で合計を計算
    result = _mm512_reduce_add_pd(sum);
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#elif defined(CALX_HAS_AVX2)
    // AVX2実装
    constexpr std::size_t simd_size = 4; // 256ビット = 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);
        // FMA命令を使用: sum += va * vb
        sum = _mm256_fmadd_pd(va, vb, sum);
    }
    
    // 水平加算で合計を計算
    __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);
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#elif defined(CALX_HAS_AVX)
    // AVX実装
    constexpr std::size_t simd_size = 4; // 256ビット = 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);
    }
    
    // 水平加算で合計を計算
    __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);
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#elif defined(CALX_HAS_SSE4_1) || defined(CALX_HAS_SSE2)
    // SSE2/SSE4.1実装
    constexpr std::size_t simd_size = 2; // 128ビット = 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);
    }
    
    // 水平加算で合計を計算
    __m128d sum_1 = _mm_add_sd(sum, _mm_unpackhi_pd(sum, sum));
    
    result = _mm_cvtsd_f64(sum_1);
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#elif defined(CALX_HAS_NEON)
    // ARM NEON実装 (注: ARMv8以降でdouble精度のNEON命令がサポートされています)
    constexpr std::size_t simd_size = 2; // 128ビット = 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);
        // 積和演算: sum += va * vb
        sum = vmlaq_f64(sum, va, vb);
    }
    
    // 水平加算で合計を計算
    result = vgetq_lane_f64(sum, 0) + vgetq_lane_f64(sum, 1);
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        result += a[i] * b[i];
    }
    
#else
    // SIMDなしの標準実装
    for (std::size_t i = 0; i < size; ++i) {
        result += a[i] * b[i];
    }
#endif
    
    return result;
}

// SIMD最適化された内積計算（型に基づいて適切な実装を選択）
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 {
        // その他の型はSIMD最適化なしで計算
        T result = numeric_traits<T>::zero();
        for (std::size_t i = 0; i < size; ++i) {
            result += a[i] * b[i];
        }
        return result;
    }
}

// SIMD最適化されたaxpy演算 (y = alpha*x + y) for float
inline void axpy_simd_float(float* y, float alpha, const float* x, std::size_t size) {
    // SIMD命令セットに基づいて最適なバージョンを選択
#if defined(CALX_HAS_AVX512)
    // AVX512実装
    constexpr std::size_t simd_size = 16; // 512ビット = 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);
        // FMA命令を使用: vy = alpha * vx + vy
        vy = _mm512_fmadd_ps(valpha, vx, vy);
        _mm512_storeu_ps(y + i, vy);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#elif defined(CALX_HAS_AVX2)
    // AVX2実装
    constexpr std::size_t simd_size = 8; // 256ビット = 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);
        // FMA命令を使用: vy = alpha * vx + vy
        vy = _mm256_fmadd_ps(valpha, vx, vy);
        _mm256_storeu_ps(y + i, vy);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#elif defined(CALX_HAS_AVX)
    // AVX実装
    constexpr std::size_t simd_size = 8; // 256ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#elif defined(CALX_HAS_SSE4_1) || defined(CALX_HAS_SSE2)
    // SSE2/SSE4.1実装
    constexpr std::size_t simd_size = 4; // 128ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#elif defined(CALX_HAS_NEON)
    // ARM NEON実装
    constexpr std::size_t simd_size = 4; // 128ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#else
    // SIMDなしの標準実装
    for (std::size_t i = 0; i < size; ++i) {
        y[i] += alpha * x[i];
    }
#endif
}

// SIMD最適化されたaxpy演算 (y = alpha*x + y) for double
inline void axpy_simd_double(double* y, double alpha, const double* x, std::size_t size) {
    // SIMD命令セットに基づいて最適なバージョンを選択
#if defined(CALX_HAS_AVX512)
    // AVX512実装
    constexpr std::size_t simd_size = 8; // 512ビット = 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);
        // FMA命令を使用: vy = alpha * vx + vy
        vy = _mm512_fmadd_pd(valpha, vx, vy);
        _mm512_storeu_pd(y + i, vy);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#elif defined(CALX_HAS_AVX2)
    // AVX2実装
    constexpr std::size_t simd_size = 4; // 256ビット = 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);
        // FMA命令を使用: vy = alpha * vx + vy
        vy = _mm256_fmadd_pd(valpha, vx, vy);
        _mm256_storeu_pd(y + i, vy);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#elif defined(CALX_HAS_AVX)
    // AVX実装
    constexpr std::size_t simd_size = 4; // 256ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#elif defined(CALX_HAS_SSE4_1) || defined(CALX_HAS_SSE2)
    // SSE2/SSE4.1実装
    constexpr std::size_t simd_size = 2; // 128ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#elif defined(CALX_HAS_NEON)
    // ARM NEON実装
    constexpr std::size_t simd_size = 2; // 128ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        y[i] += alpha * x[i];
    }
    
#else
    // SIMDなしの標準実装
    for (std::size_t i = 0; i < size; ++i) {
        y[i] += alpha * x[i];
    }
#endif
}

// SIMD最適化されたaxpy演算（型に基づいて適切な実装を選択）
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 {
        // その他の型はSIMD最適化なしで計算
        for (std::size_t i = 0; i < size; ++i) {
            y[i] += static_cast<T>(alpha) * x[i];
        }
    }
}

// SIMD最適化されたスケーリング演算 (x = alpha*x) for float
inline void scale_simd_float(float* x, float alpha, std::size_t size) {
    // SIMD命令セットに基づいて最適なバージョンを選択
#if defined(CALX_HAS_AVX512)
    // AVX512実装
    constexpr std::size_t simd_size = 16; // 512ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        x[i] *= alpha;
    }
    
#elif defined(CALX_HAS_AVX2) || defined(CALX_HAS_AVX)
    // AVX/AVX2実装
    constexpr std::size_t simd_size = 8; // 256ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        x[i] *= alpha;
    }
    
#elif defined(CALX_HAS_SSE4_1) || defined(CALX_HAS_SSE2)
    // SSE2/SSE4.1実装
    constexpr std::size_t simd_size = 4; // 128ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        x[i] *= alpha;
    }
    
#elif defined(CALX_HAS_NEON)
    // ARM NEON実装
    constexpr std::size_t simd_size = 4; // 128ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        x[i] *= alpha;
    }
    
#else
    // SIMDなしの標準実装
    for (std::size_t i = 0; i < size; ++i) {
        x[i] *= alpha;
    }
#endif
}

// SIMD最適化されたスケーリング演算 (x = alpha*x) for double
inline void scale_simd_double(double* x, double alpha, std::size_t size) {
    // SIMD命令セットに基づいて最適なバージョンを選択
#if defined(CALX_HAS_AVX512)
    // AVX512実装
    constexpr std::size_t simd_size = 8; // 512ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        x[i] *= alpha;
    }
    
#elif defined(CALX_HAS_AVX2) || defined(CALX_HAS_AVX)
    // AVX/AVX2実装
    constexpr std::size_t simd_size = 4; // 256ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        x[i] *= alpha;
    }
    
#elif defined(CALX_HAS_SSE4_1) || defined(CALX_HAS_SSE2)
    // SSE2/SSE4.1実装
    constexpr std::size_t simd_size = 2; // 128ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        x[i] *= alpha;
    }
    
#elif defined(CALX_HAS_NEON)
    // ARM NEON実装
    constexpr std::size_t simd_size = 2; // 128ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        x[i] *= alpha;
    }
    
#else
    // SIMDなしの標準実装
    for (std::size_t i = 0; i < size; ++i) {
        x[i] *= alpha;
    }
#endif
}

// SIMD最適化されたスケーリング演算（型に基づいて適切な実装を選択）
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 {
        // その他の型はSIMD最適化なしで計算
        for (std::size_t i = 0; i < size; ++i) {
            x[i] *= static_cast<T>(alpha);
        }
    }
}

// SIMD最適化されたベクトル加算 (z = x + y) for float
inline void add_simd_float(const float* x, const float* y, float* z, std::size_t size) {
    // SIMD命令セットに基づいて最適なバージョンを選択
#if defined(CALX_HAS_AVX512)
    // AVX512実装
    constexpr std::size_t simd_size = 16; // 512ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
    
#elif defined(CALX_HAS_AVX2) || defined(CALX_HAS_AVX)
    // AVX/AVX2実装
    constexpr std::size_t simd_size = 8; // 256ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
    
#elif defined(CALX_HAS_SSE4_1) || defined(CALX_HAS_SSE2)
    // SSE2/SSE4.1実装
    constexpr std::size_t simd_size = 4; // 128ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
    
#elif defined(CALX_HAS_NEON)
    // ARM NEON実装
    constexpr std::size_t simd_size = 4; // 128ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
    
#else
    // SIMDなしの標準実装
    for (std::size_t i = 0; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
#endif
}

// SIMD最適化されたベクトル加算 (z = x + y) for double
inline void add_simd_double(const double* x, const double* y, double* z, std::size_t size) {
    // SIMD命令セットに基づいて最適なバージョンを選択
#if defined(CALX_HAS_AVX512)
    // AVX512実装
    constexpr std::size_t simd_size = 8; // 512ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
    
#elif defined(CALX_HAS_AVX2) || defined(CALX_HAS_AVX)
    // AVX/AVX2実装
    constexpr std::size_t simd_size = 4; // 256ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
    
#elif defined(CALX_HAS_SSE4_1) || defined(CALX_HAS_SSE2)
    // SSE2/SSE4.1実装
    constexpr std::size_t simd_size = 2; // 128ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
    
#elif defined(CALX_HAS_NEON)
    // ARM NEON実装
    constexpr std::size_t simd_size = 2; // 128ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
    
#else
    // SIMDなしの標準実装
    for (std::size_t i = 0; i < size; ++i) {
        z[i] = x[i] + y[i];
    }
#endif
}

// SIMD最適化されたベクトル加算（型に基づいて適切な実装を選択）
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 {
        // その他の型はSIMD最適化なしで計算
        for (std::size_t i = 0; i < size; ++i) {
            z[i] = x[i] + y[i];
        }
    }
}

// SIMD最適化されたベクトル減算 (z = x - y) for float
inline void subtract_simd_float(const float* x, const float* y, float* z, std::size_t size) {
    // SIMD命令セットに基づいて最適なバージョンを選択
#if defined(CALX_HAS_AVX512)
    // AVX512実装
    constexpr std::size_t simd_size = 16; // 512ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
    
#elif defined(CALX_HAS_AVX2) || defined(CALX_HAS_AVX)
    // AVX/AVX2実装
    constexpr std::size_t simd_size = 8; // 256ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
    
#elif defined(CALX_HAS_SSE4_1) || defined(CALX_HAS_SSE2)
    // SSE2/SSE4.1実装
    constexpr std::size_t simd_size = 4; // 128ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
    
#elif defined(CALX_HAS_NEON)
    // ARM NEON実装
    constexpr std::size_t simd_size = 4; // 128ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
    
#else
    // SIMDなしの標準実装
    for (std::size_t i = 0; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
#endif
}

// SIMD最適化されたベクトル減算 (z = x - y) for double
inline void subtract_simd_double(const double* x, const double* y, double* z, std::size_t size) {
    // SIMD命令セットに基づいて最適なバージョンを選択
#if defined(CALX_HAS_AVX512)
    // AVX512実装
    constexpr std::size_t simd_size = 8; // 512ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
    
#elif defined(CALX_HAS_AVX2) || defined(CALX_HAS_AVX)
    // AVX/AVX2実装
    constexpr std::size_t simd_size = 4; // 256ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
    
#elif defined(CALX_HAS_SSE4_1) || defined(CALX_HAS_SSE2)
    // SSE2/SSE4.1実装
    constexpr std::size_t simd_size = 2; // 128ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
    
#elif defined(CALX_HAS_NEON)
    // ARM NEON実装
    constexpr std::size_t simd_size = 2; // 128ビット = 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);
    }
    
    // 残りの要素を処理
    for (std::size_t i = simd_loop_size; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
    
#else
    // SIMDなしの標準実装
    for (std::size_t i = 0; i < size; ++i) {
        z[i] = x[i] - y[i];
    }
#endif
}

// SIMD最適化されたベクトル減算（型に基づいて適切な実装を選択）
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 {
        // その他の型はSIMD最適化なしで計算
        for (std::size_t i = 0; i < size; ++i) {
            z[i] = x[i] - y[i];
        }
    }
}

// 他のSIMD最適化関数も追加可能

} // namespace simd
} // namespace computation
} // namespace calx

#endif // CALX_SIMD_BACKEND_HPP