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

// computation_policy.hpp
#pragma once

#include <cstddef>
#include <cstring>
#include <type_traits>
#include <algorithm>
#include <future>
#include <math/core/numeric_state.hpp>
#include <math/core/computation/simd_backend.hpp>
#include <math/core/simd_traits.hpp>
#include <math/core/mp/ThreadPool.hpp>

namespace calx {

    // 計算ポリシーのベースクラス
    template <typename T>
    class ComputePolicyBase {
    public:
        using value_type = T;
        using size_type = std::size_t;
    };

    // 前方宣言
    template <typename T> class DefaultComputePolicy;
    template <typename T> class SIMDComputePolicy;
    template <typename T> class MKLComputePolicy;

    // デフォルト計算ポリシー
    template <typename T>
    class DefaultComputePolicy : public ComputePolicyBase<T> {
    public:
        using value_type = T;
        using size_type = std::size_t;

        // ベクトル加算
        template <typename VecA, typename VecB, typename VecResult>
        static void vector_add(const VecA& a, const VecB& b, VecResult& result) {
            if (a.size() != b.size() || a.size() != result.size()) {
                throw std::invalid_argument("Vector dimensions mismatch for addition");
            }

            for (size_type i = 0; i < a.size(); ++i) {
                result[i] = a[i] + b[i];

                // NumericState システムの適用
                if (numeric_state_traits<T>::is_supported) {
                    if (numeric_state_traits<T>::isSpecial(a[i]) ||
                        numeric_state_traits<T>::isSpecial(b[i])) {
                        // NumericState伝播処理
                    }
                }
            }
        }

        // ベクトル減算
        template <typename VecA, typename VecB, typename VecResult>
        static void vector_subtract(const VecA& a, const VecB& b, VecResult& result) {
            if (a.size() != b.size() || a.size() != result.size()) {
                throw std::invalid_argument("Vector dimensions mismatch for subtraction");
            }

            for (size_type i = 0; i < a.size(); ++i) {
                result[i] = a[i] - b[i];

                // NumericState システムの適用
                if (numeric_state_traits<T>::is_supported) {
                    if (numeric_state_traits<T>::isSpecial(a[i]) ||
                        numeric_state_traits<T>::isSpecial(b[i])) {
                        // NumericState伝播処理
                    }
                }
            }
        }

        // ベクトル内積
        template <typename VecA, typename VecB>
        static T dot(const VecA& a, const VecB& b) {
            if (a.size() != b.size()) {
                throw std::invalid_argument("Vector dimensions mismatch for dot product");
            }

            T result = static_cast<T>(0);
            for (size_type i = 0; i < a.size(); ++i) {
                result += a[i] * b[i];

                // NumericState システムの適用
                if (numeric_state_traits<T>::is_supported) {
                    if (numeric_state_traits<T>::isSpecial(a[i]) ||
                        numeric_state_traits<T>::isSpecial(b[i])) {
                        // NumericState伝播処理
                    }
                }
            }

            return result;
        }

        // ベクトルノルム（L2ノルム）
        template <typename Vec>
        static T norm(const Vec& vec) {
            T sum_of_squares = static_cast<T>(0);

            for (size_type i = 0; i < vec.size(); ++i) {
                sum_of_squares += vec[i] * vec[i];

                // NumericState システムの適用
                if (numeric_state_traits<T>::is_supported) {
                    if (numeric_state_traits<T>::isSpecial(vec[i])) {
                        // NumericState伝播処理
                    }
                }
            }

            return std::sqrt(sum_of_squares);
        }

        // ベクトルスカラー乗算
        template <typename Vec, typename VecResult>
        static void vector_scalar_multiply(const Vec& vec, const T& scalar, VecResult& result) {
            if (vec.size() != result.size()) {
                throw std::invalid_argument("Vector dimensions mismatch for scalar multiplication");
            }

            for (size_type i = 0; i < vec.size(); ++i) {
                result[i] = vec[i] * scalar;

                // NumericState システムの適用
                if (numeric_state_traits<T>::is_supported) {
                    if (numeric_state_traits<T>::isSpecial(vec[i]) ||
                        numeric_state_traits<T>::isSpecial(scalar)) {
                        // NumericState伝播処理
                    }
                }
            }
        }

        // ベクトル要素ごとの乗算
        template <typename VecA, typename VecB, typename VecResult>
        static void elementWiseMultiply(const VecA& a, const VecB& b, VecResult& result) {
            if (a.size() != b.size() || a.size() != result.size()) {
                throw std::invalid_argument("Vector dimensions mismatch for element-wise multiplication");
            }

            for (size_type i = 0; i < a.size(); ++i) {
                result[i] = a[i] * b[i];

                // NumericState システムの適用
                if (numeric_state_traits<T>::is_supported) {
                    if (numeric_state_traits<T>::isSpecial(a[i]) ||
                        numeric_state_traits<T>::isSpecial(b[i])) {
                        // NumericState伝播処理
                    }
                }
            }
        }

        // ベクトル要素ごとの除算
        template <typename VecA, typename VecB, typename VecResult>
        static void elementWiseDivide(const VecA& a, const VecB& b, VecResult& result) {
            if (a.size() != b.size() || a.size() != result.size()) {
                throw std::invalid_argument("Vector dimensions mismatch for element-wise division");
            }

            for (size_type i = 0; i < a.size(); ++i) {
                if (b[i] == static_cast<T>(0)) {
                    // ゼロ除算を処理（特殊値を設定するか例外をスローする）
                    if (numeric_state_traits<T>::is_supported && numeric_state_traits<T>::isSpecial(a[i])) {
                        // 特殊状態を設定
                    }
                    else {
                        throw std::invalid_argument("Division by zero");
                    }
                }
                else {
                    result[i] = a[i] / b[i];

                    // NumericState システムの適用
                    if (numeric_state_traits<T>::is_supported) {
                        if (numeric_state_traits<T>::isSpecial(a[i]) ||
                            numeric_state_traits<T>::isSpecial(b[i])) {
                            // NumericState伝播処理
                        }
                    }
                }
            }
        }

