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

// vector.hpp
// ベクトルクラスの定義
//
// このファイルでは、MKL代数ライブラリのベクトルクラスを定義します。
// 固定長ベクトルや動的ベクトルの実装を含みます。

#ifndef CALX_VECTOR_HPP
#define CALX_VECTOR_HPP

#include <math/core/common.hpp>
#include <math/core/traits.hpp>
#include <math/core/expr_templates.hpp>
#include <math/core/simd_traits.hpp>
#include <math/core/vector/vector_view.hpp>

#include <algorithm>
#include <array>
#include <cassert>
#include <initializer_list>
#include <vector>
#include <cmath>

namespace calx {

    /// resize 時に値初期化をスキップするアロケータ (POD 型で未初期化を実現)
    template<typename T>
    struct DefaultInitAllocator : std::allocator<T> {
        using std::allocator<T>::allocator;
        template<typename U> struct rebind { using other = DefaultInitAllocator<U>; };

        // default-init: POD 型ではメモリに触らない
        void construct(T* p) noexcept(std::is_nothrow_default_constructible_v<T>) {
            ::new (static_cast<void*>(p)) T;
        }
        // 引数ありの場合は通常通り構築
        template<typename... Args>
        void construct(T* p, Args&&... args) {
            ::new (static_cast<void*>(p)) T(std::forward<Args>(args)...);
        }
    };

    /// 未初期化コンストラクタ用タグ型
    struct uninitialized_t {};
    /// 未初期化コンストラクタ用タグ値
    inline constexpr uninitialized_t uninitialized{};

    // basic_types.hppで前方宣言されているVectorクラスの定義
    template<typename T>
    class Vector : public VecExpr<Vector<T>> {
        using storage_type = std::vector<T, DefaultInitAllocator<T>>;
    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 storage_type::iterator;
        using const_iterator = typename storage_type::const_iterator;

        // コンストラクタ
        Vector() : data_() {}

        explicit Vector(size_type size) : data_(size) {
            std::fill(data_.begin(), data_.end(), T{});
        }

        /// 未初期化コンストラクタ: メモリ確保のみ、値は設定しない
        Vector(size_type size, uninitialized_t) : data_(size) {}

        Vector(size_type size, const T& value) : data_(size, value) {}

        Vector(std::initializer_list<T> init) : data_(init) {}

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

        // 式テンプレートからの構築
        template<typename E>
        Vector(const VecExpr<E>& expr) : data_(expr.derived().size()) {
            assignFromExpr(expr.derived());
        }

        // 式テンプレートからの代入
        template<typename E>
        Vector& operator=(const VecExpr<E>& expr) {
            const E& e = expr.derived();
            data_.resize(e.size());
            assignFromExpr(e);
            return *this;
        }

