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

// matrix.hpp
// 行列クラスの定義
//
// このファイルでは、MKL代数ライブラリの行列クラスを定義します。
// 動的サイズ行列や固定サイズ行列の実装を含みます。

#ifndef CALX_MATRIX_HPP
#define CALX_MATRIX_HPP

#include <cassert>
#include <math/core/common.hpp>
#include <math/core/traits.hpp>
#include <math/core/computation_policy.hpp>
#include <math/core/vector.hpp>

#include <algorithm>
#include <array>
#include <bit>
#include <initializer_list>
#include <memory>
#include <numeric>
#include <vector>
#include <iomanip>
#include <sstream>

namespace calx {

    // 前方宣言 (operator^ から参照)
    template<typename T> class Matrix;
    template<typename T>
    [[nodiscard]] Matrix<T> power(const Matrix<T>& base, unsigned int exponent);

    //-----------------------------------------------------------------------------
    // 行列ストレージの定義（detail名前空間内で実装）
    //-----------------------------------------------------------------------------
    namespace detail {

        // 動的サイズの行列ストレージ（std::vector ベース）
        template<typename T>
        class DynamicMatrixStorage {
        public:
            using value_type = T;
            using size_type = std::size_t;
            using reference = T&;
            using const_reference = const T&;
            using pointer = T*;
            using const_pointer = const T*;
            using iterator = typename std::vector<T>::iterator;
            using const_iterator = typename std::vector<T>::const_iterator;

            // デフォルトコンストラクタ
            DynamicMatrixStorage() : rows_(0), cols_(0), data_() {}

            // サイズ指定コンストラクタ
            DynamicMatrixStorage(size_type rows, size_type cols)
                : rows_(rows), cols_(cols), data_(rows* cols) {
            }

            // サイズと初期値指定コンストラクタ
            DynamicMatrixStorage(size_type rows, size_type cols, const T& value)
                : rows_(rows), cols_(cols), data_(rows* cols, value) {
            }

            // コピー/ムーブコンストラクタとアサインメント
            DynamicMatrixStorage(const DynamicMatrixStorage&) = default;
            DynamicMatrixStorage(DynamicMatrixStorage&&) noexcept = default;
            DynamicMatrixStorage& operator=(const DynamicMatrixStorage&) = default;
            DynamicMatrixStorage& operator=(DynamicMatrixStorage&&) noexcept = default;

            // デストラクタ
            ~DynamicMatrixStorage() = default;

            // 要素アクセス (row, col)形式 — Debug 時 assert で境界チェック
            reference operator()(size_type row, size_type col) {
                assert(row < rows_ && "Matrix::operator(): row out of range");
                assert(col < cols_ && "Matrix::operator(): col out of range");
                return data_[row * cols_ + col];
            }

            const_reference operator()(size_type row, size_type col) const {
                assert(row < rows_ && "Matrix::operator(): row out of range");
                assert(col < cols_ && "Matrix::operator(): col out of range");
                return data_[row * cols_ + col];
            }

            // 範囲チェック付き要素アクセス
            reference at(size_type row, size_type col) {
                check_index(row, col);
                return data_[row * cols_ + col];
            }

            const_reference at(size_type row, size_type col) const {
                check_index(row, col);
                return data_[row * cols_ + col];
            }

            // データポインタアクセス
            pointer data() noexcept {
                return data_.data();
            }

            const_pointer data() const noexcept {
                return data_.data();
            }

            // イテレータ
            iterator begin() noexcept {
                return data_.begin();
            }

            const_iterator begin() const noexcept {
                return data_.begin();
            }

            const_iterator cbegin() const noexcept {
                return data_.cbegin();
            }

            iterator end() noexcept {
                return data_.end();
            }

            const_iterator end() const noexcept {
                return data_.end();
            }

            const_iterator cend() const noexcept {
                return data_.cend();
            }

            // 容量
            bool empty() const noexcept {
                return data_.empty();
            }

            size_type size() const noexcept {
                return data_.size();
            }

            size_type rows() const noexcept {
                return rows_;
            }

            size_type cols() const noexcept {
                return cols_;
            }

            size_type capacity() const noexcept {
                return data_.capacity();
            }

            // メモリ確保
            void reserve(size_type new_cap) {
                data_.reserve(new_cap);
            }

            // 縮小
            void shrink_to_fit() {
                data_.shrink_to_fit();
            }

            // サイズ変更
            void resize(size_type rows, size_type cols) {
                rows_ = rows;
                cols_ = cols;
                data_.resize(rows * cols);
            }

            void resize(size_type rows, size_type cols, const value_type& value) {
                rows_ = rows;
                cols_ = cols;
                data_.resize(rows * cols, value);
            }

            // クリア
            void clear() noexcept {
                rows_ = 0;
                cols_ = 0;
                data_.clear();
            }

            // 行数・列数の変更
            void reshape(size_type rows, size_type cols) {
                if (rows * cols != rows_ * cols_) {
                    throw DimensionError("Matrix reshape: total element count must be preserved");
                }
                rows_ = rows;
                cols_ = cols;
            }

            // 正方行列かどうか
            bool is_square() const noexcept {
                return rows_ == cols_;
            }

        private:
            // インデックスチェック
            void check_index(size_type row, size_type col) const {
                if (row >= rows_) {
                    throw IndexError("Matrix row index out of range: " + std::to_string(row) +
                        " >= " + std::to_string(rows_));
                }
                if (col >= cols_) {
                    throw IndexError("Matrix column index out of range: " + std::to_string(col) +
                        " >= " + std::to_string(cols_));
                }
            }

            size_type rows_;
            size_type cols_;
            std::vector<T> data_;
        };

        // 固定サイズの行列ストレージ（std::array ベース）
        template<typename T, std::size_t Rows, std::size_t Cols>
        class StaticMatrixStorage {
        public:
            using value_type = T;
            using size_type = std::size_t;
            using reference = T&;
            using const_reference = const T&;
            using pointer = T*;
            using const_pointer = const T*;
            using iterator = pointer;
            using const_iterator = const_pointer;

            // デフォルトコンストラクタ（ゼロ初期化）
            StaticMatrixStorage() : data_{} {}

            // 値で初期化するコンストラクタ
            explicit StaticMatrixStorage(const T& value) {
                std::fill(begin(), end(), value);
            }

            // コピー/ムーブコンストラクタとアサインメント
            StaticMatrixStorage(const StaticMatrixStorage&) = default;
            StaticMatrixStorage(StaticMatrixStorage&&) noexcept = default;
            StaticMatrixStorage& operator=(const StaticMatrixStorage&) = default;
            StaticMatrixStorage& operator=(StaticMatrixStorage&&) noexcept = default;

            // デストラクタ
            ~StaticMatrixStorage() = default;

            // 要素アクセス (row, col)形式 — Debug 時 assert で境界チェック
            reference operator()(size_type row, size_type col) {
                assert(row < Rows && "StaticMatrix::operator(): row out of range");
                assert(col < Cols && "StaticMatrix::operator(): col out of range");
                return data_[row * Cols + col];
            }

            const_reference operator()(size_type row, size_type col) const {
                assert(row < Rows && "StaticMatrix::operator(): row out of range");
                assert(col < Cols && "StaticMatrix::operator(): col out of range");
                return data_[row * Cols + col];
            }

            // 範囲チェック付き要素アクセス
            reference at(size_type row, size_type col) {
                check_index(row, col);
                return data_[row * Cols + col];
            }

            const_reference at(size_type row, size_type col) const {
                check_index(row, col);
                return data_[row * Cols + col];
            }

            // データポインタアクセス
            pointer data() noexcept {
                return data_.data();
            }

            const_pointer data() const noexcept {
                return data_.data();
            }

            // イテレータ
            iterator begin() noexcept {
                return data_.data();
            }

            const_iterator begin() const noexcept {
                return data_.data();
            }

            const_iterator cbegin() const noexcept {
                return data_.data();
            }

            iterator end() noexcept {
                return data_.data() + Rows * Cols;
            }

            const_iterator end() const noexcept {
                return data_.data() + Rows * Cols;
            }

            const_iterator cend() const noexcept {
                return data_.data() + Rows * Cols;
            }

            // 容量
            constexpr bool empty() const noexcept {
                return size() == 0;
            }

            constexpr size_type size() const noexcept {
                return Rows * Cols;
            }

            constexpr size_type rows() const noexcept {
                return Rows;
            }

            constexpr size_type cols() const noexcept {
                return Cols;
            }

            // 正方行列かどうか
            constexpr bool is_square() const noexcept {
                return Rows == Cols;
            }

        private:
            // インデックスチェック
            void check_index(size_type row, size_type col) const {
                if (row >= Rows) {
                    throw IndexError("StaticMatrix row index out of range: " + std::to_string(row) +
                        " >= " + std::to_string(Rows));
                }
                if (col >= Cols) {
                    throw IndexError("StaticMatrix column index out of range: " + std::to_string(col) +
                        " >= " + std::to_string(Cols));
                }
            }

            std::array<T, Rows* Cols> data_;
        };

        // アライメント済み行列ストレージ（MKL最適化）
        template<typename T>
        class AlignedMatrixStorage {
        public:
            using value_type = T;
            using size_type = std::size_t;
            using reference = T&;
            using const_reference = const T&;
            using pointer = T*;
            using const_pointer = const T*;
            using iterator = pointer;
            using const_iterator = const_pointer;

            // アライメント値
            static constexpr size_type alignment = memory::simd_alignment();

            // デフォルトコンストラクタ
            AlignedMatrixStorage() : data_(nullptr), rows_(0), cols_(0), capacity_(0) {}

            // サイズ指定コンストラクタ
            AlignedMatrixStorage(size_type rows, size_type cols)
                : data_(nullptr), rows_(rows), cols_(cols), capacity_(rows* cols) {
                if (capacity_ > 0) {
                    data_ = memory::aligned_alloc<T>(capacity_, alignment);
                    std::uninitialized_default_construct_n(data_, capacity_);
                }
            }

            // サイズと初期値指定コンストラクタ
            AlignedMatrixStorage(size_type rows, size_type cols, const T& value)
                : data_(nullptr), rows_(rows), cols_(cols), capacity_(rows* cols) {
                if (capacity_ > 0) {
                    data_ = memory::aligned_alloc<T>(capacity_, alignment);
                    std::uninitialized_fill_n(data_, capacity_, value);
                }
            }

            // コピーコンストラクタ
            AlignedMatrixStorage(const AlignedMatrixStorage& other)
                : data_(nullptr), rows_(other.rows_), cols_(other.cols_), capacity_(other.capacity_) {
                if (capacity_ > 0) {
                    data_ = memory::aligned_alloc<T>(capacity_, alignment);
                    std::uninitialized_copy_n(other.data_, capacity_, data_);
                }
            }

            // ムーブコンストラクタ
            AlignedMatrixStorage(AlignedMatrixStorage&& other) noexcept
                : data_(other.data_), rows_(other.rows_), cols_(other.cols_), capacity_(other.capacity_) {
                other.data_ = nullptr;
                other.rows_ = 0;
                other.cols_ = 0;
                other.capacity_ = 0;
            }

            // コピー代入演算子
            AlignedMatrixStorage& operator=(const AlignedMatrixStorage& other) {
                if (this == &other) return *this;

                // 古いデータの破棄
                clear();
                if (data_) {
                    memory::aligned_free(data_);
                }

                // 新しいデータの確保
                rows_ = other.rows_;
                cols_ = other.cols_;
                capacity_ = other.capacity_;

                if (capacity_ > 0) {
                    data_ = memory::aligned_alloc<T>(capacity_, alignment);
                    std::uninitialized_copy_n(other.data_, capacity_, data_);
                }
                else {
                    data_ = nullptr;
                }

                return *this;
            }

            // ムーブ代入演算子
            AlignedMatrixStorage& operator=(AlignedMatrixStorage&& other) noexcept {
                if (this == &other) return *this;

                // 古いデータの破棄
                clear();
                if (data_) {
                    memory::aligned_free(data_);
                }

                // 他方のデータを移動
                data_ = other.data_;
                rows_ = other.rows_;
                cols_ = other.cols_;
                capacity_ = other.capacity_;

                other.data_ = nullptr;
                other.rows_ = 0;
                other.cols_ = 0;
                other.capacity_ = 0;

                return *this;
            }

            // デストラクタ
            ~AlignedMatrixStorage() {
                clear();
                if (data_) {
                    memory::aligned_free(data_);
                }
            }

            // 要素アクセス (row, col)形式
            reference operator()(size_type row, size_type col) {
                return data_[row * cols_ + col];
            }

            const_reference operator()(size_type row, size_type col) const {
                return data_[row * cols_ + col];
            }

