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

// eigenvalues.hpp
//
// 固有値問題アルゴリズム
//
// このファイルでは、行列の固有値と固有ベクトルを計算するための
// アルゴリズムを実装します。
//
// 主な機能:
// - 固有値/固有ベクトルの計算
// - 対称行列の固有値分解
// - 一般化固有値問題
// - べき乗法や逆べき乗法などの特殊なアルゴリズム

#ifndef CALX_EIGENVALUES_HPP
#define CALX_EIGENVALUES_HPP

#include <math/core/matrix.hpp>
#include <math/core/vector.hpp>
#include <math/core/Complex.hpp>
#include <math/linalg/decomposition.hpp>

#include <algorithm>
#include <complex>
#include <limits>
#include <random>
#include <utility>
#include <vector>

// 前方宣言 (eigenvalues.hpp 内の相互参照)
namespace calx { namespace algorithms {
    template<typename T>
    std::pair<Vector<T>, Matrix<T>> eigen_symmetric(const Matrix<T>& a);
    template<typename T>
    std::pair<Vector<T>, Matrix<T>>
    eigen_generalized_symmetric(const Matrix<T>& a, const Matrix<T>& b);
}}

namespace calx {
    namespace algorithms {

        //-----------------------------------------------------------------------------
        // 固有値計算の一般的な関数
        //-----------------------------------------------------------------------------

        // 行列の固有値と固有ベクトルを計算
        // 固有値は複素数として返されます（実対称行列でも）
        template<typename T>
        std::pair<Vector<Complex<T>>, Matrix<Complex<T>>>
            eigen(const Matrix<T>& a) {
            const auto n = a.rows();

            if (n != a.cols()) {
                throw DimensionError("Eigenvalue calculation requires a square matrix");
            }

            // 対称性の確認 - 関数名を統一（isSymmetric → is_symmetric）
            bool is_symmetric = true;
            for (std::size_t i = 0; i < n && is_symmetric; ++i) {
                for (std::size_t j = i + 1; j < n && is_symmetric; ++j) {
                    if (!approximately_equal(a(i, j), a(j, i))) {
                        is_symmetric = false;
                    }
                }
            }

            // 実対称行列の場合は対称行列専用の関数を使用
            if (is_symmetric) {
                // 実対称行列の固有値は実数なので、変換が必要
                auto [eigenvalues_real, eigenvectors_real] = eigen_symmetric(a);

                // 実数からの変換
                Vector<Complex<T>> eigenvalues(n);
                for (std::size_t i = 0; i < n; ++i) {
                    eigenvalues[i] = Complex<T>(eigenvalues_real[i], T(0));
                }

                // 実数行列から複素数行列への変換
                Matrix<Complex<T>> eigenvectors(n, n);
                for (std::size_t i = 0; i < n; ++i) {
                    for (std::size_t j = 0; j < n; ++j) {
                        eigenvectors(i, j) = Complex<T>(eigenvectors_real(i, j), T(0));
                    }
                }

                return { eigenvalues, eigenvectors };
            }

            // 非対称行列の場合

            // 複素数での計算が必要
            Vector<Complex<T>> eigenvalues(n);
            Matrix<Complex<T>> eigenvectors(n, n);

#if CALX_HAS_MKL
            if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double>) {
                Matrix<T> Acopy = a;
                std::vector<T> wr(n), wi(n);
                Matrix<T> vr(n, n);  // 右固有ベクトル (実数形式)
                lapack_int info;
                if constexpr (std::is_same_v<T, float>)
                    info = LAPACKE_sgeev(LAPACK_ROW_MAJOR, 'N', 'V',
                                         (lapack_int)n, Acopy.data(), (lapack_int)n,
                                         wr.data(), wi.data(),
                                         nullptr, (lapack_int)n,
                                         vr.data(), (lapack_int)n);
                else
                    info = LAPACKE_dgeev(LAPACK_ROW_MAJOR, 'N', 'V',
                                         (lapack_int)n, Acopy.data(), (lapack_int)n,
                                         wr.data(), wi.data(),
                                         nullptr, (lapack_int)n,
                                         vr.data(), (lapack_int)n);
                if (info > 0)
                    throw MathError("eigen: LAPACKE_geev did not converge");

                // LAPACKE の出力を複素数形式に変換
                for (std::size_t i = 0; i < n; ++i)
                    eigenvalues[i] = Complex<T>(wr[i], wi[i]);

                // 固有ベクトル変換: 複素共役ペアは隣接列に実部/虚部で格納される
                for (std::size_t j = 0; j < n; ++j) {
                    if (wi[j] == T(0)) {
                        // 実固有値 → 実固有ベクトル
                        for (std::size_t i = 0; i < n; ++i)
                            eigenvectors(i, j) = Complex<T>(vr(i, j), T(0));
                    } else if (j + 1 < n && wi[j] > T(0)) {
                        // 複素共役ペア: v(:,j) ± i*v(:,j+1)
                        for (std::size_t i = 0; i < n; ++i) {
                            eigenvectors(i, j)     = Complex<T>(vr(i, j),  vr(i, j + 1));
                            eigenvectors(i, j + 1) = Complex<T>(vr(i, j), -vr(i, j + 1));
                        }
                        ++j;  // 次の列をスキップ
                    }
                }
                return { eigenvalues, eigenvectors };
            }
#endif

            // 実 Schur 分解ベース: A = Q T Q^T
            // T は準上三角 (1×1 実固有値 + 2×2 複素共役ペアブロック)
            auto [T_sch, Q_sch] = schur_decomposition(a);
            const T eps = numeric_traits<T>::epsilon();

            // 準上三角行列 T から固有値を抽出
            std::size_t idx = 0;
            while (idx < n) {
                if (idx + 1 < n && std::abs(T_sch(idx + 1, idx)) > eps *
                    (std::abs(T_sch(idx, idx)) + std::abs(T_sch(idx + 1, idx + 1)) + T(1))) {
                    // 2×2 ブロック → 複素共役ペア
                    T a11 = T_sch(idx, idx), a12 = T_sch(idx, idx + 1);
                    T a21 = T_sch(idx + 1, idx), a22 = T_sch(idx + 1, idx + 1);
                    T tr = a11 + a22;
                    T det = a11 * a22 - a12 * a21;
                    T disc = tr * tr - T(4) * det;
                    if (disc < T(0)) {
                        T re = tr / T(2);
                        T im = std::sqrt(-disc) / T(2);
                        eigenvalues[idx]     = Complex<T>(re,  im);
                        eigenvalues[idx + 1] = Complex<T>(re, -im);
                    } else {
                        T sqd = std::sqrt(disc);
                        eigenvalues[idx]     = Complex<T>((tr + sqd) / T(2), T(0));
                        eigenvalues[idx + 1] = Complex<T>((tr - sqd) / T(2), T(0));
                    }
                    idx += 2;
                } else {
                    // 1×1 ブロック → 実固有値
                    eigenvalues[idx] = Complex<T>(T_sch(idx, idx), T(0));
                    idx += 1;
                }
            }

