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

// vector.hpp
// Vector class definitions
//
// This file defines the vector classes of the MKL algebra library.
// It includes implementations of both fixed-size and dynamic vectors.

#ifndef SANGI_VECTOR_HPP
#define SANGI_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 <math/core/BaseVector.hpp>

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

namespace sangi {

    /// Allocator that skips value-initialization on resize (leaves POD types uninitialized)
    template<typename T>
    struct DefaultInitAllocator : std::allocator<T> {
        using std::allocator<T>::allocator;
        template<typename U> struct rebind { using other = DefaultInitAllocator<U>; };

        // default-init: do not touch memory for POD types
        void construct(T* p) noexcept(std::is_nothrow_default_constructible_v<T>) {
            ::new (static_cast<void*>(p)) T;
        }
        // Normal construction when arguments are supplied
        template<typename... Args>
        void construct(T* p, Args&&... args) {
            ::new (static_cast<void*>(p)) T(std::forward<Args>(args)...);
        }
    };

    /// Tag type for the uninitialized constructor
    struct uninitialized_t {};
    /// Tag value for the uninitialized constructor
    inline constexpr uninitialized_t uninitialized{};

    // Definition of the Vector class forward-declared in basic_types.hpp
    // Shares sangi::BaseVector<T> as common parent and is polymorphic with StaticVector.
    template<typename T>
    class Vector : public BaseVector<T>, public VecExpr<Vector<T>> {
        using storage_type = std::vector<T, DefaultInitAllocator<T>>;
    public:
        using BVBase = BaseVector<T>;
        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;

        // Constructors
        Vector() : BVBase(), data_() {}

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

        /// Uninitialized constructor: only allocates memory; values are not set
        Vector(size_type size, uninitialized_t) : BVBase(), data_(size) {
            sync_basevector();
        }

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

        Vector(std::initializer_list<T> init) : BVBase(), data_(init) {
            sync_basevector();
        }

        // Copy/move constructors and assignments
        Vector(const Vector& other) : BVBase(), data_(other.data_), VecExpr<Vector<T>>() {
            sync_basevector();
        }
        Vector(Vector&& other) noexcept : BVBase(), data_(std::move(other.data_)),
                                          VecExpr<Vector<T>>() {
            sync_basevector();
        }
        Vector& operator=(const Vector& other) {
            if (this != &other) {
                data_ = other.data_;
                sync_basevector();
            }
            return *this;
        }
        Vector& operator=(Vector&& other) noexcept {
            if (this != &other) {
                data_ = std::move(other.data_);
                sync_basevector();
            }
            return *this;
        }

        // Construction from an expression template
        template<typename E>
        Vector(const VecExpr<E>& expr) : BVBase(), data_(expr.derived().size()) {
            assignFromExpr(expr.derived());
            sync_basevector();
        }

        // Deep-copy ctor from a BaseVector<T> view.
        // Required so that the `Vector<T> w = v;` pattern (where v is a BaseVector<T>
        // argument) works inside LinAlg function bodies. Whether the caller passes a
        // Vector<T> or a StaticVector, exactly one deep copy is performed here
        // (only when strictly necessary).
        // Note: although Vector<T> itself is derived from BaseVector<T>, Vector(const Vector&)
        // is preferred via exact-match overloading (this ctor handles StaticVector and
        // BaseVector view arguments).
        Vector(const BaseVector<T>& other)
            : BVBase(), data_(other.size()) {
            for (size_type i = 0; i < other.size(); ++i) {
                data_[i] = other[i];
            }
            sync_basevector();
        }

        // Assignment from an expression template
        template<typename E>
        Vector& operator=(const VecExpr<E>& expr) {
            const E& e = expr.derived();
            data_.resize(e.size());
            assignFromExpr(e);
            sync_basevector();
            return *this;
        }

    private:
        // Sync the BaseVector view with std::vector::data()
        void sync_basevector() noexcept {
            BVBase::set_view_nonowning(data_.data(), data_.size(), 1);
        }
    public:

    private:
        // Internal implementation of expression assignment via SIMD-packet evaluation (4x unrolled)
        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:

        // Destructor
        ~Vector() = default;

        // Element access
        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();
        }