            // 範囲チェック付き要素アクセス
            reference at(size_type row, size_type col) {
                check_index(row, col);
                return data_[row * cols_ + col];
            }

            const_reference at(size_type row, size_type col) const {
                check_index(row, col);
                return data_[row * cols_ + col];
            }

            // データポインタアクセス
            pointer data() noexcept {
                return data_;
            }

            const_pointer data() const noexcept {
                return data_;
            }

            // イテレータ
            iterator begin() noexcept {
                return data_;
            }

            const_iterator begin() const noexcept {
                return data_;
            }

            const_iterator cbegin() const noexcept {
                return data_;
            }

            iterator end() noexcept {
                return data_ + rows_ * cols_;
            }

            const_iterator end() const noexcept {
                return data_ + rows_ * cols_;
            }

            const_iterator cend() const noexcept {
                return data_ + rows_ * cols_;
            }

            // 容量
            bool empty() const noexcept {
                return rows_ == 0 || cols_ == 0;
            }

            size_type size() const noexcept {
                return rows_ * cols_;
            }

            size_type rows() const noexcept {
                return rows_;
            }

            size_type cols() const noexcept {
                return cols_;
            }

            size_type capacity() const noexcept {
                return capacity_;
            }

            // メモリ確保
            void reserve(size_type new_cap) {
                if (new_cap <= capacity_) return;

                // 新しいメモリの確保
                pointer new_data = memory::aligned_alloc<T>(new_cap, alignment);

                // 古いデータの移動
                if (size() > 0) {
                    std::uninitialized_move_n(data_, size(), new_data);
                    // 古いオブジェクトの破棄
                    std::destroy_n(data_, size());
                }

                // 古いメモリの解放
                if (data_) {
                    memory::aligned_free(data_);
                }

                // 新しいメモリの設定
                data_ = new_data;
                capacity_ = new_cap;
            }

            // クリア
            void clear() noexcept {
                if (data_) {
                    std::destroy_n(data_, size());
                }
                rows_ = 0;
                cols_ = 0;
            }

            // サイズ変更
            void resize(size_type rows, size_type cols) {
                size_type new_size = rows * cols;

                if (new_size > capacity_) {
                    // 拡大
                    reserve(new_size);
                }

                if (new_size > size()) {
                    // 要素の追加
                    std::uninitialized_default_construct_n(data_ + size(), new_size - size());
                }
                else if (new_size < size()) {
                    // 要素の削除
                    std::destroy_n(data_ + new_size, size() - new_size);
                }

                rows_ = rows;
                cols_ = cols;
            }

            void resize(size_type rows, size_type cols, const value_type& value) {
                size_type new_size = rows * cols;

                if (new_size > capacity_) {
                    // 拡大
                    reserve(new_size);
                }

                if (new_size > size()) {
                    // 要素の追加
                    std::uninitialized_fill_n(data_ + size(), new_size - size(), value);
                }
                else if (new_size < size()) {
                    // 要素の削除
                    std::destroy_n(data_ + new_size, size() - new_size);
                }

                rows_ = rows;
                cols_ = cols;
            }

            // 行数・列数の変更
            void reshape(size_type rows, size_type cols) {
                if (rows * cols != rows_ * cols_) {
                    throw DimensionError("Matrix reshape: total element count must be preserved");
                }
                rows_ = rows;
                cols_ = cols;
            }

            // 正方行列かどうか
            bool is_square() const noexcept {
                return rows_ == cols_;
            }

        private:
            // インデックスチェック
            void check_index(size_type row, size_type col) const {
                if (row >= rows_) {
                    throw IndexError("AlignedMatrix row index out of range: " + std::to_string(row) +
                        " >= " + std::to_string(rows_));
                }
                if (col >= cols_) {
                    throw IndexError("AlignedMatrix column index out of range: " + std::to_string(col) +
                        " >= " + std::to_string(cols_));
                }
            }

            pointer data_;
            size_type rows_;
            size_type cols_;
            size_type capacity_;
        };

        // 行列ビューの実装
        template<typename MatrixType>
        class MatrixView {
        public:
            // ホスト型から型定義を取得
            using value_type = typename std::remove_reference_t<MatrixType>::value_type;
            using size_type = std::size_t;
            using reference = value_type&;
            using const_reference = const value_type&;

            // 元の行列への参照型
            using matrix_reference = MatrixType&;

            // 完全な行列のビューを作成
            explicit MatrixView(matrix_reference matrix)
                : matrix_(matrix), row_offset_(0), col_offset_(0),
                num_rows_(matrix.rows()), num_cols_(matrix.cols()) {
            }

            // 行列の部分ビューを作成
            MatrixView(matrix_reference matrix,
                size_type row_offset, size_type col_offset,
                size_type num_rows, size_type num_cols)
                : matrix_(matrix), row_offset_(row_offset), col_offset_(col_offset),
                num_rows_(num_rows), num_cols_(num_cols) {

                // 範囲チェック
                if (row_offset + num_rows > matrix.rows() || col_offset + num_cols > matrix.cols()) {
                    throw IndexError("MatrixView: view range exceeds matrix bounds");
                }
            }

            // コピーコンストラクタ
            MatrixView(const MatrixView&) = default;

            // コピー代入演算子は禁止（ビュー自体のコピーは禁止）
            MatrixView& operator=(const MatrixView&) = delete;

            // 行列の内容をコピーする代入
            template<typename OtherMatrixType>
            MatrixView& operator=(const OtherMatrixType& other) {
                assign(other);
                return *this;
            }

            // 要素アクセス
            reference operator()(size_type row, size_type col) {
                return matrix_(row_offset_ + row, col_offset_ + col);
            }

            const_reference operator()(size_type row, size_type col) const {
                return matrix_(row_offset_ + row, col_offset_ + col);
            }

            // 範囲チェック付き要素アクセス
            reference at(size_type row, size_type col) {
                if (row >= num_rows_ || col >= num_cols_) {
                    throw IndexError("MatrixView::at: index out of range");
                }
                return matrix_(row_offset_ + row, col_offset_ + col);
            }

            const_reference at(size_type row, size_type col) const {
                if (row >= num_rows_ || col >= num_cols_) {
                    throw IndexError("MatrixView::at: index out of range");
                }
                return matrix_(row_offset_ + row, col_offset_ + col);
            }

            // 容量関連
            bool empty() const noexcept {
                return num_rows_ == 0 || num_cols_ == 0;
            }

            size_type rows() const noexcept {
                return num_rows_;
            }

            size_type cols() const noexcept {
                return num_cols_;
            }

            size_type size() const noexcept {
                return num_rows_ * num_cols_;
            }

            // オフセットの取得
            size_type row_offset() const noexcept {
                return row_offset_;
            }

            size_type col_offset() const noexcept {
                return col_offset_;
            }

            // 元の行列の取得
            matrix_reference parent() noexcept {
                return matrix_;
            }

            const matrix_reference parent() const noexcept {
                return matrix_;
            }

            // サブビューの作成
            MatrixView subview(size_type row_offset, size_type col_offset,
                size_type num_rows, size_type num_cols) const {
                if (row_offset + num_rows > num_rows_ || col_offset + num_cols > num_cols_) {
                    throw IndexError("MatrixView::subview: subview range exceeds view bounds");
                }

                return MatrixView(matrix_,
                    row_offset_ + row_offset,
                    col_offset_ + col_offset,
                    num_rows, num_cols);
            }

            // 行の取得
            template<typename VectorType>
            void get_row(size_type row_idx, VectorType& dest) const {
                if (row_idx >= num_rows_) {
                    throw IndexError("MatrixView::get_row: row index out of range");
                }

                if (dest.size() != num_cols_) {
                    dest.resize(num_cols_);
                }

                for (size_type j = 0; j < num_cols_; ++j) {
                    dest[j] = matrix_(row_offset_ + row_idx, col_offset_ + j);
                }
            }

            // 列の取得
            template<typename VectorType>
            void get_column(size_type col_idx, VectorType& dest) const {
                if (col_idx >= num_cols_) {
                    throw IndexError("MatrixView::get_column: column index out of range");
                }

                if (dest.size() != num_rows_) {
                    dest.resize(num_rows_);
                }

                for (size_type i = 0; i < num_rows_; ++i) {
                    dest[i] = matrix_(row_offset_ + i, col_offset_ + col_idx);
                }
            }

            // 行の設定
            template<typename VectorType>
            void set_row(size_type row_idx, const VectorType& src) {
                if (row_idx >= num_rows_) {
                    throw IndexError("MatrixView::set_row: row index out of range");
                }

                if (src.size() != num_cols_) {
                    throw DimensionError("MatrixView::set_row: vector size mismatch");
                }

                for (size_type j = 0; j < num_cols_; ++j) {
                    matrix_(row_offset_ + row_idx, col_offset_ + j) = src[j];
                }
            }

            // 列の設定
            template<typename VectorType>
            void set_column(size_type col_idx, const VectorType& src) {
                if (col_idx >= num_cols_) {
                    throw IndexError("MatrixView::set_column: column index out of range");
                }

                if (src.size() != num_rows_) {
                    throw DimensionError("MatrixView::set_column: vector size mismatch");
                }

                for (size_type i = 0; i < num_rows_; ++i) {
                    matrix_(row_offset_ + i, col_offset_ + col_idx) = src[i];
                }
            }

            // 値のコピー
            void fill(const value_type& value) {
                for (size_type i = 0; i < num_rows_; ++i) {
                    for (size_type j = 0; j < num_cols_; ++j) {
                        matrix_(row_offset_ + i, col_offset_ + j) = value;
                    }
                }
            }

            // ゼロクリア
            void zero() {
                fill(numeric_traits<value_type>::zero());
            }

            // 他の行列からの代入
            template<typename OtherMatrixType>
            void assign(const OtherMatrixType& src) {
                if (src.rows() != num_rows_ || src.cols() != num_cols_) {
                    throw DimensionError("MatrixView::assign: matrix size mismatch");
                }

                for (size_type i = 0; i < num_rows_; ++i) {
                    for (size_type j = 0; j < num_cols_; ++j) {
                        matrix_(row_offset_ + i, col_offset_ + j) = src(i, j);
                    }
                }
            }

            // 正方行列かどうか
            bool is_square() const noexcept {
                return num_rows_ == num_cols_;
            }

        private:
            matrix_reference matrix_;  // 元の行列への参照
            size_type row_offset_;     // 行オフセット
            size_type col_offset_;     // 列オフセット
            size_type num_rows_;       // ビューの行数
            size_type num_cols_;       // ビューの列数
        };

        // 読み取り専用行列ビュークラス
        template<typename MatrixType>
        class ConstMatrixView {
        public:
            // ホスト型から型定義を取得
            using value_type = typename std::remove_reference_t<MatrixType>::value_type;
            using size_type = std::size_t;
            using reference = const value_type&;  // constリファレンス
            using const_reference = const value_type&;

            // 元の行列への参照型（常にconst）
            using matrix_reference = const MatrixType&;

            // 完全な行列のビューを作成
            explicit ConstMatrixView(matrix_reference matrix)
                : matrix_(matrix), row_offset_(0), col_offset_(0),
                num_rows_(matrix.rows()), num_cols_(matrix.cols()) {
            }

            // 行列の部分ビューを作成
            ConstMatrixView(matrix_reference matrix,
                size_type row_offset, size_type col_offset,
                size_type num_rows, size_type num_cols)
                : matrix_(matrix), row_offset_(row_offset), col_offset_(col_offset),
                num_rows_(num_rows), num_cols_(num_cols) {

                // 範囲チェック
                if (row_offset + num_rows > matrix.rows() || col_offset + num_cols > matrix.cols()) {
                    throw IndexError("ConstMatrixView: view range exceeds matrix bounds");
                }
            }

            // MatrixViewからの変換コンストラクタ
            explicit ConstMatrixView(const MatrixView<MatrixType>& view)
                : matrix_(view.parent()),
                row_offset_(view.row_offset()),
                col_offset_(view.col_offset()),
                num_rows_(view.rows()),
                num_cols_(view.cols()) {
            }

            // コピーコンストラクタ
            ConstMatrixView(const ConstMatrixView&) = default;

            // コピー代入演算子は禁止（ビュー自体のコピーは禁止）
            ConstMatrixView& operator=(const ConstMatrixView&) = delete;

            // 要素アクセス（常にconst）
            const_reference operator()(size_type row, size_type col) const {
                return matrix_(row_offset_ + row, col_offset_ + col);
            }

            // 範囲チェック付き要素アクセス（常にconst）
            const_reference at(size_type row, size_type col) const {
                if (row >= num_rows_ || col >= num_cols_) {
                    throw IndexError("ConstMatrixView::at: index out of range");
                }
                return matrix_(row_offset_ + row, col_offset_ + col);
            }