    private:
        // SIMD パケット評価による式代入の内部実装 (4x アンロール)
        template<typename E>
        void assignFromExpr(const E& e) {
            const size_type n = data_.size();
            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 = data_.data();
                const size_type unroll_end = n - (n % PS4);
                const size_type vec_end = n - (n % 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));
                for (; i < n; ++i)
                    dst[i] = static_cast<T>(e[i]);
            } else {
                for (size_type i = 0; i < n; ++i)
                    data_[i] = static_cast<T>(e[i]);
            }
        }

    public:

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

        // 要素アクセス
        reference operator[](size_type index) {
            assert(index < data_.size() && "Vector::operator[]: index out of range");
            return data_[index];
        }

        const_reference operator[](size_type index) const {
            assert(index < data_.size() && "Vector::operator[]: index out of range");
            return data_[index];
        }

        reference at(size_type index) {
            return data_.at(index);
        }

        const_reference at(size_type index) const {
            return data_.at(index);
        }

        reference front() {
            if (data_.empty()) throw IndexError("Vector::front: empty vector");
            return data_.front();
        }

        const_reference front() const {
            if (data_.empty()) throw IndexError("Vector::front: empty vector");
            return data_.front();
        }

        reference back() {
            if (data_.empty()) throw IndexError("Vector::back: empty vector");
            return data_.back();
        }

        const_reference back() const {
            if (data_.empty()) throw IndexError("Vector::back: empty vector");
            return data_.back();
        }

        // イテレータ
        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 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 clear() noexcept {
            data_.clear();
        }

        void resize(size_type count) {
            data_.resize(count);
        }

        void resize(size_type count, const value_type& value) {
            data_.resize(count, value);
        }

        // data_の assign メソッドを呼び出すためのラッパー
        void assign(size_type count, const T& value) {
            data_.assign(count, value);
        }

        // データアクセス
        T* data() noexcept {
            return data_.data();
        }

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

        // SIMD パケットロード (式テンプレートのパケット評価用)
        auto packet(std::size_t i) const {
            return PacketTraits<T>::load(&data_[i]);
        }

        // ゼロベクトル
        void zero() {
            std::fill(data_.begin(), data_.end(), numeric_traits<T>::zero());
        }

        // 演算子オーバーロード
        Vector& operator+=(const Vector& rhs) {
            if (this->size() != rhs.size()) {
                throw DimensionError("Vector addition: size mismatch");
            }

            for (size_type i = 0; i < size(); ++i) {
                data_[i] += rhs[i];
            }
            return *this;
        }

        Vector& operator-=(const Vector& rhs) {
            if (this->size() != rhs.size()) {
                throw DimensionError("Vector subtraction: size mismatch");
            }

            for (size_type i = 0; i < size(); ++i) {
                data_[i] -= rhs[i];
            }
            return *this;
        }

        Vector& operator*=(const T& scalar) {
            for (size_type i = 0; i < size(); ++i) {
                data_[i] *= scalar;
            }
            return *this;
        }

        Vector& operator/=(const T& scalar) {
            // ゼロ除算チェック
            if (scalar == T{ 0 }) {
                throw std::invalid_argument("Division by zero");
            }

            for (size_type i = 0; i < size(); ++i) {
                data_[i] /= scalar;
            }
            return *this;
        }

        // 式テンプレートからの複合代入
        template<typename E>
        Vector& operator+=(const VecExpr<E>& rhs) {
            const E& r = rhs.derived();
            const size_type n = size();
            if constexpr (has_simd_packet_v<T>) {
                constexpr size_type PS = PacketTraits<T>::size;
                const size_type vec_end = n - (n % PS);
                for (size_type i = 0; i < vec_end; i += PS) {
                    auto cur = PacketTraits<T>::load(&data_[i]);
                    PacketTraits<T>::store(&data_[i], PacketTraits<T>::add(cur, r.packet(i)));
                }
                for (size_type i = vec_end; i < n; ++i)
                    data_[i] += static_cast<T>(r[i]);
            } else {
                for (size_type i = 0; i < n; ++i)
                    data_[i] += static_cast<T>(r[i]);
            }
            return *this;
        }

        template<typename E>
        Vector& operator-=(const VecExpr<E>& rhs) {
            const E& r = rhs.derived();
            const size_type n = size();
            if constexpr (has_simd_packet_v<T>) {
                constexpr size_type PS = PacketTraits<T>::size;
                const size_type vec_end = n - (n % PS);
                for (size_type i = 0; i < vec_end; i += PS) {
                    auto cur = PacketTraits<T>::load(&data_[i]);
                    PacketTraits<T>::store(&data_[i], PacketTraits<T>::sub(cur, r.packet(i)));
                }
                for (size_type i = vec_end; i < n; ++i)
                    data_[i] -= static_cast<T>(r[i]);
            } else {
                for (size_type i = 0; i < n; ++i)
                    data_[i] -= static_cast<T>(r[i]);
            }
            return *this;
        }

        // ベクトル演算
        T dot(const Vector& rhs) const {
            if (this->size() != rhs.size()) {
                throw DimensionError("Vector dot product: size mismatch");
            }

            const size_type n = size();
            if constexpr (has_simd_packet_v<T>) {
                using PT = PacketTraits<T>;
                constexpr size_type PS = PT::size;
                constexpr size_type STRIDE = PS * 4;  // 4 アキュムレータ
                const size_type unroll_end = n - (n % STRIDE);
                const size_type vec_end = n - (n % PS);
                auto acc0 = PT::set1(T{0}), acc1 = acc0, acc2 = acc0, acc3 = acc0;
                for (size_type i = 0; i < unroll_end; i += STRIDE) {
                    acc0 = PT::fmadd(PT::load(&data_[i]),          PT::load(&rhs.data_[i]),          acc0);
                    acc1 = PT::fmadd(PT::load(&data_[i + PS]),     PT::load(&rhs.data_[i + PS]),     acc1);
                    acc2 = PT::fmadd(PT::load(&data_[i + PS * 2]), PT::load(&rhs.data_[i + PS * 2]), acc2);
                    acc3 = PT::fmadd(PT::load(&data_[i + PS * 3]), PT::load(&rhs.data_[i + PS * 3]), acc3);
                }
                acc0 = PT::add(PT::add(acc0, acc1), PT::add(acc2, acc3));
                T result = PT::reduce_add(acc0);
                for (size_type i = unroll_end; i < vec_end; i += PS)
                    result += PT::reduce_add(PT::mul(PT::load(&data_[i]), PT::load(&rhs.data_[i])));
                for (size_type i = vec_end; i < n; ++i)
                    result += data_[i] * rhs.data_[i];
                return result;
            } else {
                T result = numeric_traits<T>::zero();
                for (size_type i = 0; i < n; ++i)
                    result += data_[i] * rhs[i];
                return result;
            }
        }

        T norm() const {
            const size_type n = size();
            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);
                auto acc0 = PT::set1(T{0}), acc1 = acc0, acc2 = acc0, acc3 = acc0;
                for (size_type i = 0; i < unroll_end; i += STRIDE) {
                    auto v0 = PT::load(&data_[i]),          v1 = PT::load(&data_[i + PS]);
                    auto v2 = PT::load(&data_[i + PS * 2]), v3 = PT::load(&data_[i + PS * 3]);
                    acc0 = PT::fmadd(v0, v0, acc0);
                    acc1 = PT::fmadd(v1, v1, acc1);
                    acc2 = PT::fmadd(v2, v2, acc2);
                    acc3 = PT::fmadd(v3, v3, acc3);
                }
                acc0 = PT::add(PT::add(acc0, acc1), PT::add(acc2, acc3));
                T sum_sq = PT::reduce_add(acc0);
                for (size_type i = unroll_end; i < vec_end; i += PS) {
                    auto v = PT::load(&data_[i]);
                    sum_sq += PT::reduce_add(PT::mul(v, v));
                }
                for (size_type i = vec_end; i < n; ++i)
                    sum_sq += data_[i] * data_[i];
                return static_cast<T>(std::sqrt(static_cast<double>(sum_sq)));
            } else {
                T sum_squares = T{0};
                for (size_type i = 0; i < n; ++i)
                    sum_squares += data_[i] * data_[i];
                using std::sqrt;
                return sqrt(sum_squares);  // ADL: calx::sqrt(Float) も呼ばれる
            }
        }

        /// L1 ノルム (マンハッタン距離): Σ|x_i|
        T norm_l1() const {
            T sum = T{ 0 };
            for (size_type i = 0; i < size(); ++i) {
                if constexpr (std::is_arithmetic_v<T>)
                    sum += static_cast<T>(std::abs(data_[i]));
                else {
                    using std::abs; // ADL fallback
                    sum += abs(data_[i]);
                }
            }
            return sum;
        }

        /// L∞ ノルム (チェビシェフ距離): max|x_i|
        T norm_linf() const {
            T max_val = T{ 0 };
            for (size_type i = 0; i < size(); ++i) {
                T abs_val;
                if constexpr (std::is_arithmetic_v<T>)
                    abs_val = static_cast<T>(std::abs(data_[i]));
                else {
                    using std::abs;
                    abs_val = abs(data_[i]);
                }
                if (abs_val > max_val) max_val = abs_val;
            }
            return max_val;
        }

        /// 一般 Lp ノルム: (Σ|x_i|^p)^(1/p)
        T norm_lp(double p) const {
            if (p < 1.0) {
                throw MathError("norm_lp: p must be >= 1");
            }
            if (std::isinf(p)) {
                return norm_linf();
            }
            if (p == 1.0) {
                return norm_l1();
            }
            if (p == 2.0) {
                return norm();
            }
            double sum = 0.0;
            for (size_type i = 0; i < size(); ++i) {
                sum += std::pow(std::abs(static_cast<double>(data_[i])), p);
            }
            return static_cast<T>(std::pow(sum, 1.0 / p));
        }

        Vector normalized() const {
            Vector result = *this;
            T length = norm();
            if (length > T{ 0 }) {
                result /= length;
            }
            return result;
        }

        void normalize() {
            T length = norm();
            if (length > T{ 0 }) {
                *this /= length;
            }
        }

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

        /// 先頭 n 要素のビュー (コピーなし)
        detail::VectorView<Vector> head(size_type n) {
            assert(n <= size() && "Vector::head: n exceeds size");
            return detail::VectorView<Vector>(*this, 0, n);
        }
        detail::ConstVectorView<Vector> head(size_type n) const {
            assert(n <= size() && "Vector::head: n exceeds size");
            return detail::ConstVectorView<Vector>(*this, 0, n);
        }

        /// 末尾 n 要素のビュー (コピーなし)
        detail::VectorView<Vector> tail(size_type n) {
            assert(n <= size() && "Vector::tail: n exceeds size");
            return detail::VectorView<Vector>(*this, size() - n, n);
        }
        detail::ConstVectorView<Vector> tail(size_type n) const {
            assert(n <= size() && "Vector::tail: n exceeds size");
            return detail::ConstVectorView<Vector>(*this, size() - n, n);
        }

        /// 位置 offset から count 要素のビュー (コピーなし)
        detail::VectorView<Vector> segment(size_type offset, size_type count) {
            assert(offset + count <= size() && "Vector::segment: range exceeds size");
            return detail::VectorView<Vector>(*this, offset, count);
        }
        detail::ConstVectorView<Vector> segment(size_type offset, size_type count) const {
            assert(offset + count <= size() && "Vector::segment: range exceeds size");
            return detail::ConstVectorView<Vector>(*this, offset, count);
        }

    private:
        storage_type data_;
    };

    // 二項演算子 (Vector + Vector, scalar * Vector 等) は expr_templates.hpp の
    // VecExpr 演算子に委譲。式テンプレート経由で SIMD パケット評価が有効になる。
    // 代入時 (Vector v = expr) に Vector(const VecExpr<E>&) が暗黙変換として働く。

    // 内積
    template<typename T>
    T dot(const Vector<T>& lhs, const Vector<T>& rhs) {
        return lhs.dot(rhs);
    }

    // ノルム (L2)
    template<typename T>
    T norm(const Vector<T>& vec) {
        return vec.norm();
    }

    // norm2: norm の別名 (ODE ソルバー等の汎用アルゴリズム用)
    template<typename T>
    T norm2(const Vector<T>& vec) {
        return vec.norm();
    }

    // VecExpr 版 (ET 式ノードから直接計算 — SIMD パケット評価)
    template<typename E>
    auto norm(const VecExpr<E>& expr) {
        const auto& e = expr.derived();
        using T = typename E::value_type;
        const std::size_t n = e.size();
        if constexpr (has_simd_packet_v<T>) {
            using PT = PacketTraits<T>;
            constexpr std::size_t PS = PT::size;
            const std::size_t vec_end = n - (n % PS);
            auto acc = PT::set1(T{0});
            for (std::size_t i = 0; i < vec_end; i += PS) {
                auto v = e.packet(i);
                acc = PT::fmadd(v, v, acc);
            }
            T sum = PT::reduce_add(acc);
            for (std::size_t i = vec_end; i < n; ++i) {
                T v = e[i];
                sum += v * v;
            }
            return static_cast<T>(std::sqrt(static_cast<double>(sum)));
        } else {
            T sum = T{0};
            for (std::size_t i = 0; i < n; ++i) {
                T v = e[i];
                sum += v * v;
            }
            using std::sqrt;
            return sqrt(sum);  // ADL: calx::sqrt(Float) も呼ばれる
        }
    }

    template<typename E>
    auto norm2(const VecExpr<E>& expr) {
        return norm(expr);
    }

    // L1 ノルム
    template<typename T>
    T norm_l1(const Vector<T>& vec) {
        return vec.norm_l1();
    }

    // L∞ ノルム
    template<typename T>
    T norm_linf(const Vector<T>& vec) {
        return vec.norm_linf();
    }

    // 一般 Lp ノルム
    template<typename T>
    T norm_lp(const Vector<T>& vec, double p) {
        return vec.norm_lp(p);
    }

    // 正規化
    template<typename T>
    Vector<T> normalized(const Vector<T>& vec) {
        return vec.normalized();
    }

    // 固定長ベクトル
    template<typename T, std::size_t N>
    class StaticVector : public VecExpr<StaticVector<T, N>> {
    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;

        // コンストラクタ
        StaticVector() : data_{} {}

        explicit StaticVector(const T& value) {
            std::fill(data_.begin(), data_.end(), value);
        }

        StaticVector(std::initializer_list<T> init) {
            std::size_t i = 0;
            for (const auto& value : init) {
                if (i >= N) break;
                data_[i++] = value;
            }
        }

        // コピー/ムーブ
        StaticVector(const StaticVector&) = default;
        StaticVector(StaticVector&&) noexcept = default;
        StaticVector& operator=(const StaticVector&) = default;
        StaticVector& operator=(StaticVector&&) noexcept = default;

        // 式テンプレートからの構築
        template<typename E>
        StaticVector(const VecExpr<E>& expr) : data_{} {
            assignFromExpr(expr.derived());
        }

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

    private:
        template<typename E>
        void assignFromExpr(const E& e) {
            if constexpr (has_simd_packet_v<T> && N >= PacketTraits<T>::size) {
                constexpr size_type PS = PacketTraits<T>::size;
                constexpr size_type vec_end = N - (N % PS);
                for (size_type i = 0; i < vec_end; i += PS)
                    PacketTraits<T>::store(&data_[i], e.packet(i));
                for (size_type i = vec_end; i < N; ++i)
                    data_[i] = static_cast<T>(e[i]);
            } else {
                for (size_type i = 0; i < N; ++i)
                    data_[i] = static_cast<T>(e[i]);
            }
        }

    public:

        // 要素アクセス
        reference operator[](size_type index) {
            assert(index < N && "StaticVector::operator[]: index out of range");
            return data_[index];
        }

        const_reference operator[](size_type index) const {
            assert(index < N && "StaticVector::operator[]: index out of range");
            return data_[index];
        }

        reference at(size_type index) {
            if (index >= N) throw std::out_of_range("StaticVector index out of range");
            return data_[index];
        }

        const_reference at(size_type index) const {
            if (index >= N) throw std::out_of_range("StaticVector index out of range");
            return data_[index];
        }

        // イテレータ
        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() + N; }
        const_iterator end() const noexcept { return data_.data() + N; }
        const_iterator cend() const noexcept { return data_.data() + N; }

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

        // データアクセス
        T* data() noexcept { return data_.data(); }
        const T* data() const noexcept { return data_.data(); }

        // SIMD パケットロード
        auto packet(std::size_t i) const {
            return PacketTraits<T>::load(&data_[i]);
        }

        // ゼロベクトル
        void zero() {
            std::fill(data_.begin(), data_.end(), numeric_traits<T>::zero());
        }

        // 演算子オーバーロード
        StaticVector& operator+=(const StaticVector& rhs) {
            for (size_type i = 0; i < N; ++i) {
                data_[i] += rhs[i];
            }
            return *this;
        }

        StaticVector& operator-=(const StaticVector& rhs) {
            for (size_type i = 0; i < N; ++i) {
                data_[i] -= rhs[i];
            }
            return *this;
        }

        StaticVector& operator*=(const T& scalar) {
            for (size_type i = 0; i < N; ++i) {
                data_[i] *= scalar;
            }
            return *this;
        }

        StaticVector& operator/=(const T& scalar) {
            // ゼロ除算チェック
            if (scalar == T{ 0 }) {
                throw std::invalid_argument("Division by zero");
            }
            for (size_type i = 0; i < N; ++i) {
                data_[i] /= scalar;
            }
            return *this;
        }

        // 式テンプレートからの複合代入
        template<typename E>
        StaticVector& operator+=(const VecExpr<E>& rhs) {
            const E& r = rhs.derived();
            if constexpr (has_simd_packet_v<T> && N >= PacketTraits<T>::size) {
                constexpr size_type PS = PacketTraits<T>::size;
                constexpr size_type vec_end = N - (N % PS);
                for (size_type i = 0; i < vec_end; i += PS) {
                    auto cur = PacketTraits<T>::load(&data_[i]);
                    PacketTraits<T>::store(&data_[i], PacketTraits<T>::add(cur, r.packet(i)));
                }
                for (size_type i = vec_end; i < N; ++i)
                    data_[i] += static_cast<T>(r[i]);
            } else {
                for (size_type i = 0; i < N; ++i)
                    data_[i] += static_cast<T>(r[i]);
            }
            return *this;
        }

        template<typename E>
        StaticVector& operator-=(const VecExpr<E>& rhs) {
            const E& r = rhs.derived();
            if constexpr (has_simd_packet_v<T> && N >= PacketTraits<T>::size) {
                constexpr size_type PS = PacketTraits<T>::size;
                constexpr size_type vec_end = N - (N % PS);
                for (size_type i = 0; i < vec_end; i += PS) {
                    auto cur = PacketTraits<T>::load(&data_[i]);
                    PacketTraits<T>::store(&data_[i], PacketTraits<T>::sub(cur, r.packet(i)));
                }
                for (size_type i = vec_end; i < N; ++i)
                    data_[i] -= static_cast<T>(r[i]);
            } else {
                for (size_type i = 0; i < N; ++i)
                    data_[i] -= static_cast<T>(r[i]);
            }
            return *this;
        }

        // ベクトル演算
        T dot(const StaticVector& rhs) const {
            if constexpr (has_simd_packet_v<T> && N >= PacketTraits<T>::size) {
                using PT = PacketTraits<T>;
                constexpr size_type PS = PT::size;
                constexpr size_type vec_end = N - (N % PS);
                auto acc = PT::set1(T{0});
                for (size_type i = 0; i < vec_end; i += PS)
                    acc = PT::fmadd(PT::load(&data_[i]), PT::load(&rhs.data_[i]), acc);
                T result = PT::reduce_add(acc);
                for (size_type i = vec_end; i < N; ++i)
                    result += data_[i] * rhs.data_[i];
                return result;
            } else {
                T result = T{0};
                for (size_type i = 0; i < N; ++i)
                    result += data_[i] * rhs[i];
                return result;
            }
        }

        T norm() const {
            if constexpr (has_simd_packet_v<T> && N >= PacketTraits<T>::size) {
                using PT = PacketTraits<T>;
                constexpr size_type PS = PT::size;
                constexpr size_type vec_end = N - (N % PS);
                auto acc = PT::set1(T{0});
                for (size_type i = 0; i < vec_end; i += PS) {
                    auto v = PT::load(&data_[i]);
                    acc = PT::fmadd(v, v, acc);
                }
                T sum_sq = PT::reduce_add(acc);
                for (size_type i = vec_end; i < N; ++i)
                    sum_sq += data_[i] * data_[i];
                return static_cast<T>(std::sqrt(static_cast<double>(sum_sq)));
            } else {
                T sum_squares = T{0};
                for (size_type i = 0; i < N; ++i)
                    sum_squares += data_[i] * data_[i];
                using std::sqrt;
                return sqrt(sum_squares);  // ADL: calx::sqrt(Float) も呼ばれる
            }
        }

        /// L1 ノルム: Σ|x_i|
        T norm_l1() const {
            T sum = T{ 0 };
            for (size_type i = 0; i < N; ++i) {
                if constexpr (std::is_arithmetic_v<T>)
                    sum += static_cast<T>(std::abs(data_[i]));
                else {
                    using std::abs;
                    sum += abs(data_[i]);
                }
            }
            return sum;
        }

        /// L∞ ノルム: max|x_i|
        T norm_linf() const {
            T max_val = T{ 0 };
            for (size_type i = 0; i < N; ++i) {
                T abs_val;
                if constexpr (std::is_arithmetic_v<T>)
                    abs_val = static_cast<T>(std::abs(data_[i]));
                else {
                    using std::abs;
                    abs_val = abs(data_[i]);
                }
                if (abs_val > max_val) max_val = abs_val;
            }
            return max_val;
        }

        /// 一般 Lp ノルム: (Σ|x_i|^p)^(1/p)
        T norm_lp(double p) const {
            if (p < 1.0) {
                throw MathError("norm_lp: p must be >= 1");
            }
            if (std::isinf(p)) {
                return norm_linf();
            }
            if (p == 1.0) {
                return norm_l1();
            }
            if (p == 2.0) {
                return norm();
            }
            double sum = 0.0;
            for (size_type i = 0; i < N; ++i) {
                sum += std::pow(std::abs(static_cast<double>(data_[i])), p);
            }
            return static_cast<T>(std::pow(sum, 1.0 / p));
        }

        StaticVector normalized() const {
            StaticVector result = *this;
            T length = norm();
            if (length > T{ 0 }) {
                result /= length;
            }
            return result;
        }

        void normalize() {
            T length = norm();
            if (length > T{ 0 }) {
                *this /= length;
            }
        }

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

        /// 先頭 n 要素のビュー
        detail::VectorView<StaticVector> head(size_type n) {
            assert(n <= N && "StaticVector::head: n exceeds size");
            return detail::VectorView<StaticVector>(*this, 0, n);
        }
        detail::ConstVectorView<StaticVector> head(size_type n) const {
            assert(n <= N && "StaticVector::head: n exceeds size");
            return detail::ConstVectorView<StaticVector>(*this, 0, n);
        }

        /// 末尾 n 要素のビュー
        detail::VectorView<StaticVector> tail(size_type n) {
            assert(n <= N && "StaticVector::tail: n exceeds size");
            return detail::VectorView<StaticVector>(*this, size() - n, n);
        }
        detail::ConstVectorView<StaticVector> tail(size_type n) const {
            assert(n <= N && "StaticVector::tail: n exceeds size");
            return detail::ConstVectorView<StaticVector>(*this, size() - n, n);
        }

        /// 位置 offset から count 要素のビュー
        detail::VectorView<StaticVector> segment(size_type offset, size_type count) {
            assert(offset + count <= N && "StaticVector::segment: range exceeds size");
            return detail::VectorView<StaticVector>(*this, offset, count);
        }
        detail::ConstVectorView<StaticVector> segment(size_type offset, size_type count) const {
            assert(offset + count <= N && "StaticVector::segment: range exceeds size");
            return detail::ConstVectorView<StaticVector>(*this, offset, count);
        }

    private:
        std::array<T, N> data_;
    };

    // StaticVector の二項演算子も expr_templates.hpp の VecExpr 演算子に委譲。

    // 内積（StaticVector用）
    template<typename T, std::size_t N>
    T dot(const StaticVector<T, N>& lhs, const StaticVector<T, N>& rhs) {
        return lhs.dot(rhs);
    }

    // ノルム（StaticVector用）
    template<typename T, std::size_t N>
    T norm(const StaticVector<T, N>& vec) {
        return vec.norm();
    }

    // L1 ノルム（StaticVector用）
    template<typename T, std::size_t N>
    T norm_l1(const StaticVector<T, N>& vec) {
        return vec.norm_l1();
    }

    // L∞ ノルム（StaticVector用）
    template<typename T, std::size_t N>
    T norm_linf(const StaticVector<T, N>& vec) {
        return vec.norm_linf();
    }

    // 一般 Lp ノルム（StaticVector用）
    template<typename T, std::size_t N>
    T norm_lp(const StaticVector<T, N>& vec, double p) {
        return vec.norm_lp(p);
    }

    // 正規化（StaticVector用）
    template<typename T, std::size_t N>
    StaticVector<T, N> normalized(const StaticVector<T, N>& vec) {
        return vec.normalized();
    }

    // 3D 外積 (cross product)
    template<typename T>
    StaticVector<T, 3> cross(const StaticVector<T, 3>& a, const StaticVector<T, 3>& b) {
        return StaticVector<T, 3>{
            a[1] * b[2] - a[2] * b[1],
            a[2] * b[0] - a[0] * b[2],
            a[0] * b[1] - a[1] * b[0]
        };
    }

    // =========================================================================
    // BLAS Level 1 融合関数
    // 一時オブジェクトを生成せず単一ループで処理。コンパイラが SIMD 自動ベクタライズ可能。
    // =========================================================================

    // y += alpha * x  (DAXPY)
    template<typename T>
    void axpy(T alpha, const Vector<T>& x, Vector<T>& y) {
        const auto n = x.size();
        if (n != y.size())
            throw DimensionError("axpy: size mismatch");
        const T* __restrict xp = x.data();
        T* __restrict yp = y.data();
        for (std::size_t i = 0; i < n; ++i)
            yp[i] += alpha * xp[i];
    }

    // result = alpha * x + beta * y  (out-of-place axpby)
    template<typename T>
    void axpby(T alpha, const Vector<T>& x, T beta, const Vector<T>& y, Vector<T>& result) {
        const auto n = x.size();
        if (n != y.size() || n != result.size())
            throw DimensionError("axpby: size mismatch");
        const T* __restrict xp = x.data();
        const T* __restrict yp = y.data();
        T* __restrict rp = result.data();
        for (std::size_t i = 0; i < n; ++i)
            rp[i] = alpha * xp[i] + beta * yp[i];
    }

    // y = alpha * x + beta * y  (in-place axpby)
    template<typename T>
    void axpby(T alpha, const Vector<T>& x, T beta, Vector<T>& y) {
        const auto n = x.size();
        if (n != y.size())
            throw DimensionError("axpby: size mismatch");
        const T* __restrict xp = x.data();
        T* __restrict yp = y.data();
        for (std::size_t i = 0; i < n; ++i)
            yp[i] = alpha * xp[i] + beta * yp[i];
    }

    // StaticVector 版
    template<typename T, std::size_t N>
    void axpy(T alpha, const StaticVector<T, N>& x, StaticVector<T, N>& y) {
        for (std::size_t i = 0; i < N; ++i)
            y[i] += alpha * x[i];
    }

    template<typename T, std::size_t N>
    void axpby(T alpha, const StaticVector<T, N>& x, T beta, const StaticVector<T, N>& y,
               StaticVector<T, N>& result) {
        for (std::size_t i = 0; i < N; ++i)
            result[i] = alpha * x[i] + beta * y[i];
    }

    template<typename T, std::size_t N>
    void axpby(T alpha, const StaticVector<T, N>& x, T beta, StaticVector<T, N>& y) {
        for (std::size_t i = 0; i < N; ++i)
            y[i] = alpha * x[i] + beta * y[i];
    }

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

    /**
     * @brief 外部メモリ上のベクトルビュー
     *
     * コピーなしで外部配列を Vector 互換オブジェクトとして参照する。
     * @code
     * double data[] = {1.0, 2.0, 3.0};
     * auto v = calx::VectorMap<double>(data, 3);
     * v[1] = 5.0;  // data[1] が変更される
     * @endcode
     */
    template<typename T>
    class VectorMap {
    public:
        using value_type = T;
        using size_type = std::size_t;

        VectorMap(T* data, size_type size) : data_(data), size_(size) {}

        size_type size() const { return size_; }
        T* data() { return data_; }
        const T* data() const { return data_; }

        T& operator[](size_type i) { assert(i < size_ && "VectorMap::operator[]: index out of range"); return data_[i]; }
        const T& operator[](size_type i) const { assert(i < size_ && "VectorMap::operator[]: index out of range"); return data_[i]; }

        /// Vector へのコピー
        operator Vector<T>() const {
            Vector<T> v(size_);
            for (size_type i = 0; i < size_; ++i) v[i] = data_[i];
            return v;
        }

        /// Vector からの代入
        VectorMap& operator=(const Vector<T>& v) {
            for (size_type i = 0; i < std::min(size_, v.size()); ++i)
                data_[i] = v[i];
            return *this;
        }

    private:
        T* data_;
        size_type size_;
    };

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

        ConstVectorMap(const T* data, size_type size) : data_(data), size_(size) {}

        size_type size() const { return size_; }
        const T* data() const { return data_; }
        const T& operator[](size_type i) const { assert(i < size_ && "ConstVectorMap::operator[]: index out of range"); return data_[i]; }

        operator Vector<T>() const {
            Vector<T> v(size_);
            for (size_type i = 0; i < size_; ++i) v[i] = data_[i];
            return v;
        }

    private:
        const T* data_;
        size_type size_;
    };

} // namespace calx

#endif // CALX_VECTOR_HPP