            // 固有ベクトル: 逆反復法 (inverse iteration)
            // (A - λI)x = b を繰り返し解いて λ の固有ベクトルを求める
            for (std::size_t k = 0; k < n; ++k) {
                Complex<T> lambda = eigenvalues[k];

                // 複素共役ペアの 2 番目は共役を設定
                if (k > 0 && std::abs(lambda.imag()) > eps &&
                    calx::abs(lambda - calx::conj(eigenvalues[k - 1])) < eps * T(100)) {
                    for (std::size_t i = 0; i < n; ++i)
                        eigenvectors(i, k) = calx::conj(eigenvectors(i, k - 1));
                    continue;
                }

                // (A - λI) を複素行列として構築
                Matrix<Complex<T>> B(n, n);
                for (std::size_t i = 0; i < n; ++i)
                    for (std::size_t j = 0; j < n; ++j)
                        B(i, j) = Complex<T>(a(i, j), T(0))
                                  - (i == j ? lambda : Complex<T>(T(0)));

                // 微小摂動でシフト (特異を回避)
                for (std::size_t i = 0; i < n; ++i)
                    B(i, i) += Complex<T>(eps * T(10) * calx::abs(lambda) + eps, T(0));

                // LU 分解 (部分ピボット)
                std::vector<std::size_t> piv(n);
                for (std::size_t i = 0; i < n; ++i) piv[i] = i;
                Matrix<Complex<T>> LU = B;
                for (std::size_t col = 0; col < n; ++col) {
                    // ピボット選択
                    T maxval = T(0);
                    std::size_t maxrow = col;
                    for (std::size_t row = col; row < n; ++row) {
                        T val = calx::abs(LU(row, col));
                        if (val > maxval) { maxval = val; maxrow = row; }
                    }
                    if (maxrow != col) {
                        std::swap(piv[col], piv[maxrow]);
                        for (std::size_t j = 0; j < n; ++j)
                            std::swap(LU(col, j), LU(maxrow, j));
                    }
                    if (calx::abs(LU(col, col)) < eps * eps)
                        LU(col, col) = Complex<T>(eps * eps, T(0));
                    for (std::size_t row = col + 1; row < n; ++row) {
                        LU(row, col) /= LU(col, col);
                        for (std::size_t j = col + 1; j < n; ++j)
                            LU(row, j) -= LU(row, col) * LU(col, j);
                    }
                }

                // 逆反復: x ← solve(B, x) を 3 回繰り返し
                Vector<Complex<T>> x(n);
                for (std::size_t i = 0; i < n; ++i)
                    x[i] = Complex<T>(T(1), T(0));

                for (int iter = 0; iter < 3; ++iter) {
                    // Pb を構築
                    Vector<Complex<T>> b(n);
                    for (std::size_t i = 0; i < n; ++i) b[i] = x[piv[i]];
                    // 前進代入
                    for (std::size_t i = 1; i < n; ++i)
                        for (std::size_t j = 0; j < i; ++j)
                            b[i] -= LU(i, j) * b[j];
                    // 後退代入
                    for (int i = static_cast<int>(n) - 1; i >= 0; --i) {
                        for (std::size_t j = i + 1; j < n; ++j)
                            b[i] -= LU(i, j) * b[j];
                        b[i] /= LU(i, i);
                    }
                    x = b;
                    // 正規化
                    T nrm = T(0);
                    for (std::size_t i = 0; i < n; ++i) nrm += calx::normSq(x[i]);
                    nrm = std::sqrt(nrm);
                    if (nrm > T(0))
                        for (std::size_t i = 0; i < n; ++i) x[i] /= nrm;
                }

                for (std::size_t i = 0; i < n; ++i)
                    eigenvectors(i, k) = x[i];
            }

            return { eigenvalues, eigenvectors };
        }

        // 実対称行列の固有値と固有ベクトルを計算
        // 実対称行列の固有値はすべて実数です
        template<typename T>
        std::pair<Vector<T>, Matrix<T>>
            eigen_symmetric(const Matrix<T>& a) {
            const auto n = a.rows();

            if (n != a.cols()) {
                throw DimensionError("Symmetric eigenvalue calculation requires a square matrix");
            }

            // 対称性の確認
            for (std::size_t i = 0; i < n; ++i) {
                for (std::size_t j = i + 1; j < n; ++j) {
                    if (!approximately_equal(a(i, j), a(j, i))) {
                        throw MathError("Matrix must be symmetric for eigen_symmetric");
                    }
                }
            }

            // 固有値と固有ベクトルを格納する配列
            Vector<T> eigenvalues(n);
            Matrix<T> eigenvectors(n, n);

#if CALX_HAS_MKL
            if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double>) {
                // LAPACKE_syev: 対称行列の全固有値と固有ベクトルを計算
                eigenvectors = a;  // 上書きされる
                lapack_int info;
                if constexpr (std::is_same_v<T, float>)
                    info = LAPACKE_ssyev(LAPACK_ROW_MAJOR, 'V', 'U',
                                         (lapack_int)n, eigenvectors.data(), (lapack_int)n,
                                         eigenvalues.data());
                else
                    info = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'V', 'U',
                                         (lapack_int)n, eigenvectors.data(), (lapack_int)n,
                                         eigenvalues.data());
                if (info > 0)
                    throw MathError("eigen_symmetric: LAPACKE_syev did not converge");
                // LAPACKE_syev は固有値を昇順、固有ベクトルは列として格納 (既存コードと同じ形式)
                return { eigenvalues, eigenvectors };
            }