        // Iterators
        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();
        }

        // Capacity
        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);
            sync_basevector();
        }

        void shrink_to_fit() {
            data_.shrink_to_fit();
            sync_basevector();
        }

        // Modification operations
        void clear() noexcept {
            data_.clear();
            sync_basevector();
        }

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

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

        // Wrapper to call the assign method of data_
        void assign(size_type count, const T& value) {
            data_.assign(count, value);
            sync_basevector();
        }

        // Data access
        T* data() noexcept {
            return data_.data();
        }

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

        // SIMD packet load (for packet evaluation in expression templates)
        auto packet(std::size_t i) const {
            return PacketTraits<T>::load(&data_[i]);
        }

        // Zero vector
        void zero() {
            std::fill(data_.begin(), data_.end(), numeric_traits<T>::zero());
        }

        // Operator overloads
        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) {
            // Division-by-zero check
            if (scalar == T{ 0 }) {
                throw std::invalid_argument("Division by zero");
            }

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

        // Compound assignment from an expression template
        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;
        }

        // Vector operations
        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 accumulators
                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: sangi::sqrt(Float) is also picked up
            }
        }

        /// L1 norm (Manhattan distance): sum |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-infinity norm (Chebyshev distance): 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;
        }

        /// General Lp norm: (sum |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));
        }

        // Coefficient-wise reductions (Eigen parity)

        /// Sum of all coefficients
        T sum() 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) {
                    acc0 = PT::add(acc0, PT::load(&data_[i]));
                    acc1 = PT::add(acc1, PT::load(&data_[i + PS]));
                    acc2 = PT::add(acc2, PT::load(&data_[i + PS * 2]));
                    acc3 = PT::add(acc3, PT::load(&data_[i + PS * 3]));
                }
                acc0 = PT::add(PT::add(acc0, acc1), PT::add(acc2, acc3));
                T s = PT::reduce_add(acc0);
                for (size_type i = unroll_end; i < vec_end; i += PS)
                    s += PT::reduce_add(PT::load(&data_[i]));
                for (size_type i = vec_end; i < n; ++i)
                    s += data_[i];
                return s;
            } else {
                T s = T{0};
                for (size_type i = 0; i < n; ++i) s += data_[i];
                return s;
            }
        }

        /// Arithmetic mean of all coefficients
        T mean() const {
            if (size() == 0) throw MathError("Vector::mean: empty vector");
            return sum() / static_cast<T>(size());
        }

        /// Index of the minimum coefficient (real types only)
        size_type argmin() const {
            static_assert(!numeric_traits<T>::is_complex,
                "argmin requires an ordered (non-complex) scalar type");
            if (size() == 0) throw MathError("Vector::argmin: empty vector");
            size_type idx = 0;
            for (size_type i = 1; i < size(); ++i)
                if (data_[i] < data_[idx]) idx = i;
            return idx;
        }

        /// Index of the maximum coefficient (real types only)
        size_type argmax() const {
            static_assert(!numeric_traits<T>::is_complex,
                "argmax requires an ordered (non-complex) scalar type");
            if (size() == 0) throw MathError("Vector::argmax: empty vector");
            size_type idx = 0;
            for (size_type i = 1; i < size(); ++i)
                if (data_[i] > data_[idx]) idx = i;
            return idx;
        }

        /// Minimum coefficient (real types only)
        T minCoeff() const { return data_[argmin()]; }
        /// Minimum coefficient, writing its index to `index`
        T minCoeff(size_type& index) const { index = argmin(); return data_[index]; }
        /// Maximum coefficient (real types only)
        T maxCoeff() const { return data_[argmax()]; }
        /// Maximum coefficient, writing its index to `index`
        T maxCoeff(size_type& index) const { index = argmax(); return data_[index]; }

        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;
            }
        }

        // Block operations API (Eigen style)

        /// View of the first n elements (no copy)
        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);
        }

        /// View of the last n elements (no copy)
        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);
        }

        /// View of `count` elements starting at position `offset` (no copy)
        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_;
    };

    // Binary operators (Vector + Vector, scalar * Vector, etc.) delegate to the
    // VecExpr operators in expr_templates.hpp. SIMD packet evaluation is enabled
    // via expression templates. During assignment (Vector v = expr),
    // Vector(const VecExpr<E>&) participates as an implicit conversion.

    // ============================================================================
    // dot / norm / norm2 — canonical free functions
    //
    // Concrete vectors and views (Vector, StaticVector, row/column views, or an
    // argument typed as BaseVector<T> itself) all bind to the single BaseVector<T>
    // overload below. Lazy expression nodes (a + b, alpha * v, ...) are NOT
    // BaseVectors (no storage) and are handled by the VecExpr overloads further
    // down, which evaluate SIMD packets without materializing.
    //
    // Convention matches Vector::dot / blas::dot: a symmetric (non-conjugated)
    // bilinear product. The contiguous path uses the same 4-accumulator unrolled
    // SIMD kernel as Vector::dot; strided views fall back to the stride-aware
    // operator[] loop.
    // ============================================================================
    template<typename T>
    T dot(const BaseVector<T>& lhs, const BaseVector<T>& rhs) {
        if (lhs.size() != rhs.size())
            throw DimensionError("dot: size mismatch");
        const std::size_t n = lhs.size();
        if constexpr (has_simd_packet_v<T>) {
            if (lhs.is_contiguous() && rhs.is_contiguous()) {
                using PT = PacketTraits<T>;
                constexpr std::size_t PS = PT::size;
                constexpr std::size_t STRIDE = PS * 4;  // 4 accumulators
                const std::size_t unroll_end = n - (n % STRIDE);
                const std::size_t vec_end = n - (n % PS);
                const T* a = lhs.data();
                const T* b = rhs.data();
                auto acc0 = PT::set1(T{0}), acc1 = acc0, acc2 = acc0, acc3 = acc0;
                for (std::size_t i = 0; i < unroll_end; i += STRIDE) {
                    acc0 = PT::fmadd(PT::load(&a[i]),          PT::load(&b[i]),          acc0);
                    acc1 = PT::fmadd(PT::load(&a[i + PS]),     PT::load(&b[i + PS]),     acc1);
                    acc2 = PT::fmadd(PT::load(&a[i + PS * 2]), PT::load(&b[i + PS * 2]), acc2);
                    acc3 = PT::fmadd(PT::load(&a[i + PS * 3]), PT::load(&b[i + PS * 3]), acc3);
                }
                acc0 = PT::add(PT::add(acc0, acc1), PT::add(acc2, acc3));
                T result = PT::reduce_add(acc0);
                for (std::size_t i = unroll_end; i < vec_end; i += PS)
                    result += PT::reduce_add(PT::mul(PT::load(&a[i]), PT::load(&b[i])));
                for (std::size_t i = vec_end; i < n; ++i)
                    result += a[i] * b[i];
                return result;
            }
        }
        T result = numeric_traits<T>::zero();
        for (std::size_t i = 0; i < n; ++i)
            result += lhs[i] * rhs[i];
        return result;
    }

    template<typename T>
    T norm(const BaseVector<T>& vec) {
        const T sq = dot(vec, vec);  // SIMD-accelerated sum of squares (contiguous path)
        if constexpr (numeric_traits<T>::is_complex) {
            using std::sqrt;
            return sqrt(sq);  // ADL: sangi::sqrt is also picked up
        } else {
            return static_cast<T>(std::sqrt(static_cast<double>(sq)));
        }
    }

    template<typename T>
    T norm2(const BaseVector<T>& vec) {
        return norm(vec);
    }

    // VecExpr version: pure lazy expression nodes only (a + b, alpha * v, ...),
    // computed directly via SIMD packet evaluation without materializing.
    // Concrete vectors/views derive from BaseVector too, so they are excluded here
    // and bind unambiguously to the BaseVector<T> overload above.
    template<typename E>
        requires (!std::is_base_of_v<BaseVector<typename E::value_type>, 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: sangi::sqrt(Float) is also picked up
        }
    }

    template<typename E>
        requires (!std::is_base_of_v<BaseVector<typename E::value_type>, E>)
    auto norm2(const VecExpr<E>& expr) {
        return norm(expr);
    }

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

    // L-infinity norm
    template<typename T>
    T norm_linf(const Vector<T>& vec) {
        return vec.norm_linf();
    }

    // General Lp norm
    template<typename T>
    T norm_lp(const Vector<T>& vec, double p) {
        return vec.norm_lp(p);
    }

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

    // Coefficient-wise reductions (free-function form, Eigen parity)
    template<typename T>
    T sum(const Vector<T>& vec) { return vec.sum(); }
    template<typename T>
    T mean(const Vector<T>& vec) { return vec.mean(); }
    template<typename T>
    T minCoeff(const Vector<T>& vec) { return vec.minCoeff(); }
    template<typename T>
    T maxCoeff(const Vector<T>& vec) { return vec.maxCoeff(); }
    template<typename T>
    std::size_t argmin(const Vector<T>& vec) { return vec.argmin(); }
    template<typename T>
    std::size_t argmax(const Vector<T>& vec) { return vec.argmax(); }

    // Fixed-size vector
    // Shares sangi::BaseVector<T> as common parent and is polymorphic with Vector.
    template<typename T, std::size_t N>
    class StaticVector : public BaseVector<T>, public VecExpr<StaticVector<T, N>> {
    public:
        // Compile-time size trait. Shadows BaseVector::static_size (-1) so that
        // LinAlg algorithms can static_assert dimension agreement at build time.
        static constexpr std::ptrdiff_t static_size = static_cast<std::ptrdiff_t>(N);

        using BVBase = BaseVector<T>;
        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;

        // Constructors
        StaticVector() : BVBase(), data_{} {
            sync_basevector();
        }

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

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

        // Copy / move
        StaticVector(const StaticVector& other) : BVBase(), data_(other.data_),
                                                   VecExpr<StaticVector<T, N>>() {
            sync_basevector();
        }
        StaticVector(StaticVector&& other) noexcept : BVBase(), data_(std::move(other.data_)),
                                                      VecExpr<StaticVector<T, N>>() {
            sync_basevector();
        }
        StaticVector& operator=(const StaticVector& other) {
            if (this != &other) {
                data_ = other.data_;
                // No sync needed -- std::array updates in place and data() address is stable
            }
            return *this;
        }
        StaticVector& operator=(StaticVector&&) noexcept = default;

        // Construction from an expression template
        template<typename E>
        StaticVector(const VecExpr<E>& expr) : BVBase(), data_{} {
            assignFromExpr(expr.derived());
            sync_basevector();
        }

        // Assignment from an expression template
        template<typename E>
        StaticVector& operator=(const VecExpr<E>& expr) {
            assignFromExpr(expr.derived());
            return *this;
        }

    private:
        // Sync the BaseVector view with std::array::data() (non-owning, since storage is fixed)
        void sync_basevector() noexcept {
            BVBase::set_view_nonowning(data_.data(), N, 1);
        }
    public:

    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:

        // Element access
        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];
        }

        // Iterators
        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; }

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

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

        // SIMD packet load
        auto packet(std::size_t i) const {
            return PacketTraits<T>::load(&data_[i]);
        }

        // Zero vector
        void zero() {
            std::fill(data_.begin(), data_.end(), numeric_traits<T>::zero());
        }

        // Operator overloads
        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) {
            // Division-by-zero check
            if (scalar == T{ 0 }) {
                throw std::invalid_argument("Division by zero");
            }
            for (size_type i = 0; i < N; ++i) {
                data_[i] /= scalar;
            }
            return *this;
        }

        // Compound assignment from an expression template
        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;
        }

        // Vector operations
        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: sangi::sqrt(Float) is also picked up
            }
        }

        /// L1 norm: sum |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-infinity norm: 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;
        }

        /// General Lp norm: (sum |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));
        }

        // Coefficient-wise reductions (Eigen parity)

        /// Sum of all coefficients
        T sum() const {
            T s = T{0};
            for (size_type i = 0; i < N; ++i) s += data_[i];
            return s;
        }

        /// Arithmetic mean of all coefficients
        T mean() const {
            static_assert(N > 0, "StaticVector::mean: empty vector");
            return sum() / static_cast<T>(N);
        }

        /// Index of the minimum coefficient (real types only)
        size_type argmin() const {
            static_assert(!numeric_traits<T>::is_complex,
                "argmin requires an ordered (non-complex) scalar type");
            static_assert(N > 0, "StaticVector::argmin: empty vector");
            size_type idx = 0;
            for (size_type i = 1; i < N; ++i)
                if (data_[i] < data_[idx]) idx = i;
            return idx;
        }

        /// Index of the maximum coefficient (real types only)
        size_type argmax() const {
            static_assert(!numeric_traits<T>::is_complex,
                "argmax requires an ordered (non-complex) scalar type");
            static_assert(N > 0, "StaticVector::argmax: empty vector");
            size_type idx = 0;
            for (size_type i = 1; i < N; ++i)
                if (data_[i] > data_[idx]) idx = i;
            return idx;
        }

        T minCoeff() const { return data_[argmin()]; }
        T minCoeff(size_type& index) const { index = argmin(); return data_[index]; }
        T maxCoeff() const { return data_[argmax()]; }
        T maxCoeff(size_type& index) const { index = argmax(); return data_[index]; }

        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;
            }
        }

        // Block operations API (Eigen style)

        /// View of the first n elements
        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);
        }

        /// View of the last n elements
        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);
        }

        /// View of `count` elements starting at position `offset`
        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's binary operators also delegate to the VecExpr operators in expr_templates.hpp.

    // Note: dot / norm / norm2 for StaticVector are provided by the single
    // BaseVector<T> overloads above (StaticVector derives from BaseVector<T>).

    // L1 norm (for StaticVector)
    template<typename T, std::size_t N>
    T norm_l1(const StaticVector<T, N>& vec) {
        return vec.norm_l1();
    }

    // L-infinity norm (for StaticVector)
    template<typename T, std::size_t N>
    T norm_linf(const StaticVector<T, N>& vec) {
        return vec.norm_linf();
    }

    // General Lp norm (for StaticVector)
    template<typename T, std::size_t N>
    T norm_lp(const StaticVector<T, N>& vec, double p) {
        return vec.norm_lp(p);
    }

    // Normalize (for 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]
        };
    }

    // =========================================================================
    // Fused BLAS Level 1 functions
    // Processed in a single loop without producing temporaries. The compiler
    // can auto-vectorize using 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 versions
    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 -- zero-copy view onto external memory (analogous to Eigen::Map)
    // =========================================================================

    /**
     * @brief Vector view onto external memory
     *
     * References an external array as a Vector-compatible object without copying.
     * @code
     * double data[] = {1.0, 2.0, 3.0};
     * auto v = sangi::VectorMap<double>(data, 3);
     * v[1] = 5.0;  // data[1] is modified
     * @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]; }

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

        /// Assignment from a 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_;
    };

    // Vector standard-output operator
    template<typename T>
    std::ostream& operator<<(std::ostream& os, const Vector<T>& v) {
        os << "[";
        for (std::size_t i = 0; i < v.size(); ++i) {
            if (i > 0) os << ", ";
            os << v[i];
        }
        os << "]";
        return os;
    }

    // Vector LaTeX output (column vector)
    // style: "pmatrix" (round brackets), "bmatrix" (square brackets), "vmatrix" (vertical bars)
    template<typename T>
    [[nodiscard]] std::string toLatex(const Vector<T>& v,
        const std::string& style = "pmatrix", std::size_t precision = 6) {
        std::ostringstream oss;
        std::ios::fmtflags old_flags = oss.flags();
        std::streamsize old_prec = oss.precision();
        oss << std::setprecision(precision);
        oss << "\\begin{" << style << "}\n";
        for (std::size_t i = 0; i < v.size(); ++i) {
            oss << v[i];
            if (i < v.size() - 1) oss << " \\\\";
            oss << "\n";
        }
        oss << "\\end{" << style << "}";
        oss.flags(old_flags);
        oss.precision(old_prec);
        return oss.str();
    }

} // namespace sangi

// ---------------------------------------------------------------------------
// extern template -- suppresses re-instantiation for float/double
// The definitions live in src/math/core/vector.cpp (sangi_matrix.lib).
// Other element types T are instantiated normally from the header.
// ---------------------------------------------------------------------------
#ifndef SANGI_VECTOR_EXPLICIT_INSTANTIATION
namespace sangi {
    extern template class Vector<float>;
    extern template class Vector<double>;
    extern template std::ostream& operator<<(std::ostream&, const Vector<float>&);
    extern template std::ostream& operator<<(std::ostream&, const Vector<double>&);
    extern template std::string toLatex(const Vector<float>&, const std::string&, std::size_t);
    extern template std::string toLatex(const Vector<double>&, const std::string&, std::size_t);
} // namespace sangi
#endif

#endif // SANGI_VECTOR_HPP