        // 行列操作のための関数

        // 行列加算
        template <typename MatA, typename MatB, typename MatResult>
        static void add(const MatA& a, const MatB& b, MatResult& result) {
            if (a.rows() != b.rows() || a.cols() != b.cols() ||
                a.rows() != result.rows() || a.cols() != result.cols()) {
                throw std::invalid_argument("Matrix dimensions mismatch for addition");
            }

            for (size_type i = 0; i < a.rows(); ++i) {
                for (size_type j = 0; j < a.cols(); ++j) {
                    result(i, j) = a(i, j) + b(i, j);
                }
            }
        }

        // 行列減算
        template <typename MatA, typename MatB, typename MatResult>
        static void subtract(const MatA& a, const MatB& b, MatResult& result) {
            if (a.rows() != b.rows() || a.cols() != b.cols() ||
                a.rows() != result.rows() || a.cols() != result.cols()) {
                throw std::invalid_argument("Matrix dimensions mismatch for subtraction");
            }

            for (size_type i = 0; i < a.rows(); ++i) {
                for (size_type j = 0; j < a.cols(); ++j) {
                    result(i, j) = a(i, j) - b(i, j);
                }
            }
        }

        // 行列乗算
        template <typename MatA, typename MatB, typename MatResult>
        static void multiply(const MatA& a, const MatB& b, MatResult& result) {
            if (a.cols() != b.rows() || result.rows() != a.rows() || result.cols() != b.cols()) {
                throw std::invalid_argument("Matrix dimensions mismatch for multiplication");
            }

            const size_type M = a.rows();
            const size_type N = b.cols();
            const size_type K = a.cols();

            // ゼロ初期化
            T* c_data = result.data();
            const T zero_val = numeric_traits<T>::zero();
            for (size_type i = 0; i < M * N; ++i)
                c_data[i] = zero_val;

            // raw ポインタアクセス (行優先レイアウト前提)
            const T* a_data = a.data();
            const T* b_data = b.data();
            const size_type lda = a.cols();
            const size_type ldb = b.cols();
            const size_type ldc = result.cols();

            gemm_blocked(a_data, b_data, c_data, M, N, K, lda, ldb, ldc);
        }

    private:
        // ================================================================
        // ブロック gemm: タイリング + レジスタブロッキング + B パッキング
        // ================================================================
        //
        // MR × NR マイクロカーネル: C を YMM レジスタに保持し k ループで FMA 累積
        // → C の load/store が k あたり 1 回から全体で 1 回に削減 (BK 倍の帯域削減)
        //
        // B パッキング: B パネルを NR 幅の列パネルに再配置
        // → k 方向のストライドアクセス (ldb 飛び) を連続アクセスに変換

        static constexpr size_type MR = 4;   // レジスタブロック行数

        // NR: レジスタブロック列数 (SIMD パケット 2 本分)
        static constexpr size_type gemm_NR() {
            if constexpr (has_simd_packet_v<T>)
                return PacketTraits<T>::size * 2;  // double: 8, float: 16
            else
                return 4;
        }

        static constexpr size_type BM = 64;    // 行ブロック
        static constexpr size_type BN = 64;    // 列ブロック (NR の倍数)
        static constexpr size_type BK = 128;   // 内積ブロック (レジスタブロッキングで増やせる)

        // B パッキング: B[k0:k0+bk, j0:j0+bn] を NR 幅列パネルに再配置
        // 出力レイアウト: panel p の k 番目 → packB[p * bk * NR + k * NR + 0..NR-1]
        static void pack_b_panel(const T* b, T* packB,
                                  size_type k0, size_type bk,
                                  size_type j0, size_type bn,
                                  size_type ldb) {
            const size_type NR = gemm_NR();
            const size_type n_full = bn / NR;
            // フルパネル
            for (size_type p = 0; p < n_full; ++p) {
                T* dst = packB + p * bk * NR;
                for (size_type k = 0; k < bk; ++k) {
                    const T* src = &b[(k0 + k) * ldb + j0 + p * NR];
                    std::memcpy(dst + k * NR, src, NR * sizeof(T));
                }
            }
            // 残端パネル (NR 未満) — ゼロパディング
            const size_type rem = bn % NR;
            if (rem > 0) {
                T* dst = packB + n_full * bk * NR;
                for (size_type k = 0; k < bk; ++k) {
                    const T* src = &b[(k0 + k) * ldb + j0 + n_full * NR];
                    size_type r = 0;
                    for (; r < rem; ++r)
                        dst[k * NR + r] = src[r];
                    for (; r < NR; ++r)
                        dst[k * NR + r] = T{0};
                }
            }
        }