            // 容量関連
            bool empty() const noexcept {
                return num_rows_ == 0 || num_cols_ == 0;
            }

            size_type rows() const noexcept {
                return num_rows_;
            }

            size_type cols() const noexcept {
                return num_cols_;
            }

            size_type size() const noexcept {
                return num_rows_ * num_cols_;
            }

            // オフセットの取得
            size_type row_offset() const noexcept {
                return row_offset_;
            }

            size_type col_offset() const noexcept {
                return col_offset_;
            }

            // 元の行列の取得（常にconst）
            matrix_reference parent() const noexcept {
                return matrix_;
            }

            // サブビューの作成
            ConstMatrixView subview(size_type row_offset, size_type col_offset,
                size_type num_rows, size_type num_cols) const {
                if (row_offset + num_rows > num_rows_ || col_offset + num_cols > num_cols_) {
                    throw IndexError("ConstMatrixView::subview: subview range exceeds view bounds");
                }

                return ConstMatrixView(matrix_,
                    row_offset_ + row_offset,
                    col_offset_ + col_offset,
                    num_rows, num_cols);
            }

            // 行の取得
            template<typename VectorType>
            void get_row(size_type row_idx, VectorType& dest) const {
                if (row_idx >= num_rows_) {
                    throw IndexError("ConstMatrixView::get_row: row index out of range");
                }

                if (dest.size() != num_cols_) {
                    dest.resize(num_cols_);
                }

                for (size_type j = 0; j < num_cols_; ++j) {
                    dest[j] = matrix_(row_offset_ + row_idx, col_offset_ + j);
                }
            }

            // 列の取得
            template<typename VectorType>
            void get_column(size_type col_idx, VectorType& dest) const {
                if (col_idx >= num_cols_) {
                    throw IndexError("ConstMatrixView::get_column: column index out of range");
                }

                if (dest.size() != num_rows_) {
                    dest.resize(num_rows_);
                }

                for (size_type i = 0; i < num_rows_; ++i) {
                    dest[i] = matrix_(row_offset_ + i, col_offset_ + col_idx);
                }
            }

            // 正方行列かどうか
            bool is_square() const noexcept {
                return num_rows_ == num_cols_;
            }

        private:
            matrix_reference matrix_;  // 元の行列への参照（const）
            size_type row_offset_;     // 行オフセット
            size_type col_offset_;     // 列オフセット
            size_type num_rows_;       // ビューの行数
            size_type num_cols_;       // ビューの列数
        };

        // ビュー作成のためのヘルパー関数
        template<typename MatrixType>
        MatrixView<MatrixType> make_matrix_view(MatrixType& matrix) {
            return MatrixView<MatrixType>(matrix);
        }

        template<typename MatrixType>
        MatrixView<MatrixType> make_matrix_view(MatrixType& matrix,
            std::size_t row_offset, std::size_t col_offset,
            std::size_t num_rows, std::size_t num_cols) {
            return MatrixView<MatrixType>(matrix, row_offset, col_offset, num_rows, num_cols);
        }

        template<typename MatrixType>
        ConstMatrixView<MatrixType> make_const_matrix_view(const MatrixType& matrix) {
            return ConstMatrixView<MatrixType>(matrix);
        }

        template<typename MatrixType>
        ConstMatrixView<MatrixType> make_const_matrix_view(const MatrixType& matrix,
            std::size_t row_offset, std::size_t col_offset,
            std::size_t num_rows, std::size_t num_cols) {
            return ConstMatrixView<MatrixType>(matrix, row_offset, col_offset, num_rows, num_cols);
        }

        // 行列基本クラス
        template<typename T, typename StoragePolicy>
        class MatrixBase {
        public:
            using value_type = T;
            using size_type = std::size_t;
            using reference = T&;
            using const_reference = const T&;
            using iterator = typename StoragePolicy::iterator;
            using const_iterator = typename StoragePolicy::const_iterator;

            // デフォルトコンストラクタ
            MatrixBase() = default;

            // サイズ指定コンストラクタ
            explicit MatrixBase(size_type rows, size_type cols) : storage_(rows, cols) {}

            // サイズと初期値指定コンストラクタ
            MatrixBase(size_type rows, size_type cols, const T& value) : storage_(rows, cols, value) {}

            // コピー/ムーブコンストラクタとアサインメント
            MatrixBase(const MatrixBase&) = default;
            MatrixBase(MatrixBase&&) noexcept = default;
            MatrixBase& operator=(const MatrixBase&) = default;
            MatrixBase& operator=(MatrixBase&&) noexcept = default;

            // デストラクタ
            ~MatrixBase() = default;

            // 要素アクセス
            reference operator()(size_type row, size_type col) {
                return storage_(row, col);
            }

            const_reference operator()(size_type row, size_type col) const {
                return storage_(row, col);
            }

            // 範囲チェック付き要素アクセス
            reference at(size_type row, size_type col) {
                return storage_.at(row, col);
            }

            const_reference at(size_type row, size_type col) const {
                return storage_.at(row, col);
            }

            // イテレータ
            iterator begin() noexcept {
                return storage_.begin();
            }

            const_iterator begin() const noexcept {
                return storage_.begin();
            }

            const_iterator cbegin() const noexcept {
                return storage_.cbegin();
            }

            iterator end() noexcept {
                return storage_.end();
            }

            const_iterator end() const noexcept {
                return storage_.end();
            }

            const_iterator cend() const noexcept {
                return storage_.cend();
            }

            // 容量
            bool empty() const noexcept {
                return storage_.empty();
            }

            size_type rows() const noexcept {
                return storage_.rows();
            }

            size_type cols() const noexcept {
                return storage_.cols();
            }

            size_type size() const noexcept {
                return storage_.size();
            }

            // サイズ変更
            void resize(size_type rows, size_type cols) {
                storage_.resize(rows, cols);
            }

            void resize(size_type rows, size_type cols, const value_type& value) {
                storage_.resize(rows, cols, value);
            }

            // クリア
            void clear() noexcept {
                storage_.clear();
            }

            // 行数・列数の変更
            void reshape(size_type rows, size_type cols) {
                storage_.reshape(rows, cols);
            }

            // 全要素の値設定
            void fill(const T& value) {
                for (size_type i = 0; i < rows(); ++i) {
                    for (size_type j = 0; j < cols(); ++j) {
                        storage_(i, j) = value;
                    }
                }
            }

            // ゼロクリア
            void zero() {
                fill(numeric_traits<T>::zero());
            }

            // 単位行列化
            void identity() {
                if (!is_square()) {
                    throw DimensionError("identity(): matrix must be square");
                }

                zero();

                const size_type n = rows();
                for (size_type i = 0; i < n; ++i) {
                    storage_(i, i) = numeric_traits<T>::one();
                }
            }

            // 行の取得
            template<typename VectorType>
            void get_row(size_type row_idx, VectorType& dest) const {
                if (row_idx >= rows()) {
                    throw IndexError("get_row: row index out of range");
                }

                if (dest.size() != cols()) {
                    dest.resize(cols());
                }

                for (size_type j = 0; j < cols(); ++j) {
                    dest[j] = storage_(row_idx, j);
                }
            }

            // 列の取得
            template<typename VectorType>
            void get_column(size_type col_idx, VectorType& dest) const {
                if (col_idx >= cols()) {
                    throw IndexError("get_column: column index out of range");
                }

                if (dest.size() != rows()) {
                    dest.resize(rows());
                }

                for (size_type i = 0; i < rows(); ++i) {
                    dest[i] = storage_(i, col_idx);
                }
            }

            // 行の設定
            template<typename VectorType>
            void set_row(size_type row_idx, const VectorType& src) {
                if (row_idx >= rows()) {
                    throw IndexError("set_row: row index out of range");
                }

                if (src.size() != cols()) {
                    throw DimensionError("set_row: vector size mismatch");
                }

                for (size_type j = 0; j < cols(); ++j) {
                    storage_(row_idx, j) = src[j];
                }
            }

            // 列の設定
            template<typename VectorType>
            void set_column(size_type col_idx, const VectorType& src) {
                if (col_idx >= cols()) {
                    throw IndexError("set_column: column index out of range");
                }

                if (src.size() != rows()) {
                    throw DimensionError("set_column: vector size mismatch");
                }

                for (size_type i = 0; i < rows(); ++i) {
                    storage_(i, col_idx) = src[i];
                }
            }

            // ストレージの取得（内部実装用）
            StoragePolicy& storage() noexcept {
                return storage_;
            }

            const StoragePolicy& storage() const noexcept {
                return storage_;
            }

            // 内部配列への raw ポインタ (row-major 連続メモリ)
            T* data() noexcept {
                return storage_.data();
            }

            const T* data() const noexcept {
                return storage_.data();
            }

            // 正方行列かどうか
            bool is_square() const noexcept {
                return storage_.is_square();
            }

            // 部分行列ビューの取得
            MatrixView<MatrixBase> submatrix(size_type row_offset, size_type col_offset,
                size_type num_rows, size_type num_cols) {
                return MatrixView<MatrixBase>(*this, row_offset, col_offset, num_rows, num_cols);
            }

            ConstMatrixView<MatrixBase> submatrix(size_type row_offset, size_type col_offset,
                size_type num_rows, size_type num_cols) const {
                return ConstMatrixView<MatrixBase>(*this, row_offset, col_offset, num_rows, num_cols);
            }

            // 行ビューの取得
            MatrixView<MatrixBase> row(size_type row_idx) {
                if (row_idx >= rows()) {
                    throw IndexError("row: row index out of range");
                }
                return MatrixView<MatrixBase>(*this, row_idx, 0, 1, cols());
            }

            ConstMatrixView<MatrixBase> row(size_type row_idx) const {
                if (row_idx >= rows()) {
                    throw IndexError("row: row index out of range");
                }
                return ConstMatrixView<MatrixBase>(*this, row_idx, 0, 1, cols());
            }

            // 列ビューの取得
            MatrixView<MatrixBase> column(size_type col_idx) {
                if (col_idx >= cols()) {
                    throw IndexError("column: column index out of range");
                }
                return MatrixView<MatrixBase>(*this, 0, col_idx, rows(), 1);
            }

            ConstMatrixView<MatrixBase> column(size_type col_idx) const {
                if (col_idx >= cols()) {
                    throw IndexError("column: column index out of range");
                }
                return ConstMatrixView<MatrixBase>(*this, 0, col_idx, rows(), 1);
            }

            // ブロック操作 API (Eigen スタイル)

            /// block(i, j, p, q): submatrix の別名
            MatrixView<MatrixBase> block(size_type i, size_type j, size_type p, size_type q) {
                return submatrix(i, j, p, q);
            }
            ConstMatrixView<MatrixBase> block(size_type i, size_type j, size_type p, size_type q) const {
                return submatrix(i, j, p, q);
            }

            /// 左上コーナー (0,0) から p 行 q 列
            MatrixView<MatrixBase> topLeftCorner(size_type p, size_type q) {
                return submatrix(0, 0, p, q);
            }
            ConstMatrixView<MatrixBase> topLeftCorner(size_type p, size_type q) const {
                return submatrix(0, 0, p, q);
            }

            /// 右上コーナー
            MatrixView<MatrixBase> topRightCorner(size_type p, size_type q) {
                return submatrix(0, cols() - q, p, q);
            }
            ConstMatrixView<MatrixBase> topRightCorner(size_type p, size_type q) const {
                return submatrix(0, cols() - q, p, q);
            }

            /// 左下コーナー
            MatrixView<MatrixBase> bottomLeftCorner(size_type p, size_type q) {
                return submatrix(rows() - p, 0, p, q);
            }
            ConstMatrixView<MatrixBase> bottomLeftCorner(size_type p, size_type q) const {
                return submatrix(rows() - p, 0, p, q);
            }

            /// 右下コーナー
            MatrixView<MatrixBase> bottomRightCorner(size_type p, size_type q) {
                return submatrix(rows() - p, cols() - q, p, q);
            }
            ConstMatrixView<MatrixBase> bottomRightCorner(size_type p, size_type q) const {
                return submatrix(rows() - p, cols() - q, p, q);
            }

            /// 上から n 行
            MatrixView<MatrixBase> topRows(size_type n) {
                return submatrix(0, 0, n, cols());
            }
            ConstMatrixView<MatrixBase> topRows(size_type n) const {
                return submatrix(0, 0, n, cols());
            }

            /// 下から n 行
            MatrixView<MatrixBase> bottomRows(size_type n) {
                return submatrix(rows() - n, 0, n, cols());
            }
            ConstMatrixView<MatrixBase> bottomRows(size_type n) const {
                return submatrix(rows() - n, 0, n, cols());
            }