#endif

            // ================================================================
            // ブロック三重対角化 + 暗黙的 QL 反復 + 直接逆変換
            // (LAPACK DSYEV 方式、Q を陽に形成しない最適化版)
            // ================================================================
            if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double>) {
                // --- Phase 1: ブロック Householder 三重対角化 (コンパクト形式) ---
                auto tri = tridiagonalize_compact(a);
                const std::size_t nm1 = n > 1 ? n - 1 : 0;
                std::vector<T> d_arr(n), e_arr(n, T(0));
                for (std::size_t i = 0; i < n; ++i) d_arr[i] = tri.d[i];
                for (std::size_t i = 0; i < nm1; ++i) e_arr[i] = tri.e[i];

                // --- Phase 2: 暗黙的 QL 反復 (列優先 Z でキャッシュ効率向上) ---
                // Z を列優先で格納: z_cm[j*n + k] = Z(k, j)
                std::vector<T> z_cm(n * n, T(0));
                for (std::size_t i = 0; i < n; ++i)
                    z_cm[i * n + i] = T(1);  // Z = I (列優先)

                const int max_iter = 30 * static_cast<int>(n);
                const T eps = numeric_traits<T>::epsilon();
                int total_iter = 0;

                for (std::size_t l = 0; l < n; ) {
                    std::size_t m = l;
                    while (m + 1 < n) {
                        T tst = std::abs(d_arr[m]) + std::abs(d_arr[m + 1]);
                        if (std::abs(e_arr[m]) <= eps * tst)
                            break;
                        ++m;
                    }

                    if (m == l) { ++l; continue; }

                    if (++total_iter > max_iter)
                        throw MathError("eigen_symmetric: QL iteration did not converge");

                    T g = (d_arr[l + 1] - d_arr[l]) / (T(2) * e_arr[l]);
                    T r = std::sqrt(g * g + T(1));
                    T sign_g = (g >= T(0)) ? T(1) : T(-1);
                    g = d_arr[m] - d_arr[l] + e_arr[l] / (g + sign_g * r);

                    T s = T(1), c = T(1), p = T(0);

                    bool lucky = false;
                    for (std::size_t ii = 0; ii < m - l; ++ii) {
                        std::size_t i = m - 1 - ii;
                        T f = s * e_arr[i];
                        T b = c * e_arr[i];

                        if (std::abs(f) >= std::abs(g)) {
                            c = g / f;
                            r = std::sqrt(c * c + T(1));
                            e_arr[i + 1] = f * r;
                            s = T(1) / r;
                            c *= s;
                        } else {
                            s = f / g;
                            r = std::sqrt(s * s + T(1));
                            e_arr[i + 1] = g * r;
                            c = T(1) / r;
                            s *= c;
                        }

                        if (e_arr[i + 1] == T(0)) {
                            d_arr[i + 1] -= p;
                            e_arr[m] = T(0);
                            lucky = true;
                            break;
                        }

                        g = d_arr[i + 1] - p;
                        r = (d_arr[i] - g) * s + T(2) * c * b;
                        p = s * r;
                        d_arr[i + 1] = g + p;
                        g = c * r - b;

                        // 列優先 Z: 列 i, i+1 は連続メモリ → キャッシュフレンドリー
                        T* __restrict col_i  = z_cm.data() + i * n;
                        T* __restrict col_i1 = z_cm.data() + (i + 1) * n;
                        for (std::size_t k = 0; k < n; ++k) {
                            T zi  = col_i[k];
                            T zi1 = col_i1[k];
                            col_i1[k] = s * zi + c * zi1;
                            col_i[k]  = c * zi - s * zi1;
                        }
                    }

                    if (!lucky) {
                        d_arr[l] -= p;
                        e_arr[l] = g;
                        e_arr[m] = T(0);
                    }
                }

                // --- Phase 3: 直接 Householder 逆変換 ---
                // Z ← Q * Z (Q を形成せず Householder 列を Z に直接適用)
                apply_householder_sequence<T>(
                    tri.refs, n, z_cm.data(), n);

                // 列優先 Z → 行優先 eigenvectors に変換 + 昇順ソート
                std::vector<std::size_t> idx(n);
                for (std::size_t i = 0; i < n; ++i) { eigenvalues[i] = d_arr[i]; idx[i] = i; }
                std::sort(idx.begin(), idx.end(),
                          [&](std::size_t aa, std::size_t bb) { return eigenvalues[aa] < eigenvalues[bb]; });

                Vector<T> sorted_evals(n);
                Matrix<T> sorted_evecs(n, n);
                for (std::size_t i = 0; i < n; ++i) {
                    sorted_evals[i] = eigenvalues[idx[i]];
                    // Z の列 idx[i] (列優先) → sorted_evecs の列 i (行優先)
                    const T* src_col = z_cm.data() + idx[i] * n;
                    for (std::size_t j = 0; j < n; ++j)
                        sorted_evecs(j, i) = src_col[j];
                }
                return { std::move(sorted_evals), std::move(sorted_evecs) };
            }

            // 汎用型: Schur 分解ベースの固有値計算
            auto [Sch, Q_sch] = schur_decomposition(a);
            for (std::size_t i = 0; i < n; ++i) {
                eigenvalues[i] = Sch(i, i);
                for (std::size_t j = 0; j < n; ++j) {
                    eigenvectors(j, i) = Q_sch(j, i);
                }
            }
            // 固有値を昇順にソート
            std::vector<std::size_t> idx(n);
            for (std::size_t i = 0; i < n; ++i) idx[i] = i;
            std::sort(idx.begin(), idx.end(),
                      [&](std::size_t aa, std::size_t bb) {
                          return eigenvalues[aa] < eigenvalues[bb];
                      });
            Vector<T> sorted_evals(n);
            Matrix<T> sorted_evecs(n, n);
            for (std::size_t i = 0; i < n; ++i) {
                sorted_evals[i] = eigenvalues[idx[i]];
                for (std::size_t j = 0; j < n; ++j) {
                    sorted_evecs(j, i) = eigenvectors(j, idx[i]);
                }
            }
            eigenvalues = sorted_evals;
            eigenvectors = sorted_evecs;

            return { eigenvalues, eigenvectors };
        }

        // エルミート行列（複素対称行列）の固有値と固有ベクトルを計算
        // エルミート行列の固有値はすべて実数です
        template<typename T>
        std::pair<Vector<T>, Matrix<Complex<T>>>
            eigen_hermitian(const Matrix<Complex<T>>& a) {
            const auto n = a.rows();

            if (n != a.cols()) {
                throw DimensionError("Hermitian eigenvalue calculation requires a square matrix");
            }

            // エルミート性の確認 (A^H = A)
            {
                T eps = numeric_traits<T>::epsilon() * T(100);
                for (std::size_t i = 0; i < n; ++i) {
                    for (std::size_t j = i + 1; j < n; ++j) {
                        if (calx::abs(a(i, j) - calx::conj(a(j, i))) > eps)
                            throw MathError("Matrix must be Hermitian for eigen_hermitian");
                    }
                }
            }

            // 固有値と固有ベクトルを格納する配列
            Vector<T> eigenvalues(n);
            Matrix<Complex<T>> eigenvectors(n, n);

#if CALX_HAS_MKL
            if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double>) {
                eigenvectors = a;
                lapack_int info;
                if constexpr (std::is_same_v<T, float>)
                    info = LAPACKE_cheev(LAPACK_ROW_MAJOR, 'V', 'U',
                                         (lapack_int)n,
                                         reinterpret_cast<lapack_complex_float*>(eigenvectors.data()),
                                         (lapack_int)n, eigenvalues.data());
                else
                    info = LAPACKE_zheev(LAPACK_ROW_MAJOR, 'V', 'U',
                                         (lapack_int)n,
                                         reinterpret_cast<lapack_complex_double*>(eigenvectors.data()),
                                         (lapack_int)n, eigenvalues.data());
                if (info > 0)
                    throw MathError("eigen_hermitian: LAPACKE_heev did not converge");
                return { eigenvalues, eigenvectors };
            }