        // レジスタブロック マイクロカーネル (packed B 版): C[i:i+4, 0:NR] += A * packB
        // C は 8 本の YMM レジスタで保持、k ループ全体で累積
        static void gemm_micro_kernel(const T* a, const T* packB_panel, T* c,
                                       size_type i, size_type bk,
                                       size_type lda, size_type ldc) {
            if constexpr (has_simd_packet_v<T>) {
                using PT = PacketTraits<T>;
                constexpr size_type PS = PT::size;
                constexpr size_type NR = PS * 2;

                auto c00 = PT::load(&c[(i+0)*ldc]);
                auto c01 = PT::load(&c[(i+0)*ldc + PS]);
                auto c10 = PT::load(&c[(i+1)*ldc]);
                auto c11 = PT::load(&c[(i+1)*ldc + PS]);
                auto c20 = PT::load(&c[(i+2)*ldc]);
                auto c21 = PT::load(&c[(i+2)*ldc + PS]);
                auto c30 = PT::load(&c[(i+3)*ldc]);
                auto c31 = PT::load(&c[(i+3)*ldc + PS]);

                for (size_type k = 0; k < bk; ++k) {
                    const auto b0 = PT::load(&packB_panel[k * NR]);
                    const auto b1 = PT::load(&packB_panel[k * NR + PS]);

                    const auto a0 = PT::set1(a[(i+0)*lda + k]);
                    c00 = PT::fmadd(a0, b0, c00);
                    c01 = PT::fmadd(a0, b1, c01);

                    const auto a1 = PT::set1(a[(i+1)*lda + k]);
                    c10 = PT::fmadd(a1, b0, c10);
                    c11 = PT::fmadd(a1, b1, c11);

                    const auto a2 = PT::set1(a[(i+2)*lda + k]);
                    c20 = PT::fmadd(a2, b0, c20);
                    c21 = PT::fmadd(a2, b1, c21);

                    const auto a3 = PT::set1(a[(i+3)*lda + k]);
                    c30 = PT::fmadd(a3, b0, c30);
                    c31 = PT::fmadd(a3, b1, c31);
                }

                PT::store(&c[(i+0)*ldc], c00);
                PT::store(&c[(i+0)*ldc + PS], c01);
                PT::store(&c[(i+1)*ldc], c10);
                PT::store(&c[(i+1)*ldc + PS], c11);
                PT::store(&c[(i+2)*ldc], c20);
                PT::store(&c[(i+2)*ldc + PS], c21);
                PT::store(&c[(i+3)*ldc], c30);
                PT::store(&c[(i+3)*ldc + PS], c31);
            }
        }

        // レジスタブロック マイクロカーネル (直接 B 版): パッキングなし
        // 小行列でパッキングオーバーヘッドを回避
        static void gemm_micro_kernel_direct(const T* a, const T* b, T* c,
                                              size_type i, size_type j,
                                              size_type k0, size_type bk,
                                              size_type lda, size_type ldb, size_type ldc) {
            if constexpr (has_simd_packet_v<T>) {
                using PT = PacketTraits<T>;
                constexpr size_type PS = PT::size;

                auto c00 = PT::load(&c[(i+0)*ldc + j]);
                auto c01 = PT::load(&c[(i+0)*ldc + j + PS]);
                auto c10 = PT::load(&c[(i+1)*ldc + j]);
                auto c11 = PT::load(&c[(i+1)*ldc + j + PS]);
                auto c20 = PT::load(&c[(i+2)*ldc + j]);
                auto c21 = PT::load(&c[(i+2)*ldc + j + PS]);
                auto c30 = PT::load(&c[(i+3)*ldc + j]);
                auto c31 = PT::load(&c[(i+3)*ldc + j + PS]);

                for (size_type k = k0; k < k0 + bk; ++k) {
                    const auto b0 = PT::load(&b[k * ldb + j]);
                    const auto b1 = PT::load(&b[k * ldb + j + PS]);

                    const auto a0 = PT::set1(a[(i+0)*lda + k]);
                    c00 = PT::fmadd(a0, b0, c00);
                    c01 = PT::fmadd(a0, b1, c01);

                    const auto a1 = PT::set1(a[(i+1)*lda + k]);
                    c10 = PT::fmadd(a1, b0, c10);
                    c11 = PT::fmadd(a1, b1, c11);

                    const auto a2 = PT::set1(a[(i+2)*lda + k]);
                    c20 = PT::fmadd(a2, b0, c20);
                    c21 = PT::fmadd(a2, b1, c21);

                    const auto a3 = PT::set1(a[(i+3)*lda + k]);
                    c30 = PT::fmadd(a3, b0, c30);
                    c31 = PT::fmadd(a3, b1, c31);
                }

                PT::store(&c[(i+0)*ldc + j], c00);
                PT::store(&c[(i+0)*ldc + j + PS], c01);
                PT::store(&c[(i+1)*ldc + j], c10);
                PT::store(&c[(i+1)*ldc + j + PS], c11);
                PT::store(&c[(i+2)*ldc + j], c20);
                PT::store(&c[(i+2)*ldc + j + PS], c21);
                PT::store(&c[(i+3)*ldc + j], c30);
                PT::store(&c[(i+3)*ldc + j + PS], c31);
            }
        }