            /// 左から n 列
            MatrixView<MatrixBase> leftCols(size_type n) {
                return submatrix(0, 0, rows(), n);
            }
            ConstMatrixView<MatrixBase> leftCols(size_type n) const {
                return submatrix(0, 0, rows(), n);
            }

            /// 右から n 列
            MatrixView<MatrixBase> rightCols(size_type n) {
                return submatrix(0, cols() - n, rows(), n);
            }
            ConstMatrixView<MatrixBase> rightCols(size_type n) const {
                return submatrix(0, cols() - n, rows(), n);
            }

        protected:
            StoragePolicy storage_;
        };

    } // namespace detail

    //-----------------------------------------------------------------------------
    // 三角・対称行列ビュー
    //-----------------------------------------------------------------------------

    /// 三角行列の種別
    enum class TriangularMode { Upper, Lower, UnitUpper, UnitLower, StrictlyUpper, StrictlyLower };

    /**
     * @brief TriangularView — 三角行列ビュー (Eigen 互換)
     *
     * 行列の上三角/下三角をコピーなしで参照し、三角系の演算を行う。
     * - solve(b): 三角連立方程式 Ax=b を前方/後退代入で解く
     * - operator*(v): 三角行列ベクトル積 (ゼロ要素をスキップ)
     * - toDense(): 三角部分のみをコピーした密行列を生成
     */
    template<typename MatrixType, TriangularMode Mode>
    class TriangularView {
    public:
        using T = typename std::remove_cvref_t<MatrixType>::value_type;
        using size_type = std::size_t;

        static constexpr bool is_upper = (Mode == TriangularMode::Upper ||
                                          Mode == TriangularMode::UnitUpper ||
                                          Mode == TriangularMode::StrictlyUpper);
        static constexpr bool is_unit  = (Mode == TriangularMode::UnitUpper ||
                                          Mode == TriangularMode::UnitLower);
        static constexpr bool is_strict = (Mode == TriangularMode::StrictlyUpper ||
                                           Mode == TriangularMode::StrictlyLower);

        explicit TriangularView(MatrixType& mat) : mat_(mat) {
            if (mat_.rows() != mat_.cols())
                throw DimensionError("TriangularView: matrix must be square");
        }

        size_type rows() const { return mat_.rows(); }
        size_type cols() const { return mat_.cols(); }

        /// 要素アクセス (三角領域外はゼロ, 単位対角は 1)
        T operator()(size_type i, size_type j) const {
            if constexpr (is_upper) {
                if (i > j) return T(0);
                if (is_strict && i == j) return T(0);
                if (is_unit && i == j) return T(1);
            } else {
                if (i < j) return T(0);
                if (is_strict && i == j) return T(0);
                if (is_unit && i == j) return T(1);
            }
            return mat_(i, j);
        }

        /// 三角連立方程式を解く: this * x = b → x
        /// AVX2 FMA による内積 + 行方向後退/前方代入
        Vector<T> solve(const Vector<T>& b) const {
            const auto n = rows();
            if (b.size() != n)
                throw DimensionError("TriangularView::solve: dimension mismatch");
            Vector<T> x = b;
            T* __restrict xp = x.data();
            const T* __restrict A = mat_.data();
            using PT = PacketTraits<T>;
            constexpr size_type W = PT::size;  // AVX2 double: 4

            if constexpr (is_upper) {
                for (size_type ii = 0; ii < n; ++ii) {
                    size_type i = n - 1 - ii;
                    const T* __restrict row = A + i * n;
                    size_type j = i + 1;
                    size_type len = n - j;

                    if constexpr (W > 1) {
                        auto acc0 = PT::set1(T(0));
                        auto acc1 = PT::set1(T(0));
                        size_type end2w = j + (len & ~(2 * W - 1));
                        for (; j < end2w; j += 2 * W) {
                            acc0 = PT::fmadd(PT::load(row + j),     PT::load(xp + j),     acc0);
                            acc1 = PT::fmadd(PT::load(row + j + W), PT::load(xp + j + W), acc1);
                        }
                        size_type endw = i + 1 + (len & ~(W - 1));
                        for (; j < endw; j += W)
                            acc0 = PT::fmadd(PT::load(row + j), PT::load(xp + j), acc0);
                        T sum = PT::reduce_add(PT::add(acc0, acc1));
                        for (; j < n; ++j) sum += row[j] * xp[j];
                        xp[i] -= sum;
                    } else {
                        T sum = T(0);
                        for (; j < n; ++j) sum += row[j] * xp[j];
                        xp[i] -= sum;
                    }
                    if constexpr (!is_unit && !is_strict)
                        xp[i] /= row[i];
                }
            } else {
                for (size_type i = 0; i < n; ++i) {
                    const T* __restrict row = A + i * n;
                    size_type j = 0;

                    if constexpr (W > 1) {
                        auto acc0 = PT::set1(T(0));
                        auto acc1 = PT::set1(T(0));
                        size_type end2w = (i / (2 * W)) * (2 * W);
                        for (; j < end2w; j += 2 * W) {
                            acc0 = PT::fmadd(PT::load(row + j),     PT::load(xp + j),     acc0);
                            acc1 = PT::fmadd(PT::load(row + j + W), PT::load(xp + j + W), acc1);
                        }
                        size_type endw = (i / W) * W;
                        for (; j < endw; j += W)
                            acc0 = PT::fmadd(PT::load(row + j), PT::load(xp + j), acc0);
                        T sum = PT::reduce_add(PT::add(acc0, acc1));
                        for (; j < i; ++j) sum += row[j] * xp[j];
                        xp[i] -= sum;
                    } else {
                        T sum = T(0);
                        for (; j < i; ++j) sum += row[j] * xp[j];
                        xp[i] -= sum;
                    }
                    if constexpr (!is_unit && !is_strict)
                        xp[i] /= row[i];
                }
            }
            return x;
        }

        /// 行列ベクトル積: y = this * x (三角部分のみ)
        Vector<T> operator*(const Vector<T>& x) const {
            const auto n = rows();
            if (x.size() != n)
                throw DimensionError("TriangularView::operator*: dimension mismatch");
            Vector<T> y(n, T(0));
            for (size_type i = 0; i < n; ++i) {
                size_type jstart = is_upper ? i : 0;
                size_type jend   = is_upper ? n : i + 1;
                if constexpr (is_strict) {
                    if (is_upper) ++jstart; else --jend;
                }
                for (size_type j = jstart; j < jend; ++j) {
                    if constexpr (is_unit) {
                        y[i] += (i == j ? T(1) : mat_(i, j)) * x[j];
                    } else {
                        y[i] += mat_(i, j) * x[j];
                    }
                }
            }
            return y;
        }

        /// 三角部分のみの密行列を生成
        Matrix<T> toDense() const {
            const auto n = rows();
            Matrix<T> result(n, n, T(0));
            for (size_type i = 0; i < n; ++i)
                for (size_type j = 0; j < n; ++j)
                    result(i, j) = (*this)(i, j);
            return result;
        }

        /// 元の行列に三角部分を書き戻す (対称行列の生成などに使用)
        template<typename SrcMatrix>
        void setFrom(const SrcMatrix& src) requires (!std::is_const_v<MatrixType>) {
            const auto n = rows();
            for (size_type i = 0; i < n; ++i) {
                size_type jstart = is_upper ? i : 0;
                size_type jend   = is_upper ? n : i + 1;
                for (size_type j = jstart; j < jend; ++j)
                    mat_(i, j) = src(i, j);
            }
        }

    private:
        MatrixType& mat_;
    };

    /**
     * @brief SelfAdjointView — 対称行列ビュー (Eigen 互換)
     *
     * 上三角または下三角から対称行列全体を参照する。
     * 格納されていない側の要素は転置で返す。
     */
    template<typename MatrixType, TriangularMode Mode>
    class SelfAdjointView {
        static_assert(Mode == TriangularMode::Upper || Mode == TriangularMode::Lower,
                      "SelfAdjointView requires Upper or Lower mode");
    public:
        using T = typename std::remove_cvref_t<MatrixType>::value_type;
        using size_type = std::size_t;

        static constexpr bool from_upper = (Mode == TriangularMode::Upper);

        explicit SelfAdjointView(MatrixType& mat) : mat_(mat) {
            if (mat_.rows() != mat_.cols())
                throw DimensionError("SelfAdjointView: matrix must be square");
        }

        size_type rows() const { return mat_.rows(); }
        size_type cols() const { return mat_.cols(); }

        /// 要素アクセス (格納されていない側は転置)
        T operator()(size_type i, size_type j) const {
            if constexpr (from_upper) {
                return (i <= j) ? mat_(i, j) : mat_(j, i);
            } else {
                return (i >= j) ? mat_(i, j) : mat_(j, i);
            }
        }

        /// 対称行列ベクトル積: y = this * x
        /// AVX2 FMA による内積 + axpy の同時実行
        Vector<T> operator*(const Vector<T>& x) const {
            const auto n = rows();
            if (x.size() != n)
                throw DimensionError("SelfAdjointView::operator*: dimension mismatch");
            Vector<T> y(n, T(0));
            const T* __restrict A = mat_.data();
            const T* __restrict xp = x.data();
            T* __restrict yp = y.data();
            using PT = PacketTraits<T>;
            constexpr size_type W = PT::size;

            if constexpr (from_upper) {
                for (size_type i = 0; i < n; ++i) {
                    const T* __restrict row = A + i * n;
                    T xi = xp[i];
                    size_type j = i + 1;
                    size_type len = n - j;

                    if constexpr (W > 1) {
                        auto vsum = PT::set1(T(0));
                        auto vxi = PT::set1(xi);
                        size_type endw = j + (len & ~(W - 1));
                        for (; j < endw; j += W) {
                            auto va = PT::load(row + j);
                            vsum = PT::fmadd(va, PT::load(xp + j), vsum);
                            PT::store(yp + j, PT::fmadd(va, vxi, PT::load(yp + j)));
                        }
                        T sum = row[i] * xi + PT::reduce_add(vsum);
                        for (; j < n; ++j) {
                            T aij = row[j];
                            sum += aij * xp[j];
                            yp[j] += aij * xi;
                        }
                        yp[i] += sum;
                    } else {
                        T sum = row[i] * xi;
                        for (; j < n; ++j) {
                            T aij = row[j];
                            sum += aij * xp[j];
                            yp[j] += aij * xi;
                        }
                        yp[i] += sum;
                    }
                }
            } else {
                for (size_type i = 0; i < n; ++i) {
                    const T* __restrict row = A + i * n;
                    T xi = xp[i];
                    size_type j = 0;

                    if constexpr (W > 1) {
                        auto vsum = PT::set1(T(0));
                        auto vxi = PT::set1(xi);
                        size_type endw = (i / W) * W;
                        for (; j < endw; j += W) {
                            auto va = PT::load(row + j);
                            vsum = PT::fmadd(va, PT::load(xp + j), vsum);
                            PT::store(yp + j, PT::fmadd(va, vxi, PT::load(yp + j)));
                        }
                        T sum = row[i] * xi + PT::reduce_add(vsum);
                        for (; j < i; ++j) {
                            T aij = row[j];
                            sum += aij * xp[j];
                            yp[j] += aij * xi;
                        }
                        yp[i] += sum;
                    } else {
                        T sum = row[i] * xi;
                        for (; j < i; ++j) {
                            T aij = row[j];
                            sum += aij * xp[j];
                            yp[j] += aij * xi;
                        }
                        yp[i] += sum;
                    }
                }
            }
            return y;
        }

        /// 対称行列の密行列コピーを生成
        Matrix<T> toDense() const {
            const auto n = rows();
            Matrix<T> result(n, n);
            for (size_type i = 0; i < n; ++i)
                for (size_type j = 0; j < n; ++j)
                    result(i, j) = (*this)(i, j);
            return result;
        }

    private:
        MatrixType& mat_;
    };

    //-----------------------------------------------------------------------------
    // 行列操作関数
    //-----------------------------------------------------------------------------

    // 行列加算: C = A + B
    template<typename MatrixA, typename MatrixB, typename MatrixC,
        typename ComputePolicy = DefaultComputePolicy<typename MatrixA::value_type>>
        void add(const MatrixA& a, const MatrixB& b, MatrixC& c) {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        // 次元チェック
        if (a.rows() != b.rows() || a.cols() != b.cols() ||
            a.rows() != c.rows() || a.cols() != c.cols()) {
            throw DimensionError("Matrix addition: dimension mismatch");
        }

        // コンピュートポリシーに処理を委譲
        ComputePolicy::add(a, b, c);
    }