#endif
            // エルミート行列 A = X + iY を 2n×2n 実対称行列に変換して解く
            // A(u+iv) = λ(u+iv) ⇔ Xu-Yv = λu, Yu+Xv = λv
            // M = [[X, -Y], [Y, X]] の固有値は A の固有値の 2 重コピー
            // (X は対称, Y は反対称なので M は対称)
            Matrix<T> M(2 * n, 2 * n, T(0));
            for (std::size_t i = 0; i < n; ++i) {
                for (std::size_t j = 0; j < n; ++j) {
                    T re = a(i, j).real();
                    T im = a(i, j).imag();
                    M(i, j) = re;           // X
                    M(i, n + j) = -im;      // -Y
                    M(n + i, j) = im;       // Y
                    M(n + i, n + j) = re;   // X
                }
            }

            auto [evals_2n, evecs_2n] = eigen_symmetric(M);

            // 固有値は 2 重に現れるので、1 つおきに取得 (昇順)
            for (std::size_t i = 0; i < n; ++i)
                eigenvalues[i] = evals_2n[2 * i];

            // 固有ベクトルを Complex に変換
            // M の固有ベクトル [u; v] から A の固有ベクトル u + iv を構成
            for (std::size_t k = 0; k < n; ++k) {
                std::size_t k2 = 2 * k; // 対応する 2n 側のインデックス
                T norm_sq = T(0);
                for (std::size_t i = 0; i < n; ++i) {
                    T re = evecs_2n(i, k2);
                    T im = evecs_2n(n + i, k2);
                    eigenvectors(i, k) = Complex<T>(re, im);
                    norm_sq += re * re + im * im;
                }
                // 正規化
                T inv_norm = T(1) / std::sqrt(norm_sq);
                for (std::size_t i = 0; i < n; ++i)
                    eigenvectors(i, k) = eigenvectors(i, k) * Complex<T>(inv_norm);
            }

            return { eigenvalues, eigenvectors };
        }

        //-----------------------------------------------------------------------------
        // 一般化固有値問題
        //-----------------------------------------------------------------------------

        // 一般化固有値問題 Ax = λBx の解
        // (A,B)は一般化固有対と呼ばれます
        template<typename T>
        std::pair<Vector<Complex<T>>, Matrix<Complex<T>>>
            eigen_generalized(const Matrix<T>& a, const Matrix<T>& b) {
            const auto n = a.rows();

            if (n != a.cols() || b.rows() != n || b.cols() != n) {
                throw DimensionError("Generalized eigenvalue problem requires square matrices of same size");
            }

            // 対称性の確認
            bool is_symmetric_a = true;
            bool is_symmetric_b = true;

            for (std::size_t i = 0; i < n && (is_symmetric_a || is_symmetric_b); ++i) {
                for (std::size_t j = i + 1; j < n && (is_symmetric_a || is_symmetric_b); ++j) {
                    if (is_symmetric_a && !approximately_equal(a(i, j), a(j, i))) {
                        is_symmetric_a = false;
                    }
                    if (is_symmetric_b && !approximately_equal(b(i, j), b(j, i))) {
                        is_symmetric_b = false;
                    }
                }
            }

            // 両方とも対称行列の場合は対称一般化固有値問題
            if (is_symmetric_a && is_symmetric_b) {
                auto [eigenvalues_real, eigenvectors_real] = eigen_generalized_symmetric(a, b);

                // 実数値を複素数に変換
                Vector<Complex<T>> eigenvalues(n);
                Matrix<Complex<T>> eigenvectors(n, n);

                for (std::size_t i = 0; i < n; ++i) {
                    eigenvalues[i] = Complex<T>(eigenvalues_real[i], T(0));
                    for (std::size_t j = 0; j < n; ++j) {
                        eigenvectors(i, j) = Complex<T>(eigenvectors_real(i, j), T(0));
                    }
                }

                return { eigenvalues, eigenvectors };
            }

            // 非対称の一般化固有値問題

            // 固有値と固有ベクトルを格納する配列
            Vector<Complex<T>> eigenvalues(n);
            Matrix<Complex<T>> eigenvectors(n, n);

            // QZ 分解による一般化固有値問題の解法
            // (MKL LAPACKE パスは Matrix::data() 未対応のため将来実装予定)
            auto [S, Tri, Qm, Zm] = qz_decomposition(a, b);

            const T eps_eig = numeric_traits<T>::epsilon();

            // (S, Tri) のブロック対角から一般化固有値を抽出
            std::size_t i = 0;
            while (i < n) {
                if (i + 1 < n && std::abs(S(i + 1, i)) > eps_eig *
                    (std::abs(S(i, i)) + std::abs(S(i + 1, i + 1)) + T(1))) {
                    // 2×2 ブロック → 複素共役対
                    T a11 = S(i, i), a12 = S(i, i + 1);
                    T a21 = S(i + 1, i), a22 = S(i + 1, i + 1);
                    T b11 = Tri(i, i), b12 = Tri(i, i + 1), b22 = Tri(i + 1, i + 1);

                    // 2×2 一般化固有値: det(S22 - λ T22) = 0
                    // (a11-λb11)(a22-λb22) - (a12-λb12)*a21 = 0
                    // 係数: b11*b22*λ² - (a11*b22+a22*b11-a21*b12)*λ + (a11*a22-a12*a21) = 0
                    T qa = b11 * b22;
                    T qb = -(a11 * b22 + a22 * b11 - a21 * b12);
                    T qc = a11 * a22 - a12 * a21;

                    if (std::abs(qa) < eps_eig) {
                        // 退化ケース
                        eigenvalues[i] = Complex<T>(std::numeric_limits<T>::infinity(), T(0));
                        eigenvalues[i + 1] = Complex<T>(std::numeric_limits<T>::infinity(), T(0));
                    } else {
                        T disc = qb * qb - T(4) * qa * qc;
                        T real_part = -qb / (T(2) * qa);
                        if (disc < T(0)) {
                            T imag_part = std::sqrt(-disc) / (T(2) * qa);
                            eigenvalues[i] = Complex<T>(real_part, imag_part);
                            eigenvalues[i + 1] = Complex<T>(real_part, -imag_part);
                        } else {
                            T sq = std::sqrt(disc);
                            eigenvalues[i] = Complex<T>((-qb + sq) / (T(2) * qa), T(0));
                            eigenvalues[i + 1] = Complex<T>((-qb - sq) / (T(2) * qa), T(0));
                        }
                    }
                    i += 2;
                } else {
                    // 1×1 ブロック
                    T alpha = S(i, i);
                    T beta_val = Tri(i, i);
                    if (std::abs(beta_val) < eps_eig) {
                        eigenvalues[i] = Complex<T>(
                            std::numeric_limits<T>::infinity(), T(0));
                    } else {
                        eigenvalues[i] = Complex<T>(alpha / beta_val, T(0));
                    }
                    ++i;
                }
            }

            // 固有ベクトル: 逆変換で近似的に取得
            // (S - λ_i T) z = 0 を解いて Zm * z で元空間に変換
            for (std::size_t col = 0; col < n; ++col) {
                if (std::abs(eigenvalues[col].imag()) > eps_eig) {
                    // 複素固有値: 近似的に Z の対応列を使用
                    for (std::size_t row = 0; row < n; ++row) {
                        eigenvectors(row, col) = Complex<T>(Zm(row, col), T(0));
                    }
                    continue;
                }
                // 実固有値: 逆反復で固有ベクトルを精密化
                T lam = eigenvalues[col].real();
                // (A - λB)v = 0 の近似解を逆反復で求める
                // 初期ベクトル: Z の col 列
                Matrix<T> C(n, n);
                for (std::size_t r = 0; r < n; ++r) {
                    for (std::size_t c = 0; c < n; ++c) {
                        C(r, c) = a(r, c) - lam * b(r, c);
                    }
                }
                // C に微小摂動 (特異を避ける)
                for (std::size_t r = 0; r < n; ++r) {
                    C(r, r) += eps_eig * (std::abs(C(r, r)) + T(1));
                }
                // 逆反復 (最大 5 ステップ)
                Vector<T> v(n);
                for (std::size_t r = 0; r < n; ++r) v[r] = Zm(r, col);
                try {
                    auto [lu_mat, pivots] = lu_decomposition(C);
                    for (int iter = 0; iter < 5; ++iter) {
                        // ピボット置換
                        Vector<T> rhs = v;
                        for (std::size_t ii = 0; ii < n; ++ii)
                            if (pivots[ii] != ii) std::swap(rhs[ii], rhs[pivots[ii]]);
                        // 前進代入
                        for (std::size_t ii = 1; ii < n; ++ii)
                            for (std::size_t jj = 0; jj < ii; ++jj)
                                rhs[ii] -= lu_mat(ii, jj) * rhs[jj];
                        // 後退代入
                        for (std::size_t ii = n; ii-- > 0; ) {
                            for (std::size_t jj = ii + 1; jj < n; ++jj)
                                rhs[ii] -= lu_mat(ii, jj) * rhs[jj];
                            rhs[ii] /= lu_mat(ii, ii);
                        }
                        v = rhs;
                        // 正規化
                        T nrm = T(0);
                        for (std::size_t r = 0; r < n; ++r) nrm += v[r] * v[r];
                        nrm = std::sqrt(nrm);
                        if (nrm > eps_eig)
                            for (std::size_t r = 0; r < n; ++r) v[r] /= nrm;
                    }
                    for (std::size_t r = 0; r < n; ++r)
                        eigenvectors(r, col) = Complex<T>(v[r], T(0));
                } catch (...) {
                    for (std::size_t r = 0; r < n; ++r) {
                        eigenvectors(r, col) = Complex<T>(Zm(r, col), T(0));
                    }
                }
            }

            return { eigenvalues, eigenvectors };
        }

        // 対称一般化固有値問題 Ax = λBx の解
        // (A,B)は対称行列のペアで、Bは正定値
        template<typename T>
        std::pair<Vector<T>, Matrix<T>>
            eigen_generalized_symmetric(const Matrix<T>& a, const Matrix<T>& b) {
            const auto n = a.rows();

            if (n != a.cols() || b.rows() != n || b.cols() != n) {
                throw DimensionError("Symmetric generalized eigenvalue problem requires square matrices of same size");
            }

            // 対称性の確認
            for (std::size_t i = 0; i < n; ++i) {
                for (std::size_t j = i + 1; j < n; ++j) {
                    if (!approximately_equal(a(i, j), a(j, i)) || !approximately_equal(b(i, j), b(j, i))) {
                        throw MathError("Matrices must be symmetric for eigen_generalized_symmetric");
                    }
                }
            }

            // 固有値と固有ベクトルを格納する配列
            Vector<T> eigenvalues(n);
            Matrix<T> eigenvectors(n, n);

            // Cholesky + 標準固有値変換による対称一般化固有値問題の解法
            // (MKL LAPACKE パスは Matrix::data() 未対応のため将来実装予定)

            // 方法: Cholesky分解を使用してL^T * L = B とし、
            // C = L^(-1) * A * L^(-T) の固有値問題に変換

            try {
                // B の Cholesky 分解: B = L * L^T
                Matrix<T> L = cholesky_decomposition(b);

                // C = L^{-1} A L^{-T} を三角ソルブで構築 (逆行列を明示的に作らない)
                // Step 1: Y = L^{-1} A  (前進代入で各列を解く)
                Matrix<T> Y(n, n);
                for (std::size_t col = 0; col < n; ++col) {
                    // L * y_col = A_col を解く
                    for (std::size_t i = 0; i < n; ++i) {
                        T sum = a(i, col);
                        for (std::size_t k = 0; k < i; ++k)
                            sum -= L(i, k) * Y(k, col);
                        Y(i, col) = sum / L(i, i);
                    }
                }

                // Step 2: C = Y * L^{-T}
                // C^T = L^{-1} Y^T なので、C^T の各列 j (= C の各行 j) について
                // L * z_j = Y^T の j 列 (= Y の j 行) を前進代入で解く
                // → z_j = C^T の j 列 = C の j 行
                Matrix<T> C(n, n);
                for (std::size_t row = 0; row < n; ++row) {
                    // L * z = Y(row, :)^T を前進代入で解く → z = C(row, :)^T
                    for (std::size_t i = 0; i < n; ++i) {
                        T sum = Y(row, i);
                        for (std::size_t k = 0; k < i; ++k)
                            sum -= L(i, k) * C(row, k);
                        C(row, i) = sum / L(i, i);
                    }
                }

                // C の固有値と固有ベクトルを計算
                auto [c_eigenvalues, c_eigenvectors] = eigen_symmetric(C);

                eigenvalues = c_eigenvalues;

                // 固有ベクトルを変換: v = L^{-T} * c_v
                // L^T * v = c_v を後退代入で解く
                for (std::size_t j = 0; j < n; ++j) {
                    for (std::size_t ii = 0; ii < n; ++ii) {
                        std::size_t i = n - 1 - ii;
                        T sum = c_eigenvectors(i, j);
                        for (std::size_t k = i + 1; k < n; ++k)
                            sum -= L(k, i) * eigenvectors(k, j);
                        eigenvectors(i, j) = sum / L(i, i);
                    }
                }
            }
            catch (const std::exception&) {
                throw MathError("Symmetric generalized eigenvalue calculation failed: B matrix is not positive definite");
            }

            return { eigenvalues, eigenvectors };
        }

        //-----------------------------------------------------------------------------
        // 特殊な固有値アルゴリズム
        //-----------------------------------------------------------------------------

        // べき乗法（最大絶対値の固有値とその固有ベクトルを計算）
        template<typename T>
        std::pair<T, Vector<T>> power_method(const Matrix<T>& a,
            T tolerance = T(1e-10),
            std::size_t max_iterations = 1000) {
            const auto n = a.rows();

            if (n != a.cols()) {
                throw DimensionError("Power method requires a square matrix");
            }

            // 初期ベクトルはランダムまたは単位ベクトル
            Vector<T> x(n, T(1));
            x[0] = T(1);  // 先頭要素を1に設定

            // 初期ベクトルの正規化
            const T x_norm = std::sqrt(dot(x, x));
            for (std::size_t i = 0; i < n; ++i) {
                x[i] /= x_norm;
            }

            // 固有値の初期推定値
            T lambda = T(0);
            T lambda_prev;

            // べき乗法の反復
            for (std::size_t iter = 0; iter < max_iterations; ++iter) {
                // 前回の固有値を保存
                lambda_prev = lambda;

                // y = A * x
                Vector<T> y(n, T(0));
                for (std::size_t i = 0; i < n; ++i) {
                    for (std::size_t j = 0; j < n; ++j) {
                        y[i] += a(i, j) * x[j];
                    }
                }

                // レイリー商 λ = (x^T * A * x) / (x^T * x)
                lambda = dot(x, y);

                // 新しいベクトルの正規化
                const T y_norm = std::sqrt(dot(y, y));
                for (std::size_t i = 0; i < n; ++i) {
                    x[i] = y[i] / y_norm;
                }

                // 収束判定
                if (std::abs(lambda - lambda_prev) < tolerance * std::abs(lambda)) {
                    break;
                }

                // 最大反復回数の確認
                if (iter == max_iterations - 1) {
                    throw std::runtime_error("Power method did not converge within "
                        + std::to_string(max_iterations) + " iterations");
                }
            }

            return { lambda, x };
        }

        // 逆べき乗法（特定の値に最も近い固有値とその固有ベクトルを計算）
        template<typename T>
        std::pair<T, Vector<T>> inverse_power_method(const Matrix<T>& a,
            T sigma = T(0),  // シフト（目標値）
            T tolerance = T(1e-10),
            std::size_t max_iterations = 1000) {
            const auto n = a.rows();

            if (n != a.cols()) {
                throw DimensionError("Inverse power method requires a square matrix");
            }

            // (A - σI)の計算
            Matrix<T> a_shifted = a;
            for (std::size_t i = 0; i < n; ++i) {
                a_shifted(i, i) -= sigma;
            }

            // (A - σI)のLU分解
            auto [lu, pivots] = lu_decomposition(a_shifted);

            // 初期ベクトル: 全成分が非零 (対角行列でも全固有空間にアクセス可能)
            Vector<T> x(n, T(1));

            // 初期ベクトルの正規化
            {
                const T x_norm = std::sqrt(dot(x, x));
                for (std::size_t i = 0; i < n; ++i) {
                    x[i] /= x_norm;
                }
            }

            // 固有値の初期推定値
            T mu = T(0);
            T mu_prev;

            // 逆べき乗法の反復
            for (std::size_t iter = 0; iter < max_iterations; ++iter) {
                // 前回の固有値を保存
                mu_prev = mu;

                // (A - σI)y = x を解く（LU分解済みの行列を使用）
                Vector<T> y = x;
                // ピボット置換
                for (typename Matrix<T>::size_type pi = 0; pi < n; ++pi) {
                    if (pivots[pi] != pi) {
                        std::swap(y[pi], y[pivots[pi]]);
                    }
                }
                // 前進代入
                for (typename Matrix<T>::size_type fi = 1; fi < n; ++fi) {
                    for (typename Matrix<T>::size_type fj = 0; fj < fi; ++fj) {
                        y[fi] -= lu(fi, fj) * y[fj];
                    }
                }
                // 後退代入
                for (typename Matrix<T>::size_type bi = n; bi-- > 0; ) {
                    for (typename Matrix<T>::size_type bj = bi + 1; bj < n; ++bj) {
                        y[bi] -= lu(bi, bj) * y[bj];
                    }
                    y[bi] /= lu(bi, bi);
                }

                // (A-σI)^{-1} のレイリー商 μ = x^T * y (x は正規化済み)
                mu = dot(x, y);

                // 新しいベクトルの正規化
                const T y_norm = std::sqrt(dot(y, y));
                for (std::size_t i = 0; i < n; ++i) {
                    x[i] = y[i] / y_norm;
                }

                // 収束判定
                if (std::abs(mu - mu_prev) < tolerance * std::abs(mu)) {
                    break;
                }

                // 最大反復回数の確認
                if (iter == max_iterations - 1) {
                    throw std::runtime_error("Inverse power method did not converge within "
                        + std::to_string(max_iterations) + " iterations");
                }
            }

            // 実際の固有値を計算（シフトを戻す）
            const T lambda = sigma + T(1) / mu;

            return { lambda, x };
        }

        // レイリー商反復法（高速な収束を持つ固有値計算）
        template<typename T>
        std::pair<T, Vector<T>> rayleigh_quotient_iteration(const Matrix<T>& a,
            const Vector<T>& initial_guess,
            T tolerance = T(1e-10),
            std::size_t max_iterations = 100) {
            const auto n = a.rows();

            if (n != a.cols()) {
                throw DimensionError("Rayleigh quotient iteration requires a square matrix");
            }

            if (initial_guess.size() != n) {
                throw DimensionError("Initial guess vector size must match matrix dimension");
            }

            // 初期ベクトルのコピーと正規化
            Vector<T> x = initial_guess;
            const T x_norm = std::sqrt(dot(x, x));
            for (std::size_t i = 0; i < n; ++i) {
                x[i] /= x_norm;
            }

            // 初期レイリー商（固有値の推定値）
            T rho = T(0);
            for (std::size_t i = 0; i < n; ++i) {
                for (std::size_t j = 0; j < n; ++j) {
                    rho += x[i] * a(i, j) * x[j];
                }
            }

            // レイリー商反復
            for (std::size_t iter = 0; iter < max_iterations; ++iter) {
                // (A - ρI)の計算
                Matrix<T> a_shifted = a;
                for (std::size_t i = 0; i < n; ++i) {
                    a_shifted(i, i) -= rho;
                }

                // (A - ρI)のLU分解
                Matrix<T> lu;
                std::vector<typename Matrix<T>::size_type> pivots;

                try {
                    auto decomp = lu_decomposition(a_shifted);
                    lu = std::move(decomp.first);
                    pivots = std::move(decomp.second);
                }
                catch (const std::exception&) {
                    // LU分解が失敗した場合（特異行列）
                    // ρが固有値に非常に近い場合に発生する可能性がある
                    // 微小な摂動を加えて再試行
                    const T epsilon = numeric_traits<T>::epsilon() * T(100);
                    rho += epsilon;
                    continue;
                }

                // (A - ρI)y = x を解く（LU分解済み行列を使用）
                Vector<T> y = x;
                for (typename Matrix<T>::size_type pi = 0; pi < n; ++pi) {
                    if (pivots[pi] != pi) {
                        std::swap(y[pi], y[pivots[pi]]);
                    }
                }
                for (typename Matrix<T>::size_type fi = 1; fi < n; ++fi) {
                    for (typename Matrix<T>::size_type fj = 0; fj < fi; ++fj) {
                        y[fi] -= lu(fi, fj) * y[fj];
                    }
                }
                for (typename Matrix<T>::size_type bi = n; bi-- > 0; ) {
                    for (typename Matrix<T>::size_type bj = bi + 1; bj < n; ++bj) {
                        y[bi] -= lu(bi, bj) * y[bj];
                    }
                    y[bi] /= lu(bi, bi);
                }

                // 新しいベクトルの正規化
                const T y_norm = std::sqrt(dot(y, y));
                for (std::size_t i = 0; i < n; ++i) {
                    y[i] /= y_norm;
                }

                // 新しいレイリー商の計算
                T rho_new = T(0);
                for (std::size_t i = 0; i < n; ++i) {
                    for (std::size_t j = 0; j < n; ++j) {
                        rho_new += y[i] * a(i, j) * y[j];
                    }
                }

                // 収束判定
                if (std::abs(rho_new - rho) < tolerance * std::abs(rho_new)) {
                    rho = rho_new;
                    x = y;
                    break;
                }

                // 更新
                rho = rho_new;
                x = y;

                // 最大反復回数の確認
                if (iter == max_iterations - 1) {
                    throw std::runtime_error("Rayleigh quotient iteration did not converge within "
                        + std::to_string(max_iterations) + " iterations");
                }
            }

            return { rho, x };
        }

        //-----------------------------------------------------------------------------
        // ゲルシュゴリン円定理 (Gershgorin Circle Theorem)
        //-----------------------------------------------------------------------------
        // 各固有値は少なくとも 1 つのゲルシュゴリン円 D(a_ii, R_i) に含まれる。
        // D(a_ii, R_i) = { z ∈ C : |z - a_ii| ≤ R_i }
        // R_i = Σ_{j≠i} |a_ij|

        // ゲルシュゴリン円: 中心と半径のペア
        template<typename T>
        struct GershgorinDisc {
            T center;   // a_ii (対角要素)
            T radius;   // Σ_{j≠i} |a_ij|
        };

        // 全ゲルシュゴリン円を返す
        template<typename T>
        std::vector<GershgorinDisc<T>> gershgorin_discs(const Matrix<T>& A) {
            const auto n = A.rows();
            if (n != A.cols())
                throw DimensionError("gershgorin_discs: matrix must be square");

            std::vector<GershgorinDisc<T>> discs(n);
            for (std::size_t i = 0; i < n; ++i) {
                discs[i].center = A(i, i);
                T r = T(0);
                for (std::size_t j = 0; j < n; ++j)
                    if (j != i) r += std::abs(A(i, j));
                discs[i].radius = r;
            }
            return discs;
        }

        // 固有値の存在範囲 [min, max] を返す (行ベースと列ベースの共通部分)
        template<typename T>
        std::pair<T, T> gershgorin_bounds(const Matrix<T>& A) {
            const auto n = A.rows();
            if (n != A.cols())
                throw DimensionError("gershgorin_bounds: matrix must be square");

            // 行ベースの円 (初期値は最初の行で設定)
            T r0 = T(0);
            for (std::size_t j = 1; j < n; ++j) r0 += std::abs(A(0, j));
            T lo = A(0, 0) - r0;
            T hi = A(0, 0) + r0;
            for (std::size_t i = 0; i < n; ++i) {
                T r = T(0);
                for (std::size_t j = 0; j < n; ++j)
                    if (j != i) r += std::abs(A(i, j));
                lo = std::min(lo, A(i, i) - r);
                hi = std::max(hi, A(i, i) + r);
            }

            // 列ベースの円 (A^T のゲルシュゴリン円)
            for (std::size_t j = 0; j < n; ++j) {
                T r = T(0);
                for (std::size_t i = 0; i < n; ++i)
                    if (i != j) r += std::abs(A(i, j));
                lo = std::max(lo, A(j, j) - r);  // 共通部分 → max of lower bounds
                hi = std::min(hi, A(j, j) + r);  // 共通部分 → min of upper bounds
            }

            // 共通部分が空の場合は行ベースに戻す
            if (lo > hi) {
                lo = A(0, 0) - r0;
                hi = A(0, 0) + r0;
                for (std::size_t i = 0; i < n; ++i) {
                    T r = T(0);
                    for (std::size_t j = 0; j < n; ++j)
                        if (j != i) r += std::abs(A(i, j));
                    lo = std::min(lo, A(i, i) - r);
                    hi = std::max(hi, A(i, i) + r);
                }
            }

            return {lo, hi};
        }

        //-----------------------------------------------------------------------------
        // 複素行列の固有値ソルバー (ComplexEigenSolver)
        //-----------------------------------------------------------------------------
        // 一般複素行列 A ∈ C^{n×n} の固有値・固有ベクトルを計算。
        // 方法: 複素 Hessenberg 分解 → 複素 QR 反復 (単一シフト)

        template<typename T>
        std::pair<Vector<Complex<T>>, Matrix<Complex<T>>>
        eigen_complex(const Matrix<Complex<T>>& A) {
            using C = Complex<T>;
            const auto n = A.rows();
            if (n != A.cols())
                throw DimensionError("eigen_complex: matrix must be square");

            if (n == 0) return {{}, {}};

            if (n == 1) {
                Vector<C> evals(1);
                evals[0] = A(0, 0);
                Matrix<C> evecs(1, 1, C(T(1), T(0)));
                return {evals, evecs};
            }

            // Hermitian チェック: A = A* なら実固有値 (高速パス)
            bool is_hermitian = true;
            for (std::size_t i = 0; i < n && is_hermitian; ++i)
                for (std::size_t j = i + 1; j < n && is_hermitian; ++j)
                    if (calx::abs(A(i, j) - calx::conj(A(j, i))) >
                        numeric_traits<T>::epsilon() * T(100))
                        is_hermitian = false;

            // Step 1: 複素 Hessenberg 分解 A = Q H Q*
            // Householder 変換で上 Hessenberg 形式に変換
            Matrix<C> H = A;
            Matrix<C> Q(n, n, C(T(0)));
            for (std::size_t i = 0; i < n; ++i) Q(i, i) = C(T(1));

            for (std::size_t k = 0; k + 2 < n; ++k) {
                // v = H(k+1:n, k) の Householder ベクトルを計算
                std::size_t m = n - k - 1;
                std::vector<C> v(m);
                for (std::size_t i = 0; i < m; ++i)
                    v[i] = H(k + 1 + i, k);

                T nrm = T(0);
                for (auto& vi : v) nrm += calx::normSq(vi);
                nrm = std::sqrt(nrm);

                if (nrm < numeric_traits<T>::epsilon()) continue;

                // v[0] にフェーズを合わせる
                T abs_v0 = calx::abs(v[0]);
                C alpha = (abs_v0 > T(0)) ?
                    C(-nrm) * v[0] / C(abs_v0) : C(nrm, T(0));
                v[0] -= alpha;

                // v を正規化
                T vnrm = T(0);
                for (auto& vi : v) vnrm += calx::normSq(vi);
                vnrm = std::sqrt(vnrm);
                if (vnrm < numeric_traits<T>::epsilon()) continue;
                for (auto& vi : v) vi /= vnrm;

                // H ← (I - 2vv*) H (I - 2vv*)
                // 左から: H(k+1:n, :) -= 2 v (v* H(k+1:n, :))
                for (std::size_t j = k; j < n; ++j) {
                    C dot(T(0));
                    for (std::size_t i = 0; i < m; ++i)
                        dot += calx::conj(v[i]) * H(k + 1 + i, j);
                    dot *= C(T(2));
                    for (std::size_t i = 0; i < m; ++i)
                        H(k + 1 + i, j) -= v[i] * dot;
                }
                // 右から: H(:, k+1:n) -= 2 (H(:, k+1:n) v) v*
                for (std::size_t i = 0; i < n; ++i) {
                    C dot(T(0));
                    for (std::size_t j = 0; j < m; ++j)
                        dot += H(i, k + 1 + j) * v[j];
                    dot *= C(T(2));
                    for (std::size_t j = 0; j < m; ++j)
                        H(i, k + 1 + j) -= dot * calx::conj(v[j]);
                }
                // Q 蓄積: Q(:, k+1:n) -= 2 (Q(:, k+1:n) v) v*
                for (std::size_t i = 0; i < n; ++i) {
                    C dot(T(0));
                    for (std::size_t j = 0; j < m; ++j)
                        dot += Q(i, k + 1 + j) * v[j];
                    dot *= C(T(2));
                    for (std::size_t j = 0; j < m; ++j)
                        Q(i, k + 1 + j) -= dot * calx::conj(v[j]);
                }
            }

            // サブ対角要素のクリーンアップ
            for (std::size_t i = 0; i + 2 < n; ++i)
                for (std::size_t j = 0; j < i; ++j)
                    H(i + 2, j) = C(T(0));

            // Step 2: 複素 QR 反復 (単一シフト + Givens 回転)
            const int max_iter = 300 * static_cast<int>(n);
            const T eps = numeric_traits<T>::epsilon();
            int ihi = static_cast<int>(n);
            int total_iter = 0;

            while (ihi > 1 && total_iter < max_iter) {
                // デフレーション
                int ilo = ihi - 1;
                while (ilo > 0) {
                    T threshold = eps * (calx::abs(H(ilo - 1, ilo - 1)) +
                                         calx::abs(H(ilo, ilo)));
                    if (threshold == T(0)) threshold = eps;
                    if (calx::abs(H(ilo, ilo - 1)) <= threshold) {
                        H(ilo, ilo - 1) = C(T(0));
                        break;
                    }
                    --ilo;
                }

                if (ilo == ihi - 1) {
                    --ihi;
                    continue;
                }

                // Wilkinson シフト: H(ihi-1, ihi-1) に近い固有値
                C sigma = H(ihi - 1, ihi - 1);

                // QR ステップ: Givens 回転で H - σI = QR → H ← RQ + σI
                for (int i = ilo; i < ihi - 1; ++i) {
                    C x = H(i, i) - sigma;
                    C y = H(i + 1, i);
                    T r = std::sqrt(calx::normSq(x) + calx::normSq(y));
                    if (r < eps) { ++total_iter; continue; }
                    C c = x / C(r);
                    C s = y / C(r);

                    // 左から回転: H([i, i+1], :) を更新
                    for (int j = ilo; j < static_cast<int>(n); ++j) {
                        C t1 = H(i, j);
                        C t2 = H(i + 1, j);
                        H(i, j)     = calx::conj(c) * t1 + calx::conj(s) * t2;
                        H(i + 1, j) = -s * t1 + c * t2;
                    }
                    // 右から回転: H(:, [i, i+1]) を更新
                    int limit = std::min(i + 3, ihi);
                    for (int j = 0; j < limit; ++j) {
                        C t1 = H(j, i);
                        C t2 = H(j, i + 1);
                        H(j, i)     = c * t1 + s * t2;
                        H(j, i + 1) = -calx::conj(s) * t1 + calx::conj(c) * t2;
                    }
                    // Q 蓄積
                    for (std::size_t j = 0; j < n; ++j) {
                        C t1 = Q(j, i);
                        C t2 = Q(j, i + 1);
                        Q(j, i)     = c * t1 + s * t2;
                        Q(j, i + 1) = -calx::conj(s) * t1 + calx::conj(c) * t2;
                    }
                }
                ++total_iter;
            }

            // Step 3: 固有値を対角から抽出
            Vector<C> eigenvalues(n);
            for (std::size_t i = 0; i < n; ++i)
                eigenvalues[i] = H(i, i);

            // Step 4: 固有ベクトル (後退代入 + Q 変換)
            Matrix<C> eigenvectors(n, n, C(T(0)));
            for (std::size_t k = 0; k < n; ++k) {
                C lambda = eigenvalues[k];
                std::vector<C> v(n, C(T(0)));
                v[k] = C(T(1));

                for (int i = static_cast<int>(k) - 1; i >= 0; --i) {
                    C sum(T(0));
                    for (std::size_t j = i + 1; j <= k; ++j)
                        sum += H(i, j) * v[j];
                    C diag = H(i, i) - lambda;
                    if (calx::abs(diag) > eps * T(100))
                        v[i] = -sum / diag;
                    else
                        v[i] = -sum / C(eps * T(100));
                }

                // Q * v
                for (std::size_t i = 0; i < n; ++i) {
                    C sum(T(0));
                    for (std::size_t j = 0; j < n; ++j)
                        sum += Q(i, j) * v[j];
                    eigenvectors(i, k) = sum;
                }

                // 正規化
                T norm_sq = T(0);
                for (std::size_t i = 0; i < n; ++i)
                    norm_sq += calx::normSq(eigenvectors(i, k));
                T nrm = std::sqrt(norm_sq);
                if (nrm > T(0))
                    for (std::size_t i = 0; i < n; ++i)
                        eigenvectors(i, k) = eigenvectors(i, k) / nrm;
            }

            return {eigenvalues, eigenvectors};
        }

    } // namespace algorithms
} // namespace calx

#endif // CALX_EIGENVALUES_HPP