        // 小行列パス: パッキングなし、レジスタブロッキングのみ
        // パッキングのオーバーヘッドが演算量に対して大きい場合に使用
        static void gemm_direct(const T* a, const T* b, T* c,
                                 size_type M, size_type N, size_type K,
                                 size_type lda, size_type ldb, size_type ldc) {
            if constexpr (has_simd_packet_v<T>) {
                using PT = PacketTraits<T>;
                constexpr size_type PS = PT::size;
                constexpr size_type NR = PS * 2;

                // MR×NR レジスタブロック (B を直接読み出し)
                size_type i = 0;
                for (; i + MR <= M; i += MR) {
                    size_type j = 0;
                    for (; j + NR <= N; j += NR) {
                        gemm_micro_kernel_direct(a, b, c, i, j, 0, K, lda, ldb, ldc);
                    }
                    // 残端列: 1 パケット幅
                    for (; j + PS <= N; j += PS) {
                        auto c00 = PT::load(&c[(i+0)*ldc + j]);
                        auto c10 = PT::load(&c[(i+1)*ldc + j]);
                        auto c20 = PT::load(&c[(i+2)*ldc + j]);
                        auto c30 = PT::load(&c[(i+3)*ldc + j]);
                        for (size_type k = 0; k < K; ++k) {
                            const auto bv = PT::load(&b[k * ldb + j]);
                            c00 = PT::fmadd(PT::set1(a[(i+0)*lda + k]), bv, c00);
                            c10 = PT::fmadd(PT::set1(a[(i+1)*lda + k]), bv, c10);
                            c20 = PT::fmadd(PT::set1(a[(i+2)*lda + k]), bv, c20);
                            c30 = PT::fmadd(PT::set1(a[(i+3)*lda + k]), bv, c30);
                        }
                        PT::store(&c[(i+0)*ldc + j], c00);
                        PT::store(&c[(i+1)*ldc + j], c10);
                        PT::store(&c[(i+2)*ldc + j], c20);
                        PT::store(&c[(i+3)*ldc + j], c30);
                    }
                    // スカラー残端
                    for (; j < N; ++j) {
                        T s0 = c[(i+0)*ldc+j], s1 = c[(i+1)*ldc+j];
                        T s2 = c[(i+2)*ldc+j], s3 = c[(i+3)*ldc+j];
                        for (size_type k = 0; k < K; ++k) {
                            T bkj = b[k*ldb+j];
                            s0 += a[(i+0)*lda+k] * bkj; s1 += a[(i+1)*lda+k] * bkj;
                            s2 += a[(i+2)*lda+k] * bkj; s3 += a[(i+3)*lda+k] * bkj;
                        }
                        c[(i+0)*ldc+j] = s0; c[(i+1)*ldc+j] = s1;
                        c[(i+2)*ldc+j] = s2; c[(i+3)*ldc+j] = s3;
                    }
                }
                // 残端行
                for (; i < M; ++i) {
                    for (size_type k = 0; k < K; ++k) {
                        const auto a_ik = PT::set1(a[i*lda + k]);
                        T a_val = a[i*lda + k];
                        size_type j = 0;
                        const size_type j_vec_end = (N / PS) * PS;
                        for (; j < j_vec_end; j += PS) {
                            auto cv = PT::load(&c[i*ldc + j]);
                            PT::store(&c[i*ldc + j], PT::fmadd(a_ik, PT::load(&b[k*ldb + j]), cv));
                        }
                        for (; j < N; ++j)
                            c[i*ldc + j] += a_val * b[k*ldb + j];
                    }
                }
            } else {
                for (size_type i = 0; i < M; ++i)
                    for (size_type k = 0; k < K; ++k) {
                        T a_val = a[i*lda + k];
                        for (size_type j = 0; j < N; ++j)
                            c[i*ldc + j] += a_val * b[k*ldb + j];
                    }
            }
        }

        // 大行列パス: B パッキング + レジスタブロッキング
        static void gemm_packed(const T* a, const T* b, T* c,
                                 size_type M, size_type N, size_type K,
                                 size_type lda, size_type ldb, size_type ldc) {
            if constexpr (has_simd_packet_v<T>) {
                using PT = PacketTraits<T>;
                constexpr size_type PS = PT::size;
                constexpr size_type NR = PS * 2;

                constexpr size_type PACK_B_SIZE = (BN / NR + 1) * BK * NR;
                alignas(32) T packB[PACK_B_SIZE];

                for (size_type jj = 0; jj < N; jj += BN) {
                    const size_type bn = ((jj + BN) < N) ? BN : (N - jj);
                    for (size_type kk = 0; kk < K; kk += BK) {
                        const size_type bk = ((kk + BK) < K) ? BK : (K - kk);

                        pack_b_panel(b, packB, kk, bk, jj, bn, ldb);

                        for (size_type ii = 0; ii < M; ii += BM) {
                            const size_type iend = ((ii + BM) < M) ? (ii + BM) : M;

                            size_type i = ii;
                            for (; i + MR <= iend; i += MR) {
                                const size_type n_full = bn / NR;
                                for (size_type p = 0; p < n_full; ++p) {
                                    gemm_micro_kernel(
                                        a + kk, packB + p * bk * NR,
                                        c + jj + p * NR,
                                        i, bk, lda, ldc);
                                }
                                const size_type j_rem_start = jj + n_full * NR;
                                for (size_type j = j_rem_start; j < jj + bn; ++j) {
                                    T sum = c[i*ldc+j];      T sum1 = c[(i+1)*ldc+j];
                                    T sum2 = c[(i+2)*ldc+j]; T sum3 = c[(i+3)*ldc+j];
                                    for (size_type k = 0; k < bk; ++k) {
                                        T bkj = b[(kk+k)*ldb+j];
                                        sum  += a[ i     *lda+kk+k] * bkj;
                                        sum1 += a[(i+1)*lda+kk+k] * bkj;
                                        sum2 += a[(i+2)*lda+kk+k] * bkj;
                                        sum3 += a[(i+3)*lda+kk+k] * bkj;
                                    }
                                    c[i*ldc+j] = sum;      c[(i+1)*ldc+j] = sum1;
                                    c[(i+2)*ldc+j] = sum2; c[(i+3)*ldc+j] = sum3;
                                }
                            }
                            for (; i < iend; ++i) {
                                for (size_type k = 0; k < bk; ++k) {
                                    const auto a_ik = PT::set1(a[i*lda+kk+k]);
                                    T a_val = a[i*lda+kk+k];
                                    const T* b_row = &b[(kk+k)*ldb+jj];
                                    T* c_row = &c[i*ldc+jj];
                                    size_type j = 0;
                                    const size_type j_vec_end = (bn / PS) * PS;
                                    for (; j < j_vec_end; j += PS) {
                                        auto cv = PT::load(&c_row[j]);
                                        PT::store(&c_row[j], PT::fmadd(a_ik, PT::load(&b_row[j]), cv));
                                    }
                                    for (; j < bn; ++j)
                                        c_row[j] += a_val * b_row[j];
                                }
                            }
                        }
                    }
                }
            } else {
                for (size_type jj = 0; jj < N; jj += BN) {
                    const size_type jend = ((jj + BN) < N) ? (jj + BN) : N;
                    for (size_type kk = 0; kk < K; kk += BK) {
                        const size_type kend = ((kk + BK) < K) ? (kk + BK) : K;
                        for (size_type ii = 0; ii < M; ii += BM) {
                            const size_type iend = ((ii + BM) < M) ? (ii + BM) : M;
                            for (size_type i = ii; i < iend; ++i)
                                for (size_type k = kk; k < kend; ++k) {
                                    T a_val = a[i*lda+k];
                                    for (size_type j = jj; j < jend; ++j)
                                        c[i*ldc+j] += a_val * b[k*ldb+j];
                                }
                        }
                    }
                }
            }
        }