    // 行列減算: C = A - B
    template<typename MatrixA, typename MatrixB, typename MatrixC,
        typename ComputePolicy = DefaultComputePolicy<typename MatrixA::value_type>>
        void subtract(const MatrixA& a, const MatrixB& b, MatrixC& c) {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        // 次元チェック
        if (a.rows() != b.rows() || a.cols() != b.cols() ||
            a.rows() != c.rows() || a.cols() != c.cols()) {
            throw DimensionError("Matrix subtraction: dimension mismatch");
        }

        const size_type rows = a.rows();
        const size_type cols = a.cols();

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

    // 行列スカラー乗算: B = alpha * A
    template<typename MatrixA, typename Scalar, typename MatrixB,
        typename ComputePolicy = DefaultComputePolicy<typename MatrixA::value_type>>
        void scale(const MatrixA& a, const Scalar& alpha, MatrixB& b) {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        // 次元チェック
        if (a.rows() != b.rows() || a.cols() != b.cols()) {
            throw DimensionError("Matrix scaling: dimension mismatch");
        }

        const size_type rows = a.rows();
        const size_type cols = a.cols();

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

    // 行列-行列乗算: C = A * B
    template<typename MatrixA, typename MatrixB, typename MatrixC,
        typename ComputePolicy = DefaultComputePolicy<typename MatrixA::value_type>>
        void multiply(const MatrixA& a, const MatrixB& b, MatrixC& c) {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        // 次元チェック
        if (a.cols() != b.rows() ||
            c.rows() != a.rows() || c.cols() != b.cols()) {
            throw DimensionError("Matrix multiplication: dimension mismatch");
        }

        // コンピュートポリシーに処理を委譲
        ComputePolicy::multiply(a, b, c);
    }

    // 行列-ベクトル乗算: y = A * x
    // 小行列 (m<=256): 2行同時処理で x ベクトルのキャッシュ再利用
    // 大行列: 4-way アンロール FMA
    template<typename MatrixA, typename VectorX, typename VectorY,
        typename ComputePolicy = DefaultComputePolicy<typename MatrixA::value_type>>
        void multiply_vector(const MatrixA& a, const VectorX& x, VectorY& y) {
        using T = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        if (a.cols() != x.size() || y.size() != a.rows()) {
            throw DimensionError("Matrix-vector multiplication: dimension mismatch");
        }

        const size_type m = a.rows();
        const size_type n = a.cols();
        const T* a_data = a.data();
        const T* x_data = x.data();
        const size_type lda = a.cols();

        if constexpr (has_simd_packet_v<T>) {
            using PT = PacketTraits<T>;
            constexpr size_type PS = PT::size;
            constexpr size_type STRIDE = PS * 4;
            const size_type unroll_end = n - (n % STRIDE);
            const size_type vec_end = n - (n % PS);

            // 小行列パス: 2行同時処理 (x ベクトルの L1 キャッシュ再利用)
            size_type i = 0;
            if (m <= 256) {
                for (; i + 1 < m; i += 2) {
                    const T* row0 = &a_data[i * lda];
                    const T* row1 = &a_data[(i + 1) * lda];
                    auto a0 = PT::set1(T{0}), a1 = a0, b0 = a0, b1 = a0;
                    auto c0 = a0, c1 = a0, d0 = a0, d1 = a0;
                    size_type j = 0;
                    for (; j < unroll_end; j += STRIDE) {
                        auto x0 = PT::load(&x_data[j]);
                        auto x1 = PT::load(&x_data[j + PS]);
                        auto x2 = PT::load(&x_data[j + PS * 2]);
                        auto x3 = PT::load(&x_data[j + PS * 3]);
                        a0 = PT::fmadd(PT::load(&row0[j]),          x0, a0);
                        a1 = PT::fmadd(PT::load(&row0[j + PS]),     x1, a1);
                        b0 = PT::fmadd(PT::load(&row0[j + PS * 2]), x2, b0);
                        b1 = PT::fmadd(PT::load(&row0[j + PS * 3]), x3, b1);
                        c0 = PT::fmadd(PT::load(&row1[j]),          x0, c0);
                        c1 = PT::fmadd(PT::load(&row1[j + PS]),     x1, c1);
                        d0 = PT::fmadd(PT::load(&row1[j + PS * 2]), x2, d0);
                        d1 = PT::fmadd(PT::load(&row1[j + PS * 3]), x3, d1);
                    }
                    auto sum0_p = PT::add(PT::add(a0, a1), PT::add(b0, b1));
                    auto sum1_p = PT::add(PT::add(c0, c1), PT::add(d0, d1));
                    // 端数をパケットアキュムレータに追加 (reduce_add は最後に 1 回だけ)
                    for (; j < vec_end; j += PS) {
                        auto xp = PT::load(&x_data[j]);
                        sum0_p = PT::fmadd(PT::load(&row0[j]), xp, sum0_p);
                        sum1_p = PT::fmadd(PT::load(&row1[j]), xp, sum1_p);
                    }
                    T s0 = PT::reduce_add(sum0_p);
                    T s1 = PT::reduce_add(sum1_p);
                    for (; j < n; ++j) {
                        s0 += row0[j] * x_data[j];
                        s1 += row1[j] * x_data[j];
                    }
                    y[i] = s0;
                    y[i + 1] = s1;
                }
            }

            // 残り (奇数行 or 大行列)
            for (; i < m; ++i) {
                const T* a_row = &a_data[i * lda];
                auto acc0 = PT::set1(T{0}), acc1 = acc0, acc2 = acc0, acc3 = acc0;
                size_type j = 0;
                for (; j < unroll_end; j += STRIDE) {
                    acc0 = PT::fmadd(PT::load(&a_row[j]),          PT::load(&x_data[j]),          acc0);
                    acc1 = PT::fmadd(PT::load(&a_row[j + PS]),     PT::load(&x_data[j + PS]),     acc1);
                    acc2 = PT::fmadd(PT::load(&a_row[j + PS * 2]), PT::load(&x_data[j + PS * 2]), acc2);
                    acc3 = PT::fmadd(PT::load(&a_row[j + PS * 3]), PT::load(&x_data[j + PS * 3]), acc3);
                }
                auto acc = PT::add(PT::add(acc0, acc1), PT::add(acc2, acc3));
                // 端数もパケットアキュムレータに追加
                for (; j < vec_end; j += PS)
                    acc = PT::fmadd(PT::load(&a_row[j]), PT::load(&x_data[j]), acc);
                T sum = PT::reduce_add(acc);
                for (; j < n; ++j)
                    sum += a_row[j] * x_data[j];
                y[i] = sum;
            }
        } else {
            for (size_type i = 0; i < m; ++i) {
                T sum = T{0};
                for (size_type j = 0; j < n; ++j)
                    sum += a(i, j) * x[j];
                y[i] = sum;
            }
        }
    }

    // 行列転置: B = A^T (computation_policy の最適化転置に委譲)
    template<typename MatrixA, typename MatrixB,
        typename ComputePolicy = DefaultComputePolicy<typename MatrixA::value_type>>
        void transpose(const MatrixA& a, MatrixB& b) {
        ComputePolicy::transpose(a, b);
    }

    // 行列の対角要素抽出: diag = diag(A)
    template<typename MatrixA, typename VectorDiag>
    void extract_diagonal(const MatrixA& a, VectorDiag& diag) {
        using size_type = typename MatrixA::size_type;

        const size_type min_dim = std::min(a.rows(), a.cols());

        // サイズ調整
        if (diag.size() != min_dim) {
            diag.resize(min_dim);
        }

        // 対角要素の抽出
        for (size_type i = 0; i < min_dim; ++i) {
            diag[i] = a(i, i);
        }
    }

    // 対角行列の生成: A = diag(v)
    template<typename VectorDiag, typename MatrixA>
    void diagonal_matrix(const VectorDiag& diag, MatrixA& a) {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        const size_type n = diag.size();

        // 次元チェック
        if (a.rows() != n || a.cols() != n) {
            throw DimensionError("diagonal_matrix: dimension mismatch");
        }

        // ゼロ初期化
        a.zero();

        // 対角要素の設定
        for (size_type i = 0; i < n; ++i) {
            a(i, i) = diag[i];
        }
    }

    // トレース計算: tr(A)
    template<typename MatrixA>
    auto trace(const MatrixA& a) -> typename MatrixA::value_type {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        if (!a.is_square()) {
            throw DimensionError("trace: matrix must be square");
        }

        const size_type n = a.rows();
        value_type result = numeric_traits<value_type>::zero();

        for (size_type i = 0; i < n; ++i) {
            result += a(i, i);
        }

        return result;
    }

    // フロベニウスノルム: ||A||_F
    template<typename MatrixA>
    auto frobenius_norm(const MatrixA& a) -> typename MatrixA::value_type {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        const size_type rows = a.rows();
        const size_type cols = a.cols();

        value_type sum_squares = numeric_traits<value_type>::zero();

        for (size_type i = 0; i < rows; ++i) {
            for (size_type j = 0; j < cols; ++j) {
                sum_squares += a(i, j) * a(i, j);
            }
        }

        return static_cast<value_type>(std::sqrt(static_cast<double>(sum_squares)));
    }

    // L1 行列ノルム (最大列和): max_j Σ_i |a(i,j)|
    template<typename MatrixA>
    auto matrix_norm_l1(const MatrixA& a) -> typename MatrixA::value_type {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        const size_type rows = a.rows();
        const size_type cols = a.cols();

        value_type max_col_sum = numeric_traits<value_type>::zero();

        for (size_type j = 0; j < cols; ++j) {
            value_type col_sum = numeric_traits<value_type>::zero();
            for (size_type i = 0; i < rows; ++i) {
                col_sum += static_cast<value_type>(
                    std::abs(static_cast<double>(a(i, j))));
            }
            if (col_sum > max_col_sum) max_col_sum = col_sum;
        }

        return max_col_sum;
    }

    // L∞ 行列ノルム (最大行和): max_i Σ_j |a(i,j)|
    template<typename MatrixA>
    auto matrix_norm_linf(const MatrixA& a) -> typename MatrixA::value_type {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        const size_type rows = a.rows();
        const size_type cols = a.cols();

        value_type max_row_sum = numeric_traits<value_type>::zero();

        for (size_type i = 0; i < rows; ++i) {
            value_type row_sum = numeric_traits<value_type>::zero();
            for (size_type j = 0; j < cols; ++j) {
                row_sum += static_cast<value_type>(
                    std::abs(static_cast<double>(a(i, j))));
            }
            if (row_sum > max_row_sum) max_row_sum = row_sum;
        }

        return max_row_sum;
    }

    // 行列式の計算（基本的な実装、小サイズ行列用）
    template<typename MatrixA>
    auto determinant(const MatrixA& a) -> typename MatrixA::value_type {
        using value_type = typename MatrixA::value_type;
        using size_type = typename MatrixA::size_type;

        if (!a.is_square()) {
            throw DimensionError("determinant: matrix must be square");
        }

        const size_type n = a.rows();

        // 1x1 行列
        if (n == 1) {
            return a(0, 0);
        }

        // 2x2 行列
        if (n == 2) {
            return a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
        }

        // 3x3 行列
        if (n == 3) {
            return a(0, 0) * (a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1)) -
                a(0, 1) * (a(1, 0) * a(2, 2) - a(1, 2) * a(2, 0)) +
                a(0, 2) * (a(1, 0) * a(2, 1) - a(1, 1) * a(2, 0));
        }

        // 4x4以上: 型に応じてアルゴリズムを選択
        if constexpr (numeric_traits<value_type>::is_integer) {
            // 整数型: Bareiss アルゴリズム (除算なしで行列式を計算)
            // M[i][j] = (M[i][j]*M[k][k] - M[i][k]*M[k][j]) / prev_pivot
            // この除算は常に割り切れることが保証される (Sylvester の恒等式)
            Matrix<value_type> M(n, n);
            for (size_type i = 0; i < n; ++i)
                for (size_type j = 0; j < n; ++j)
                    M(i, j) = a(i, j);
            int sign = 1;
            value_type prev_pivot = numeric_traits<value_type>::one();

            for (size_type k = 0; k < n; ++k) {
                // 部分ピボット選択
                size_type pivot_row = k;
                for (size_type i = k + 1; i < n; ++i) {
                    if (M(i, k) != numeric_traits<value_type>::zero() &&
                        (M(pivot_row, k) == numeric_traits<value_type>::zero() ||
                         numeric_traits<value_type>::pivotBetter(M(i, k), M(pivot_row, k))))
                        pivot_row = i;
                }

                if (pivot_row != k) {
                    for (size_type j = 0; j < n; ++j)
                        std::swap(M(k, j), M(pivot_row, j));
                    sign = -sign;
                }

                if (M(k, k) == numeric_traits<value_type>::zero())
                    return numeric_traits<value_type>::zero();

                for (size_type i = k + 1; i < n; ++i) {
                    for (size_type j = k + 1; j < n; ++j) {
                        M(i, j) = (M(i, j) * M(k, k) - M(i, k) * M(k, j)) / prev_pivot;
                    }
                }
                // i > k の行の k 列目をゼロにする (次のステップで不要だが明示的に)
                for (size_type i = k + 1; i < n; ++i)
                    M(i, k) = numeric_traits<value_type>::zero();

                prev_pivot = M(k, k);
            }
            return M(n - 1, n - 1) * value_type(sign);
        } else {
            // 浮動小数点・有理数: ガウス消去法で O(n³) に行列式を計算
            Matrix<value_type> lu(n, n);
            for (size_type i = 0; i < n; ++i)
                for (size_type j = 0; j < n; ++j)
                    lu(i, j) = a(i, j);
            value_type det = numeric_traits<value_type>::one();
            int sign = 1;

            for (size_type k = 0; k < n; ++k) {
                // 部分ピボット選択 (pivotBetter で型に応じた最適な基準を使用)
                size_type pivot_row = k;
                for (size_type i = k + 1; i < n; ++i) {
                    if (lu(i, k) != numeric_traits<value_type>::zero() &&
                        (lu(pivot_row, k) == numeric_traits<value_type>::zero() ||
                         numeric_traits<value_type>::pivotBetter(lu(i, k), lu(pivot_row, k))))
                        pivot_row = i;
                }

                if (pivot_row != k) {
                    for (size_type j = 0; j < n; ++j)
                        std::swap(lu(k, j), lu(pivot_row, j));
                    sign = -sign;
                }

                if (lu(k, k) == numeric_traits<value_type>::zero())
                    return numeric_traits<value_type>::zero(); // 特異

                for (size_type i = k + 1; i < n; ++i) {
                    value_type factor = lu(i, k) / lu(k, k);
                    for (size_type j = k + 1; j < n; ++j)
                        lu(i, j) -= factor * lu(k, j);
                }
            }

            for (size_type i = 0; i < n; ++i)
                det *= lu(i, i);
            if (sign < 0) det = -det;
            return det;
        }
    }

    // 逆行列の計算（Gauss-Jordan 消去法）
    // 除算が必要なため、体 (field) の要素型でのみ使用可能
    template<typename T>
    Matrix<T> inverse(const Matrix<T>& a) {
        using size_type = typename Matrix<T>::size_type;

        if (!a.is_square()) {
            throw DimensionError("inverse: matrix must be square");
        }

        const size_type n = a.rows();

        if (n == 0) {
            throw DimensionError("inverse: matrix is empty");
        }

        // 1x1 行列
        if (n == 1) {
            T val = a(0, 0);
            if (val == numeric_traits<T>::zero()) {
                throw std::runtime_error("inverse: singular matrix");
            }
            Matrix<T> result(1, 1);
            result(0, 0) = numeric_traits<T>::one() / val;
            return result;
        }

        // 2x2 行列 (直接公式)
        if (n == 2) {
            T det = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
            if (det == numeric_traits<T>::zero()) {
                throw std::runtime_error("inverse: singular matrix");
            }
            T inv_det = numeric_traits<T>::one() / det;
            Matrix<T> result(2, 2);
            result(0, 0) = a(1, 1) * inv_det;
            result(0, 1) = -a(0, 1) * inv_det;
            result(1, 0) = -a(1, 0) * inv_det;
            result(1, 1) = a(0, 0) * inv_det;
            return result;
        }

        // 一般: Gauss-Jordan 消去法
        // [A | I] → [I | A^{-1}]
        Matrix<T> aug(n, 2 * n);
        for (size_type i = 0; i < n; ++i) {
            for (size_type j = 0; j < n; ++j) {
                aug(i, j) = a(i, j);
            }
            aug(i, n + i) = numeric_traits<T>::one();
        }

        for (size_type col = 0; col < n; ++col) {
            // 部分ピボット選択 (pivotBetter で型に応じた最適な基準を使用)
            size_type pivot_row = col;
            for (size_type i = col + 1; i < n; ++i) {
                if (!(aug(i, col) == numeric_traits<T>::zero()) &&
                    (aug(pivot_row, col) == numeric_traits<T>::zero() ||
                     numeric_traits<T>::pivotBetter(aug(i, col), aug(pivot_row, col))))
                    pivot_row = i;
            }

            // ピボット行の交換
            if (pivot_row != col) {
                for (size_type j = 0; j < 2 * n; ++j) {
                    T tmp = aug(col, j);
                    aug(col, j) = aug(pivot_row, j);
                    aug(pivot_row, j) = tmp;
                }
            }

            T pivot = aug(col, col);
            if (pivot == numeric_traits<T>::zero()) {
                throw std::runtime_error("inverse: singular matrix");
            }

            // ピボット行をスケーリング
            T inv_pivot = numeric_traits<T>::one() / pivot;
            for (size_type j = 0; j < 2 * n; ++j) {
                aug(col, j) *= inv_pivot;
            }

            // 他の行から消去
            for (size_type i = 0; i < n; ++i) {
                if (i != col) {
                    T factor = aug(i, col);
                    for (size_type j = 0; j < 2 * n; ++j) {
                        aug(i, j) -= factor * aug(col, j);
                    }
                }
            }
        }

        // 結果を抽出
        Matrix<T> result(n, n);
        for (size_type i = 0; i < n; ++i) {
            for (size_type j = 0; j < n; ++j) {
                result(i, j) = aug(i, n + j);
            }
        }
        return result;
    }

    //-----------------------------------------------------------------------------
    // 主要な行列クラス
    //-----------------------------------------------------------------------------

    // 動的サイズ行列
    template<typename T>
    class Matrix : public detail::MatrixBase<T, detail::DynamicMatrixStorage<T>>,
                   public MatExpr<Matrix<T>> {
    public:
        using Base = detail::MatrixBase<T, detail::DynamicMatrixStorage<T>>;
        using value_type = typename Base::value_type;
        using size_type = typename Base::size_type;
        using reference = typename Base::reference;
        using const_reference = typename Base::const_reference;
        using iterator = typename Base::iterator;
        using const_iterator = typename Base::const_iterator;

        // コンストラクタ
        Matrix() : Base() {}

        // 行数と列数を指定するコンストラクタ
        Matrix(size_type rows, size_type cols) : Base(rows, cols) {}

        // 行数と列数と初期値を指定するコンストラクタ
        Matrix(size_type rows, size_type cols, const T& value) : Base(rows, cols, value) {}

        // 初期化リストを使用するコンストラクタ (2次元)
        Matrix(std::initializer_list<std::initializer_list<T>> init) : Base() {
            size_type rows = init.size();
            size_type cols = rows > 0 ? init.begin()->size() : 0;

            this->resize(rows, cols);

            size_type i = 0;
            for (const auto& row : init) {
                size_type j = 0;
                for (const auto& val : row) {
                    if (j < cols) {
                        (*this)(i, j) = val;
                        ++j;
                    }
                }
                ++i;
            }
        }

        // コピー/ムーブコンストラクタとアサインメント
        Matrix(const Matrix& other) = default;
        Matrix(Matrix&& other) noexcept = default;
        Matrix& operator=(const Matrix& other) = default;
        Matrix& operator=(Matrix&& other) noexcept = default;

        // SIMD パケットロード (線形インデックス, 行優先連続メモリ)
        auto packet(std::size_t idx) const {
            return PacketTraits<T>::load(this->data() + idx);
        }

        // 式テンプレートからの構築
        template<typename E>
        Matrix(const MatExpr<E>& expr) : Base(expr.derived().rows(), expr.derived().cols()) {
            assignFromMatExpr(expr.derived());
        }

        // 式テンプレートからの代入
        template<typename E>
        Matrix& operator=(const MatExpr<E>& expr) {
            const E& e = expr.derived();
            this->resize(e.rows(), e.cols());
            assignFromMatExpr(e);
            return *this;
        }

    private:
        template<typename E>
        void assignFromMatExpr(const E& e) {
            const size_type total = this->rows() * this->cols();
            if constexpr (has_simd_packet_v<T>) {
                using PT = PacketTraits<T>;
                constexpr size_type PS = PT::size;
                constexpr size_type PS4 = PS * 4;
                T* __restrict dst = this->data();
                const size_type unroll_end = total - (total % PS4);
                const size_type vec_end = total - (total % PS);
                size_type i = 0;
                for (; i < unroll_end; i += PS4) {
                    PT::store(dst + i,          e.packet(i));
                    PT::store(dst + i + PS,     e.packet(i + PS));
                    PT::store(dst + i + PS * 2, e.packet(i + PS * 2));
                    PT::store(dst + i + PS * 3, e.packet(i + PS * 3));
                }
                for (; i < vec_end; i += PS)
                    PT::store(dst + i, e.packet(i));
                // 残端はスカラー (線形→(row,col)変換)
                for (; i < total; ++i)
                    dst[i] = static_cast<T>(e(i / this->cols(), i % this->cols()));
            } else {
                for (size_type i = 0; i < this->rows(); ++i)
                    for (size_type j = 0; j < this->cols(); ++j)
                        (*this)(i, j) = static_cast<T>(e(i, j));
            }
        }

    public:

        // 他のMatrixからのコピー構築
        template<typename OtherT, typename OtherStorage>
        Matrix(const detail::MatrixBase<OtherT, OtherStorage>& other)
            : Base(other.rows(), other.cols()) {

            for (size_type i = 0; i < this->rows(); ++i) {
                for (size_type j = 0; j < this->cols(); ++j) {
                    (*this)(i, j) = static_cast<T>(other(i, j));
                }
            }
        }

        // MatrixViewからのコピー構築
        template<typename MatrixType>
        Matrix(const detail::MatrixView<MatrixType>& view)
            : Base(view.rows(), view.cols()) {

            for (size_type i = 0; i < this->rows(); ++i) {
                for (size_type j = 0; j < this->cols(); ++j) {
                    (*this)(i, j) = static_cast<T>(view(i, j));
                }
            }
        }

        // ConstMatrixViewからのコピー構築
        template<typename MatrixType>
        Matrix(const detail::ConstMatrixView<MatrixType>& view)
            : Base(view.rows(), view.cols()) {

            for (size_type i = 0; i < this->rows(); ++i) {
                for (size_type j = 0; j < this->cols(); ++j) {
                    (*this)(i, j) = static_cast<T>(view(i, j));
                }
            }
        }

        // デストラクタ
        ~Matrix() = default;

        /// 三角行列ビュー
        template<TriangularMode Mode>
        TriangularView<Matrix, Mode> triangularView() {
            return TriangularView<Matrix, Mode>(*this);
        }
        template<TriangularMode Mode>
        TriangularView<const Matrix, Mode> triangularView() const {
            return TriangularView<const Matrix, Mode>(const_cast<const Matrix&>(*this));
        }

        /// 対称行列ビュー
        template<TriangularMode Mode>
        SelfAdjointView<Matrix, Mode> selfAdjointView() {
            return SelfAdjointView<Matrix, Mode>(*this);
        }
        template<TriangularMode Mode>
        SelfAdjointView<const Matrix, Mode> selfAdjointView() const {
            return SelfAdjointView<const Matrix, Mode>(const_cast<const Matrix&>(*this));
        }

        // 加算演算子
        Matrix& operator+=(const Matrix& rhs) {
            if (this->rows() != rhs.rows() || this->cols() != rhs.cols()) {
                throw DimensionError("Matrix addition: dimension mismatch");
            }

            for (size_type i = 0; i < this->rows(); ++i) {
                for (size_type j = 0; j < this->cols(); ++j) {
                    (*this)(i, j) += rhs(i, j);
                }
            }

            return *this;
        }

        // 減算演算子
        Matrix& operator-=(const Matrix& rhs) {
            if (this->rows() != rhs.rows() || this->cols() != rhs.cols()) {
                throw DimensionError("Matrix subtraction: dimension mismatch");
            }

            for (size_type i = 0; i < this->rows(); ++i) {
                for (size_type j = 0; j < this->cols(); ++j) {
                    (*this)(i, j) -= rhs(i, j);
                }
            }

            return *this;
        }

        // スカラー乗算演算子
        Matrix& operator*=(const T& scalar) {
            for (size_type i = 0; i < this->rows(); ++i) {
                for (size_type j = 0; j < this->cols(); ++j) {
                    (*this)(i, j) *= scalar;
                }
            }

            return *this;
        }

        // スカラー除算演算子
        Matrix& operator/=(const T& scalar) {
            // ゼロ除算チェック
            if (scalar == T{ 0 }) {
                throw std::invalid_argument("Division by zero");
            }

            for (size_type i = 0; i < this->rows(); ++i) {
                for (size_type j = 0; j < this->cols(); ++j) {
                    (*this)(i, j) /= scalar;
                }
            }

            return *this;
        }

        // 式テンプレートからの複合代入
        template<typename E>
        Matrix& operator+=(const MatExpr<E>& rhs) {
            const E& r = rhs.derived();
            for (size_type i = 0; i < this->rows(); ++i)
                for (size_type j = 0; j < this->cols(); ++j)
                    (*this)(i, j) += static_cast<T>(r(i, j));
            return *this;
        }

        template<typename E>
        Matrix& operator-=(const MatExpr<E>& rhs) {
            const E& r = rhs.derived();
            for (size_type i = 0; i < this->rows(); ++i)
                for (size_type j = 0; j < this->cols(); ++j)
                    (*this)(i, j) -= static_cast<T>(r(i, j));
            return *this;
        }

        // 転置行列の取得 (新規 Matrix を確保して返す)
        // 大行列ではヒープ確保+ゼロ初期化のオーバーヘッドあり
        // 高性能が必要なら transpose_into() で既存バッファに書き込む
        Matrix transpose() const {
            Matrix result(this->cols(), this->rows());
            DefaultComputePolicy<T>::transpose(*this, result);
            return result;
        }

        // 既存バッファへの転置 (allocation-free)
        // dst は (cols × rows) サイズに事前確保済みであること
        void transpose_into(Matrix& dst) const {
            DefaultComputePolicy<T>::transpose(*this, dst);
        }

        // 単位行列の生成（静的メソッド）
        static Matrix identity(size_type size) {
            Matrix result(size, size);

            // ゼロで初期化
            result.zero();

            // 対角成分を1に設定
            for (size_type i = 0; i < size; ++i) {
                result(i, i) = numeric_traits<T>::one();
            }

            return result;
        }

        // 行列式の計算
        T determinant() const {
            return calx::determinant(*this);
        }

        // フロベニウスノルムの計算
        T norm() const {
            return calx::frobenius_norm(*this);
        }

        // L1 行列ノルム (最大列和)
        T norm_l1() const {
            return calx::matrix_norm_l1(*this);
        }

        // L∞ 行列ノルム (最大行和)
        T norm_linf() const {
            return calx::matrix_norm_linf(*this);
        }

        // トレースの計算
        T trace() const {
            return calx::trace(*this);
        }

        // 逆行列の計算
        [[nodiscard]] Matrix inverse() const {
            return calx::inverse(*this);
        }

        /// operator^: A^'T' = 転置, A^'H' = エルミート転置, A^(-1) = 逆行列, A^n = 冪乗
        ///
        /// 用例:
        ///   auto At = A ^ 'T';    // 転置
        ///   auto Ah = A ^ 'H';    // エルミート転置 (共役転置)
        ///   auto Ainv = A ^ (-1); // 逆行列
        ///   auto A3 = A ^ 3;      // A*A*A
        [[nodiscard]] Matrix operator^(int n) const {
            if (n == 'T') return this->transpose();
            if (n == 'H') {
                // エルミート転置: conj(A^T)
                Matrix result = this->transpose();
                for (size_type i = 0; i < result.rows(); ++i)
                    for (size_type j = 0; j < result.cols(); ++j)
                        result(i, j) = numeric_traits<T>::conj(result(i, j));
                return result;
            }
            if (n == -1) return this->inverse();
            if (n < 0)
                return calx::power(this->inverse(), static_cast<unsigned int>(-n));
            return calx::power(*this, static_cast<unsigned int>(n));
        }
   };

   // 固定サイズ行列
   template<typename T, std::size_t Rows, std::size_t Cols>
   class StaticMatrix : public detail::MatrixBase<T, detail::StaticMatrixStorage<T, Rows, Cols>>,
                        public MatExpr<StaticMatrix<T, Rows, Cols>> {
   public:
       using Base = detail::MatrixBase<T, detail::StaticMatrixStorage<T, Rows, Cols>>;
       using value_type = typename Base::value_type;
       using size_type = typename Base::size_type;
       using reference = typename Base::reference;
       using const_reference = typename Base::const_reference;
       using iterator = typename Base::iterator;
       using const_iterator = typename Base::const_iterator;

       // コンストラクタ
       StaticMatrix() : Base() {}

       // 初期値を指定するコンストラクタ
       explicit StaticMatrix(const T& value) {
           this->fill(value);
       }

       // 初期化リストを使用するコンストラクタ (2次元)
       StaticMatrix(std::initializer_list<std::initializer_list<T>> init) : Base() {
           size_type i = 0;
           for (const auto& row : init) {
               if (i >= Rows) break;

               size_type j = 0;
               for (const auto& val : row) {
                   if (j >= Cols) break;
                   (*this)(i, j) = val;
                   ++j;
               }
               ++i;
           }
       }

       // コピー/ムーブコンストラクタとアサインメント
       StaticMatrix(const StaticMatrix&) = default;
       StaticMatrix(StaticMatrix&&) noexcept = default;
       StaticMatrix& operator=(const StaticMatrix&) = default;
       StaticMatrix& operator=(StaticMatrix&&) noexcept = default;

       // SIMD パケットロード (線形インデックス)
       auto packet(std::size_t idx) const {
           return PacketTraits<T>::load(this->data() + idx);
       }

       // 式テンプレートからの構築
       template<typename E>
       StaticMatrix(const MatExpr<E>& expr) : Base() {
           assignFromMatExpr(expr.derived());
       }

       // 式テンプレートからの代入
       template<typename E>
       StaticMatrix& operator=(const MatExpr<E>& expr) {
           assignFromMatExpr(expr.derived());
           return *this;
       }

    private:
       template<typename E>
       void assignFromMatExpr(const E& e) {
           constexpr size_type total = Rows * Cols;
           if constexpr (has_simd_packet_v<T> && total >= PacketTraits<T>::size) {
               constexpr size_type PS = PacketTraits<T>::size;
               constexpr size_type vec_end = total - (total % PS);
               T* __restrict dst = this->data();
               for (size_type i = 0; i < vec_end; i += PS)
                   PacketTraits<T>::store(dst + i, e.packet(i));
               for (size_type i = vec_end; i < total; ++i)
                   dst[i] = static_cast<T>(e(i / Cols, i % Cols));
           } else {
               for (size_type i = 0; i < Rows; ++i)
                   for (size_type j = 0; j < Cols; ++j)
                       (*this)(i, j) = static_cast<T>(e(i, j));
           }
       }

    public:

       // 他のMatrixからのコピー構築
       template<typename OtherT, typename OtherStorage>
       StaticMatrix(const detail::MatrixBase<OtherT, OtherStorage>& other) {
           if (other.rows() != Rows || other.cols() != Cols) {
               throw DimensionError("StaticMatrix: dimension mismatch in copy constructor");
           }

           for (size_type i = 0; i < Rows; ++i) {
               for (size_type j = 0; j < Cols; ++j) {
                   (*this)(i, j) = static_cast<T>(other(i, j));
               }
           }
       }

       // MatrixViewからのコピー構築
       template<typename MatrixType>
       StaticMatrix(const detail::MatrixView<MatrixType>& view) {
           if (view.rows() != Rows || view.cols() != Cols) {
               throw DimensionError("StaticMatrix: dimension mismatch in copy constructor from view");
           }

           for (size_type i = 0; i < Rows; ++i) {
               for (size_type j = 0; j < Cols; ++j) {
                   (*this)(i, j) = static_cast<T>(view(i, j));
               }
           }
       }

       // ConstMatrixViewからのコピー構築
       template<typename MatrixType>
       StaticMatrix(const detail::ConstMatrixView<MatrixType>& view) {
           if (view.rows() != Rows || view.cols() != Cols) {
               throw DimensionError("StaticMatrix: dimension mismatch in copy constructor from const view");
           }

           for (size_type i = 0; i < Rows; ++i) {
               for (size_type j = 0; j < Cols; ++j) {
                   (*this)(i, j) = static_cast<T>(view(i, j));
               }
           }
       }

       // デストラクタ
       ~StaticMatrix() = default;

       // 加算演算子
       StaticMatrix& operator+=(const StaticMatrix& rhs) {
           for (size_type i = 0; i < Rows; ++i) {
               for (size_type j = 0; j < Cols; ++j) {
                   (*this)(i, j) += rhs(i, j);
               }
           }

           return *this;
       }

       // 減算演算子
       StaticMatrix& operator-=(const StaticMatrix& rhs) {
           for (size_type i = 0; i < Rows; ++i) {
               for (size_type j = 0; j < Cols; ++j) {
                   (*this)(i, j) -= rhs(i, j);
               }
           }

           return *this;
       }

       // スカラー乗算演算子
       StaticMatrix& operator*=(const T& scalar) {
           for (size_type i = 0; i < Rows; ++i) {
               for (size_type j = 0; j < Cols; ++j) {
                   (*this)(i, j) *= scalar;
               }
           }

           return *this;
       }

       // スカラー除算演算子
       StaticMatrix& operator/=(const T& scalar) {
           // ゼロ除算チェック
           if (scalar == T{ 0 }) {
               throw std::invalid_argument("Division by zero");
           }

           for (size_type i = 0; i < Rows; ++i) {
               for (size_type j = 0; j < Cols; ++j) {
                   (*this)(i, j) /= scalar;
               }
           }

           return *this;
       }

       // 式テンプレートからの複合代入
       template<typename E>
       StaticMatrix& operator+=(const MatExpr<E>& rhs) {
           const E& r = rhs.derived();
           for (size_type i = 0; i < Rows; ++i)
               for (size_type j = 0; j < Cols; ++j)
                   (*this)(i, j) += static_cast<T>(r(i, j));
           return *this;
       }

       template<typename E>
       StaticMatrix& operator-=(const MatExpr<E>& rhs) {
           const E& r = rhs.derived();
           for (size_type i = 0; i < Rows; ++i)
               for (size_type j = 0; j < Cols; ++j)
                   (*this)(i, j) -= static_cast<T>(r(i, j));
           return *this;
       }

       // 転置行列の取得
       StaticMatrix<T, Cols, Rows> transpose() const {
           StaticMatrix<T, Cols, Rows> result;

           for (size_type i = 0; i < Rows; ++i) {
               for (size_type j = 0; j < Cols; ++j) {
                   result(j, i) = (*this)(i, j);
               }
           }

           return result;
       }

       // 正方行列用の特殊化メソッド
       T determinant() const requires (Rows == Cols) {
           return calx::determinant(*this);
       }

       T trace() const requires (Rows == Cols) {
           return calx::trace(*this);
       }

       static StaticMatrix identity() requires (Rows == Cols) {
           StaticMatrix result;
           result.zero();  // まずゼロ初期化
           for (size_type i = 0; i < Rows; ++i) {
               result(i, i) = numeric_traits<T>::one();  // 対角要素を1に
           }
           return result;
       }

       // フロベニウスノルムの計算
       T norm() const {
           return calx::frobenius_norm(*this);
       }

       // L1 行列ノルム (最大列和)
       T norm_l1() const {
           return calx::matrix_norm_l1(*this);
       }

       // L∞ 行列ノルム (最大行和)
       T norm_linf() const {
           return calx::matrix_norm_linf(*this);
       }

       // コンパイル時に行数と列数を取得
       static constexpr size_type static_rows() { return Rows; }
       static constexpr size_type static_cols() { return Cols; }
   };

   //-----------------------------------------------------------------------------
   // グローバル演算子
   //-----------------------------------------------------------------------------

   // 要素単位演算 (+, -, *scalar, /scalar, 単項-) は
   // expr_templates.hpp の MatExpr 演算子に完全委譲 (ET チェーン有効)
   // → 2.0*A + B - C/3.0 のような複合式でも 1 パスで SIMD 評価

   // 行列-行列乗算
   template<typename T>
   Matrix<T> operator*(const Matrix<T>& lhs, const Matrix<T>& rhs) {
       Matrix<T> result(lhs.rows(), rhs.cols());
       multiply(lhs, rhs, result);
       return result;
   }

   // 行列-行列乗算 (MatExpr オーバーロード: 式を評価してから乗算)
   template<typename T, typename E>
   Matrix<T> operator*(const Matrix<T>& lhs, const MatExpr<E>& rhs) {
       Matrix<T> rhs_eval = rhs;
       return lhs * rhs_eval;
   }

   template<typename T, typename E>
   Matrix<T> operator*(const MatExpr<E>& lhs, const Matrix<T>& rhs) {
       Matrix<T> lhs_eval = lhs;
       return lhs_eval * rhs;
   }

   // 行列-ベクトル乗算
   template<typename T>
   Vector<T> operator*(const Matrix<T>& mat, const Vector<T>& vec) {
       Vector<T> result(mat.rows());
       multiply_vector(mat, vec, result);
       return result;
   }

   // StaticMatrix 要素単位演算も ET に完全委譲 (MatExpr 基底経由)

   template<typename T, std::size_t M, std::size_t N, std::size_t P>
   StaticMatrix<T, M, P> operator*(const StaticMatrix<T, M, N>& lhs,
       const StaticMatrix<T, N, P>& rhs) {
       StaticMatrix<T, M, P> result;
       multiply(lhs, rhs, result);
       return result;
   }

   // 行列-ベクトル乗算（StaticMatrix * 動的 Vector）
   template<typename T, std::size_t Rows, std::size_t Cols>
   Vector<T> operator*(const StaticMatrix<T, Rows, Cols>& mat, const Vector<T>& vec) {
       Vector<T> result(Rows);
       multiply_vector(mat, vec, result);
       return result;
   }

   // 行列-ベクトル乗算（StaticMatrix * StaticVector）
   template<typename T, std::size_t Rows, std::size_t Cols>
   StaticVector<T, Rows> operator*(const StaticMatrix<T, Rows, Cols>& mat,
       const StaticVector<T, Cols>& vec) {
       StaticVector<T, Rows> result;
       for (std::size_t i = 0; i < Rows; ++i) {
           T sum = T(0);
           for (std::size_t j = 0; j < Cols; ++j) {
               sum += mat(i, j) * vec[j];
           }
           result[i] = sum;
       }
       return result;
   }

   // StaticMatrix 逆行列 (2x2 直接公式)
   template<typename T>
   [[nodiscard]] StaticMatrix<T, 2, 2> inverse(const StaticMatrix<T, 2, 2>& m) {
       T det = m(0, 0) * m(1, 1) - m(0, 1) * m(1, 0);
       if (std::abs(det) < std::numeric_limits<T>::epsilon()) {
           throw LinearAlgebraError("inverse: singular 2x2 matrix");
       }
       T inv_det = T(1) / det;
       StaticMatrix<T, 2, 2> result;
       result(0, 0) =  m(1, 1) * inv_det;
       result(0, 1) = -m(0, 1) * inv_det;
       result(1, 0) = -m(1, 0) * inv_det;
       result(1, 1) =  m(0, 0) * inv_det;
       return result;
   }

   // StaticMatrix 逆行列 (3x3 余因子行列)
   template<typename T>
   [[nodiscard]] StaticMatrix<T, 3, 3> inverse(const StaticMatrix<T, 3, 3>& m) {
       T det = m(0, 0) * (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1))
             - m(0, 1) * (m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0))
             + m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0));
       if (std::abs(det) < std::numeric_limits<T>::epsilon()) {
           throw LinearAlgebraError("inverse: singular 3x3 matrix");
       }
       T inv_det = T(1) / det;
       StaticMatrix<T, 3, 3> result;
       result(0, 0) = (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1)) * inv_det;
       result(0, 1) = (m(0, 2) * m(2, 1) - m(0, 1) * m(2, 2)) * inv_det;
       result(0, 2) = (m(0, 1) * m(1, 2) - m(0, 2) * m(1, 1)) * inv_det;
       result(1, 0) = (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) * inv_det;
       result(1, 1) = (m(0, 0) * m(2, 2) - m(0, 2) * m(2, 0)) * inv_det;
       result(1, 2) = (m(0, 2) * m(1, 0) - m(0, 0) * m(1, 2)) * inv_det;
       result(2, 0) = (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0)) * inv_det;
       result(2, 1) = (m(0, 1) * m(2, 0) - m(0, 0) * m(2, 1)) * inv_det;
       result(2, 2) = (m(0, 0) * m(1, 1) - m(0, 1) * m(1, 0)) * inv_det;
       return result;
   }

   //-----------------------------------------------------------------------------
   // ユーティリティ関数
   //-----------------------------------------------------------------------------

   // ゼロ行列の生成
   template<typename T>
   Matrix<T> zeros(std::size_t rows, std::size_t cols) {
       Matrix<T> result(rows, cols);
       result.zero();
       return result;
   }

   // 単位行列の生成
   template<typename T>
   Matrix<T> eye(std::size_t size) {
       return Matrix<T>::identity(size);
   }

   // 単位行列の生成 (identity<T>(n) 形式)
   template<typename T>
   Matrix<T> identity(std::size_t size) {
       return Matrix<T>::identity(size);
   }

   // 全ての要素が1の行列の生成
   template<typename T>
   Matrix<T> ones(std::size_t rows, std::size_t cols) {
       return Matrix<T>(rows, cols, T{ 1 });
   }

   // 特殊な2x2行列 - 回転行列
   template<typename T>
   Matrix<T> rotation2D(T angle) {
       Matrix<T> result(2, 2);

       const T cos_val = std::cos(angle);
       const T sin_val = std::sin(angle);

       result(0, 0) = cos_val;
       result(0, 1) = -sin_val;
       result(1, 0) = sin_val;
       result(1, 1) = cos_val;

       return result;
   }

   // 特殊な3x3行列 - 3Dの各軸周りの回転行列
   template<typename T>
   Matrix<T> rotationX(T angle) {
       Matrix<T> result(3, 3);
       result.zero();
       for (size_t i = 0; i < 3; ++i) {
           result(i, i) = T{ 1 };
       }

       const T cos_val = std::cos(angle);
       const T sin_val = std::sin(angle);

       result(1, 1) = cos_val;
       result(1, 2) = -sin_val;
       result(2, 1) = sin_val;
       result(2, 2) = cos_val;

       return result;
   }

   template<typename T>
   Matrix<T> rotationY(T angle) {
       Matrix<T> result(3, 3);
       result.zero();
       for (size_t i = 0; i < 3; ++i) {
           result(i, i) = T{ 1 };
       }

       const T cos_val = std::cos(angle);
       const T sin_val = std::sin(angle);

       result(0, 0) = cos_val;
       result(0, 2) = sin_val;
       result(2, 0) = -sin_val;
       result(2, 2) = cos_val;

       return result;
   }

   template<typename T>
   Matrix<T> rotationZ(T angle) {
       Matrix<T> result(3, 3);
       result.zero();
       for (size_t i = 0; i < 3; ++i) {
           result(i, i) = T{ 1 };
       }

       const T cos_val = std::cos(angle);
       const T sin_val = std::sin(angle);

       result(0, 0) = cos_val;
       result(0, 1) = -sin_val;
       result(1, 0) = sin_val;
       result(1, 1) = cos_val;

       return result;
   }

   // 行列の連結 - 横方向
   template<typename T>
   Matrix<T> hconcat(const Matrix<T>& a, const Matrix<T>& b) {
       if (a.rows() != b.rows()) {
           throw DimensionError("hconcat: matrices must have the same number of rows");
       }

       Matrix<T> result(a.rows(), a.cols() + b.cols());

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

       return result;
   }

   // 行列の連結 - 縦方向
   template<typename T>
   Matrix<T> vconcat(const Matrix<T>& a, const Matrix<T>& b) {
       if (a.cols() != b.cols()) {
           throw DimensionError("vconcat: matrices must have the same number of columns");
       }

       Matrix<T> result(a.rows() + b.rows(), a.cols());

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

       return result;
   }

   // 行列の整形 (reshape)
   template<typename T>
   Matrix<T> reshape(const Matrix<T>& mat, std::size_t new_rows, std::size_t new_cols) {
       if (mat.cols() == 0 || mat.rows() == 0) {
           if (new_rows == 0 || new_cols == 0) {
               return Matrix<T>(new_rows, new_cols);
           }
           throw DimensionError("reshape: total element count must be preserved");
       }
       if (mat.rows() * mat.cols() != new_rows * new_cols) {
           throw DimensionError("reshape: total element count must be preserved");
       }

       Matrix<T> result(new_rows, new_cols);

       for (std::size_t i = 0; i < mat.rows() * mat.cols(); ++i) {
           const std::size_t old_row = i / mat.cols();
           const std::size_t old_col = i % mat.cols();
           const std::size_t new_row = i / new_cols;
           const std::size_t new_col = i % new_cols;

           result(new_row, new_col) = mat(old_row, old_col);
       }

       return result;
   }

   //-----------------------------------------------------------------------------
   // 行列形式での入出力
   //-----------------------------------------------------------------------------

   // 行列の標準出力への出力
   template<typename T>
   void print_matrix(const Matrix<T>& mat, std::ostream& os = std::cout,
       std::size_t width = 12, std::size_t precision = 6) {
       // 以前の出力設定を保存
       std::ios::fmtflags old_flags = os.flags();
       std::streamsize old_precision = os.precision();

       // 出力形式の設定
       os << std::fixed << std::setprecision(precision);

       const std::size_t rows = mat.rows();
       const std::size_t cols = mat.cols();

       for (std::size_t i = 0; i < rows; ++i) {
           for (std::size_t j = 0; j < cols; ++j) {
               os << std::setw(width) << mat(i, j);
               if (j < cols - 1) os << " ";
           }
           os << "\n";
       }

       // 出力設定を元に戻す
       os.flags(old_flags);
       os.precision(old_precision);
   }

   // 簡単なフォーマット出力
   template<typename T>
   std::string to_string(const Matrix<T>& mat, std::size_t precision = 4) {
       std::ostringstream ss;
       print_matrix(mat, ss, 8, precision);
       return ss.str();
   }

   // 演算子<<のオーバーロード
   template<typename T>
   std::ostream& operator<<(std::ostream& os, const Matrix<T>& mat) {
       print_matrix(mat, os);
       return os;
   }

   // 行列の冪乗（右向きバイナリ法）
   // 正方行列のみ。power.hpp の汎用 power() を Matrix 向けに特殊化。
   template<typename T>
   [[nodiscard]] Matrix<T> power(const Matrix<T>& base, unsigned int exponent) {
       if (base.rows() != base.cols()) {
           throw DimensionError("power(): matrix must be square");
       }
       if (exponent == 0) return Matrix<T>::identity(base.rows());
       if (exponent == 1) return base;

       Matrix<T> y = Matrix<T>::identity(base.rows());
       unsigned int bit = std::bit_width(exponent) - 1;
       while (true) {
           y = y * y;
           if ((exponent >> bit) & 1u) {
               y = y * base;
           }
           if (bit == 0) break;
           --bit;
       }
       return y;
   }

    // =========================================================================
    // MatrixMap — 外部メモリへのゼロコピービュー (Eigen::Map 相当)
    // =========================================================================

    /**
     * @brief 外部メモリ上の行列ビュー (行優先)
     *
     * @code
     * double data[6] = {1,2,3,4,5,6};
     * auto M = calx::MatrixMap<double>(data, 2, 3);
     * M(0, 1) = 10.0;  // data[1] が変更される
     * @endcode
     */
    template<typename T>
    class MatrixMap {
    public:
        using value_type = T;
        using size_type = std::size_t;

        MatrixMap(T* data, size_type rows, size_type cols)
            : data_(data), rows_(rows), cols_(cols) {}

        size_type rows() const { return rows_; }
        size_type cols() const { return cols_; }
        T* data() { return data_; }
        const T* data() const { return data_; }

        T& operator()(size_type i, size_type j) { return data_[i * cols_ + j]; }
        const T& operator()(size_type i, size_type j) const { return data_[i * cols_ + j]; }

        /// Matrix へのコピー
        operator Matrix<T>() const {
            Matrix<T> m(rows_, cols_);
            for (size_type i = 0; i < rows_; ++i)
                for (size_type j = 0; j < cols_; ++j)
                    m(i, j) = data_[i * cols_ + j];
            return m;
        }

        /// Matrix からの代入
        MatrixMap& operator=(const Matrix<T>& m) {
            for (size_type i = 0; i < std::min(rows_, m.rows()); ++i)
                for (size_type j = 0; j < std::min(cols_, m.cols()); ++j)
                    data_[i * cols_ + j] = m(i, j);
            return *this;
        }

    private:
        T* data_;
        size_type rows_, cols_;
    };

    template<typename T>
    class ConstMatrixMap {
    public:
        using value_type = T;
        using size_type = std::size_t;

        ConstMatrixMap(const T* data, size_type rows, size_type cols)
            : data_(data), rows_(rows), cols_(cols) {}

        size_type rows() const { return rows_; }
        size_type cols() const { return cols_; }
        const T* data() const { return data_; }

        const T& operator()(size_type i, size_type j) const { return data_[i * cols_ + j]; }

        operator Matrix<T>() const {
            Matrix<T> m(rows_, cols_);
            for (size_type i = 0; i < rows_; ++i)
                for (size_type j = 0; j < cols_; ++j)
                    m(i, j) = data_[i * cols_ + j];
            return m;
        }

    private:
        const T* data_;
        size_type rows_, cols_;
    };

} // namespace calx

#endif // CALX_MATRIX_HPP