        // ================================================================
        // パッキング判定: M*N*K 単独では不十分
        //
        // パッキングコスト = O(K*N) (B パネルのコピー)
        // 演算量          = O(M*N*K)
        // 償却回数        = M / MR (行方向のマイクロカーネル呼び出し回数)
        //
        // → M が小さいとパッキングのコストが演算量に対して大きい
        // → N < 2*NR ではフルパネルが 1 つ以下でパッキングの意味が薄い
        // → ただし B が L2 を超える場合 (K*N > 64K) はキャッシュミスが支配的で
        //   M が小さくてもパッキングが有利
        // ================================================================

        static bool should_pack(size_type M, size_type N, size_type K) {
            const size_type NR = gemm_NR();
            // N が 2*NR 未満: パネルが 1 つ以下、パッキング不要
            if (N < NR * 2) return false;
            // M が十分大きく演算量があればパッキング
            // 2M 未満 (100³=1M) では B が L2 に収まりパッキング不要
            if (M >= 32 && M * N * K >= 2'000'000) return true;
            // B が大きく L2 に収まらない場合、M が小さくてもパッキング
            if (K * N >= 65536) return true;
            return false;
        }

        // ================================================================
        // 並列化判定: M 方向分割のためのヒューリスティック
        //
        // スレッドプールの起動・同期オーバーヘッド ≈ 300-400 us
        // → single-thread 実行時間がこれを超えないと並列化の恩恵なし
        // → M が小さいとスレッドあたりの行数が不足し効率低下
        //
        // ベンチマーク結果 (128 HW threads):
        //   192³ (7M):       ratio 1.33 → PARALLEL 勝ち (クロスオーバー)
        //   150³ (3.4M):     ratio 0.84 → SINGLE 勝ち
        //   32×500×500 (8M): ratio 0.63 → SINGLE (M 不足)
        //   64×500×500 (16M):ratio 1.30 → PARALLEL (M=64 は最低ライン)
        // ================================================================
        static bool should_parallel(size_type M, size_type N, size_type K) {
            // M が 64 未満: スレッドあたりの行数が不足
            if (M < 64) return false;
            // 総演算量 6M FLOP 以上で並列化
            if (M * N * K < size_type(6'000'000)) return false;
            // スレッドあたりの行数 × K が十分か検証
            // 各スレッドの仕事量 = (M/n_threads) × N × K FLOPs
            // スレッド起動コスト ~50us に見合うには ~500K FLOPs/thread 必要
            auto hw = std::thread::hardware_concurrency();
            if (hw <= 1) return false;
            size_type rows_per_thread = M / hw;
            if (rows_per_thread < MR) return false;  // MR 行未満はマイクロカーネル効率が悪い
            return rows_per_thread * N * K >= size_type(500'000);
        }

        // 並列 gemm: M 方向を分割して各スレッドに配分
        // 各スレッドは C の異なる行範囲を担当 → 同期不要
        static void gemm_parallel(const T* a, const T* b, T* c,
                                   size_type M, size_type N, size_type K,
                                   size_type lda, size_type ldb, size_type ldc) {
            auto& pool = threadPool();
            const unsigned hw = std::thread::hardware_concurrency();
            const unsigned n_threads = (hw > 1) ? hw : 1;

            // MR 単位で行分割 (レジスタブロッキングのアライメント)
            const size_type rows_per_thread = ((M / n_threads + MR - 1) / MR) * MR;

            std::vector<std::future<void>> futures;
            futures.reserve(n_threads);

            for (unsigned t = 0; t < n_threads; ++t) {
                const size_type i_start = t * rows_per_thread;
                if (i_start >= M) break;
                const size_type i_end = (i_start + rows_per_thread < M)
                                         ? i_start + rows_per_thread : M;
                const size_type sub_M = i_end - i_start;

                futures.push_back(pool.submit([=]() {
                    gemm_packed(a + i_start * lda, b,
                                c + i_start * ldc,
                                sub_M, N, K, lda, ldb, ldc);
                }));
            }

            for (auto& f : futures)
                f.get();
        }

        static void gemm_blocked(const T* a, const T* b, T* c,
                                  size_type M, size_type N, size_type K,
                                  size_type lda, size_type ldb, size_type ldc) {
            if (!should_pack(M, N, K))
                gemm_direct(a, b, c, M, N, K, lda, ldb, ldc);
            else if (!should_parallel(M, N, K))
                gemm_packed(a, b, c, M, N, K, lda, ldb, ldc);
            else
                gemm_parallel(a, b, c, M, N, K, lda, ldb, ldc);
        }

    public:

        // raw ポインタ gemm (C += A * B, C は事前ゼロ初期化が必要)
        static void gemm(const T* a, const T* b, T* c,
                         size_type M, size_type N, size_type K,
                         size_type lda, size_type ldb, size_type ldc) {
            gemm_blocked(a, b, c, M, N, K, lda, ldb, ldc);
        }

        // スカラー乗算
        template <typename Mat, typename MatResult>
        static void scalarMultiply(const Mat& mat, const T& scalar, MatResult& result) {
            if (mat.rows() != result.rows() || mat.cols() != result.cols()) {
                throw std::invalid_argument("Matrix dimensions mismatch for scalar multiplication");
            }

            for (size_type i = 0; i < mat.rows(); ++i) {
                for (size_type j = 0; j < mat.cols(); ++j) {
                    result(i, j) = mat(i, j) * scalar;
                }
            }
        }

        // 転置 (AVX2 マイクロカーネル + ブロック転置 + 大規模並列化)
        //
        // double: 4×4 AVX2 カーネル (256-bit = 4 doubles)
        // float:  8×8 AVX2 カーネル (256-bit = 8 floats)
        // その他: スカラーブロック転置
        // 大規模行列 (>= 512×512): スレッドプール並列化
        template <typename Mat, typename MatResult>
        static void transpose(const Mat& mat, MatResult& result) {
            if (mat.rows() != result.cols() || mat.cols() != result.rows()) {
                throw std::invalid_argument("Matrix dimensions mismatch for transpose");
            }

            const size_type R = mat.rows();
            const size_type C = mat.cols();

            if constexpr ((std::is_same_v<T, double> || std::is_same_v<T, float>)
                          && has_simd_packet_v<T>) {
                transpose_block_avx2_dispatch(mat.data(), result.data(), R, C);
            } else if constexpr (std::is_arithmetic_v<T>) {
                transpose_block_scalar(mat.data(), result.data(), R, C);
            } else {
                for (size_type i = 0; i < R; ++i)
                    for (size_type j = 0; j < C; ++j)
                        result(j, i) = mat(i, j);
            }
        }

    private:
        // 並列化閾値: 要素数がこれ以上なら並列化
        static constexpr size_type TRANSPOSE_PARALLEL_THRESHOLD = 512 * 512;

        // スカラー版ブロック転置
        static void transpose_block_scalar(const T* src, T* dst,
                                            size_type R, size_type C) {
            constexpr size_type BS = 32;
            for (size_type ii = 0; ii < R; ii += BS) {
                const size_type iend = (ii + BS < R) ? ii + BS : R;
                for (size_type jj = 0; jj < C; jj += BS) {
                    const size_type jend = (jj + BS < C) ? jj + BS : C;
                    for (size_type i = ii; i < iend; ++i)
                        for (size_type j = jj; j < jend; ++j)
                            dst[j * R + i] = src[i * C + j];
                }
            }
        }

#ifdef __AVX2__
        // =====================================================================
        // AVX2 4×4 double in-register transpose
        // =====================================================================
        static void transpose_4x4d_kernel(const double* src, size_type src_stride,
                                            double* dst, size_type dst_stride) {
            auto r0 = _mm256_loadu_pd(src);
            auto r1 = _mm256_loadu_pd(src + src_stride);
            auto r2 = _mm256_loadu_pd(src + 2 * src_stride);
            auto r3 = _mm256_loadu_pd(src + 3 * src_stride);

            auto t0 = _mm256_unpacklo_pd(r0, r1);
            auto t1 = _mm256_unpackhi_pd(r0, r1);
            auto t2 = _mm256_unpacklo_pd(r2, r3);
            auto t3 = _mm256_unpackhi_pd(r2, r3);

            _mm256_storeu_pd(dst,                    _mm256_permute2f128_pd(t0, t2, 0x20));
            _mm256_storeu_pd(dst + dst_stride,       _mm256_permute2f128_pd(t1, t3, 0x20));
            _mm256_storeu_pd(dst + 2 * dst_stride,   _mm256_permute2f128_pd(t0, t2, 0x31));
            _mm256_storeu_pd(dst + 3 * dst_stride,   _mm256_permute2f128_pd(t1, t3, 0x31));
        }

        // double ブロック転置 (1ブロック帯: 行 [ii, iend) × 全列)
        static void transpose_block_avx2_double_strip(
                const double* src, double* dst,
                size_type R, size_type C,
                size_type ii, size_type iend) {
            constexpr size_type BS = 32;
            for (size_type jj = 0; jj < C; jj += BS) {
                const size_type jend = (jj + BS < C) ? jj + BS : C;
                size_type i = ii;
                for (; i + 4 <= iend; i += 4) {
                    size_type j = jj;
                    for (; j + 4 <= jend; j += 4) {
                        transpose_4x4d_kernel(
                            src + i * C + j, C,
                            dst + j * R + i, R);
                    }
                    for (; j < jend; ++j)
                        for (size_type k = i; k < i + 4; ++k)
                            dst[j * R + k] = src[k * C + j];
                }
                for (; i < iend; ++i)
                    for (size_type j = jj; j < jend; ++j)
                        dst[j * R + i] = src[i * C + j];
            }
        }

        // =====================================================================
        // AVX2 8×8 float in-register transpose
        // =====================================================================
        static void transpose_8x8f_kernel(const float* src, size_type src_stride,
                                            float* dst, size_type dst_stride) {
            // 8 行をロード
            __m256 r0 = _mm256_loadu_ps(src);
            __m256 r1 = _mm256_loadu_ps(src + src_stride);
            __m256 r2 = _mm256_loadu_ps(src + 2 * src_stride);
            __m256 r3 = _mm256_loadu_ps(src + 3 * src_stride);
            __m256 r4 = _mm256_loadu_ps(src + 4 * src_stride);
            __m256 r5 = _mm256_loadu_ps(src + 5 * src_stride);
            __m256 r6 = _mm256_loadu_ps(src + 6 * src_stride);
            __m256 r7 = _mm256_loadu_ps(src + 7 * src_stride);

            // Phase 1: 128-bit レーン内 4×4 転置 (unpacklo/hi + shuffle)
            __m256 t0 = _mm256_unpacklo_ps(r0, r1);  // a0b0 a1b1 a4b4 a5b5
            __m256 t1 = _mm256_unpackhi_ps(r0, r1);  // a2b2 a3b3 a6b6 a7b7
            __m256 t2 = _mm256_unpacklo_ps(r2, r3);
            __m256 t3 = _mm256_unpackhi_ps(r2, r3);
            __m256 t4 = _mm256_unpacklo_ps(r4, r5);
            __m256 t5 = _mm256_unpackhi_ps(r4, r5);
            __m256 t6 = _mm256_unpacklo_ps(r6, r7);
            __m256 t7 = _mm256_unpackhi_ps(r6, r7);

            // Phase 2: shuffle_ps で 2×2 ブロック化
            __m256 s0 = _mm256_shuffle_ps(t0, t2, 0x44); // row0-lo, row0-hi (128-bit lanes)
            __m256 s1 = _mm256_shuffle_ps(t0, t2, 0xEE);
            __m256 s2 = _mm256_shuffle_ps(t1, t3, 0x44);
            __m256 s3 = _mm256_shuffle_ps(t1, t3, 0xEE);
            __m256 s4 = _mm256_shuffle_ps(t4, t6, 0x44);
            __m256 s5 = _mm256_shuffle_ps(t4, t6, 0xEE);
            __m256 s6 = _mm256_shuffle_ps(t5, t7, 0x44);
            __m256 s7 = _mm256_shuffle_ps(t5, t7, 0xEE);

            // Phase 3: permute2f128 で 128-bit レーンを入れ替え
            _mm256_storeu_ps(dst,                    _mm256_permute2f128_ps(s0, s4, 0x20));
            _mm256_storeu_ps(dst + dst_stride,       _mm256_permute2f128_ps(s1, s5, 0x20));
            _mm256_storeu_ps(dst + 2 * dst_stride,   _mm256_permute2f128_ps(s2, s6, 0x20));
            _mm256_storeu_ps(dst + 3 * dst_stride,   _mm256_permute2f128_ps(s3, s7, 0x20));
            _mm256_storeu_ps(dst + 4 * dst_stride,   _mm256_permute2f128_ps(s0, s4, 0x31));
            _mm256_storeu_ps(dst + 5 * dst_stride,   _mm256_permute2f128_ps(s1, s5, 0x31));
            _mm256_storeu_ps(dst + 6 * dst_stride,   _mm256_permute2f128_ps(s2, s6, 0x31));
            _mm256_storeu_ps(dst + 7 * dst_stride,   _mm256_permute2f128_ps(s3, s7, 0x31));
        }

        // float ブロック転置 (1ブロック帯: 行 [ii, iend) × 全列)
        static void transpose_block_avx2_float_strip(
                const float* src, float* dst,
                size_type R, size_type C,
                size_type ii, size_type iend) {
            constexpr size_type BS = 32;
            for (size_type jj = 0; jj < C; jj += BS) {
                const size_type jend = (jj + BS < C) ? jj + BS : C;
                size_type i = ii;
                for (; i + 8 <= iend; i += 8) {
                    size_type j = jj;
                    for (; j + 8 <= jend; j += 8) {
                        transpose_8x8f_kernel(
                            src + i * C + j, C,
                            dst + j * R + i, R);
                    }
                    for (; j < jend; ++j)
                        for (size_type k = i; k < i + 8; ++k)
                            dst[j * R + k] = src[k * C + j];
                }
                for (; i < iend; ++i)
                    for (size_type j = jj; j < jend; ++j)
                        dst[j * R + i] = src[i * C + j];
            }
        }

        // AVX2 ディスパッチ (double/float 共通エントリポイント)
        static void transpose_block_avx2_dispatch(const T* src, T* dst,
                                                    size_type R, size_type C) {
            constexpr size_type BS = 32;
            const size_type total = R * C;

            // 大規模行列: スレッドプール並列化 (行ブロック帯単位)
            if (total >= TRANSPOSE_PARALLEL_THRESHOLD) {
                auto& pool = threadPool();
                const unsigned hw = std::thread::hardware_concurrency();
                const size_type nthreads = (hw > 1) ? hw : 1;
                if (nthreads > 1) {
                    const size_type n_strips = (R + BS - 1) / BS;
                    const size_type strips_per_thread = (n_strips + nthreads - 1) / nthreads;

                    std::vector<std::future<void>> futures;
                    futures.reserve(nthreads);

                    for (size_type t = 0; t < nthreads; ++t) {
                        const size_type strip_begin = t * strips_per_thread;
                        if (strip_begin >= n_strips) break;
                        const size_type strip_end = std::min((t + 1) * strips_per_thread, n_strips);
                        const size_type ii = strip_begin * BS;
                        const size_type iend = std::min(strip_end * BS, R);

                        futures.push_back(pool.submit([=]() {
                            if constexpr (std::is_same_v<T, double>) {
                                transpose_block_avx2_double_strip(src, dst, R, C, ii, iend);
                            } else {
                                transpose_block_avx2_float_strip(src, dst, R, C, ii, iend);
                            }
                        }));
                    }
                    for (auto& f : futures) f.get();
                    return;
                }
            }

            // シングルスレッド
            if constexpr (std::is_same_v<T, double>) {
                transpose_block_avx2_double_strip(src, dst, R, C, 0, R);
            } else {
                transpose_block_avx2_float_strip(src, dst, R, C, 0, R);
            }
        }

#else
        // AVX2 未対応の場合はスカラーブロック転置にフォールバック
        static void transpose_block_avx2_dispatch(const T* src, T* dst,
                                                    size_type R, size_type C) {
            transpose_block_scalar(src, dst, R, C);
        }
#endif

    public:
        // 行列式
        template <typename Mat>
        static T determinant(const Mat& mat) {
            if (mat.rows() != mat.cols()) {
                throw std::invalid_argument("Matrix must be square to calculate determinant");
            }

            const size_type n = mat.rows();

            // 2x2行列の場合の高速パス
            if (n == 2) {
                return mat(0, 0) * mat(1, 1) - mat(0, 1) * mat(1, 0);
            }

            // 3x3行列の場合の高速パス
            if (n == 3) {
                return mat(0, 0) * (mat(1, 1) * mat(2, 2) - mat(1, 2) * mat(2, 1)) -
                    mat(0, 1) * (mat(1, 0) * mat(2, 2) - mat(1, 2) * mat(2, 0)) +
                    mat(0, 2) * (mat(1, 0) * mat(2, 1) - mat(1, 1) * mat(2, 0));
            }

            // より大きな行列の場合は余因子展開を使用
            T det = static_cast<T>(0);
            int sign = 1;

            for (size_type j = 0; j < n; ++j) {
                // 余因子行列を作成
                Mat minor_matrix(n - 1, n - 1);
                for (size_type i = 1; i < n; ++i) {
                    size_type col_idx = 0;
                    for (size_type k = 0; k < n; ++k) {
                        if (k != j) {
                            minor_matrix(i - 1, col_idx++) = mat(i, k);
                        }
                    }
                }

                det += sign * mat(0, j) * determinant(minor_matrix);
                sign = -sign;
            }

            return det;
        }

        // 逆行列
        template <typename Mat, typename MatResult>
        static bool inverse(const Mat& mat, MatResult& result) {
            if (mat.rows() != mat.cols() || result.rows() != mat.rows() || result.cols() != mat.cols()) {
                throw std::invalid_argument("Matrix must be square for inverse calculation");
            }

            const size_type n = mat.rows();

            // 拡張行列を作成 [A|I]
            Mat augmented(n, 2 * n);

            // 拡張行列を初期化
            for (size_type i = 0; i < n; ++i) {
                for (size_type j = 0; j < n; ++j) {
                    augmented(i, j) = mat(i, j);
                    augmented(i, j + n) = (i == j) ? static_cast<T>(1) : static_cast<T>(0);
                }
            }

            // Gauss-Jordan消去法
            for (size_type i = 0; i < n; ++i) {
                // ピボット要素を探す
                size_type pivot_row = i;
                T max_val = std::abs(augmented(i, i));

                for (size_type k = i + 1; k < n; ++k) {
                    if (std::abs(augmented(k, i)) > max_val) {
                        max_val = std::abs(augmented(k, i));
                        pivot_row = k;
                    }
                }

                // 行列が特異の場合
                if (max_val < std::numeric_limits<T>::epsilon()) {
                    return false;  // 逆行列は存在しない
                }

                // 行を交換
                if (pivot_row != i) {
                    for (size_type j = 0; j < 2 * n; ++j) {
                        std::swap(augmented(i, j), augmented(pivot_row, j));
                    }
                }

                // ピボット要素を1にする
                T pivot = augmented(i, i);
                for (size_type j = 0; j < 2 * n; ++j) {
                    augmented(i, j) /= pivot;
                }

                // 他の行のピボット列成分を0にする
                for (size_type k = 0; k < n; ++k) {
                    if (k != i) {
                        T factor = augmented(k, i);
                        for (size_type j = 0; j < 2 * n; ++j) {
                            augmented(k, j) -= factor * augmented(i, j);
                        }
                    }
                }
            }

            // 逆行列部分を結果にコピー
            for (size_type i = 0; i < n; ++i) {
                for (size_type j = 0; j < n; ++j) {
                    result(i, j) = augmented(i, j + n);
                }
            }

            return true;
        }
    };

    // SIMD計算ポリシー
    template <typename T>
    class SIMDComputePolicy : public DefaultComputePolicy<T> {
    public:
        using value_type = T;
        using size_type = std::size_t;

        // ベクトル加算 (SIMD)
        template <typename VecA, typename VecB, typename VecResult>
        static void vector_add(const VecA& a, const VecB& b, VecResult& result) {
            if (a.size() != b.size() || a.size() != result.size()) {
                throw std::invalid_argument("Vector dimensions mismatch for addition");
            }
            if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double>) {
                computation::simd::add_simd<T>(&a[0], &b[0], &result[0], a.size());
            } else {
                DefaultComputePolicy<T>::vector_add(a, b, result);
            }
        }

        // ベクトル減算 (SIMD)
        template <typename VecA, typename VecB, typename VecResult>
        static void vector_subtract(const VecA& a, const VecB& b, VecResult& result) {
            if (a.size() != b.size() || a.size() != result.size()) {
                throw std::invalid_argument("Vector dimensions mismatch for subtraction");
            }
            if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double>) {
                computation::simd::subtract_simd<T>(&a[0], &b[0], &result[0], a.size());
            } else {
                DefaultComputePolicy<T>::vector_subtract(a, b, result);
            }
        }

        // ベクトル内積 (SIMD)
        template <typename VecA, typename VecB>
        static T dot(const VecA& a, const VecB& b) {
            if (a.size() != b.size()) {
                throw std::invalid_argument("Vector dimensions mismatch for dot product");
            }
            if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double>) {
                return computation::simd::dot_product_simd<T>(&a[0], &b[0], a.size());
            } else {
                return DefaultComputePolicy<T>::dot(a, b);
            }
        }

        // ベクトルノルム (SIMD)
        template <typename Vec>
        static T norm(const Vec& vec) {
            if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double>) {
                T sum = computation::simd::dot_product_simd<T>(&vec[0], &vec[0], vec.size());
                return std::sqrt(sum);
            } else {
                return DefaultComputePolicy<T>::norm(vec);
            }
        }

        // ベクトルスカラー乗算 (SIMD)
        template <typename Vec, typename VecResult>
        static void vector_scalar_multiply(const Vec& vec, const T& scalar, VecResult& result) {
            if (vec.size() != result.size()) {
                throw std::invalid_argument("Vector dimensions mismatch for scalar multiplication");
            }
            if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double>) {
                // コピーしてスケール
                std::memcpy(&result[0], &vec[0], vec.size() * sizeof(T));
                computation::simd::scale_simd<T>(&result[0], scalar, result.size());
            } else {
                DefaultComputePolicy<T>::vector_scalar_multiply(vec, scalar, result);
            }
        }
    };

    // MKL計算ポリシー（将来的に実装）
    template <typename T>
    class MKLComputePolicy : public DefaultComputePolicy<T> {
        // Intel MKLを使った最適化された実装をここに追加
    };

    // ポリシーの自動選択
    template <typename T>
    struct OptimalPolicy {
#ifdef USE_MKL
        using type = MKLComputePolicy<T>;
#elif defined(CALX_HAS_AVX512) || defined(CALX_HAS_AVX2) || defined(CALX_HAS_AVX) || \
      defined(CALX_HAS_SSE2) || defined(CALX_HAS_NEON) || defined(USE_SIMD)
        using type = SIMDComputePolicy<T>;
#else
        using type = DefaultComputePolicy<T>;
#endif
    };

    template <typename T>
    using OptimalPolicyType = typename OptimalPolicy<T>::type;

} // namespace calx