diff --git a/source/hpr/container.cpp b/source/hpr/container.cpp
deleted file mode 100644
index aa966c8..0000000
--- a/source/hpr/container.cpp
+++ /dev/null
@@ -1,2 +0,0 @@
-
-#include <hpr/container.hpp>
diff --git a/source/hpr/container.hpp b/source/hpr/container.hpp
index 106210c..52ae007 100644
--- a/source/hpr/container.hpp
+++ b/source/hpr/container.hpp
@@ -1,4 +1,7 @@
 #pragma once
 
-#include <iterator>
-
+#include <hpr/container/iterator.hpp>
+#include <hpr/container/sequence.hpp>
+#include <hpr/container/static_array.hpp>
+#include <hpr/container/dynamic_array.hpp>
+#include <hpr/container/tree_node.hpp>
diff --git a/source/hpr/container/dynamic_array.hpp b/source/hpr/container/dynamic_array.hpp
new file mode 100644
index 0000000..8b3abe8
--- /dev/null
+++ b/source/hpr/container/dynamic_array.hpp
@@ -0,0 +1,307 @@
+#pragma once
+
+#include <hpr/container/sequence.hpp>
+
+#include <functional>
+#include <limits>
+
+
+namespace hpr
+{
+    // forward declaration
+
+    template <typename T>
+    class DynamicArray;
+
+    // type traits
+
+    template <typename T>
+    struct is_sequence<DynamicArray<T>> : public std::true_type
+    {};
+
+    // aliases
+
+    template <typename T>
+    using darray = DynamicArray<T>;
+
+    // class definition
+
+    template <typename T>
+    class DynamicArray : public Sequence<T>
+    {
+
+        using base = Sequence<T>;
+
+    public:
+
+        using difference_type = std::ptrdiff_t;
+        using value_type = T;
+        using size_type = Size;
+        using pointer = T*;
+        using reference = T&;
+        using iterator = Iterator<T>;
+        using const_pointer = T const*;
+        using const_reference = T const&;
+        using const_iterator = Iterator<const T>;
+
+    public:
+
+        friend constexpr
+        void swap(DynamicArray& main, DynamicArray& other) noexcept
+        {
+            using std::swap;
+            swap(static_cast<base&>(main), static_cast<base&>(other));
+        }
+
+        constexpr
+        DynamicArray() noexcept :
+            base {}
+        {}
+
+        constexpr explicit
+        DynamicArray(const base& b) noexcept :
+            base {b}
+        {}
+
+        constexpr explicit
+        DynamicArray(base&& b) noexcept :
+            base {std::forward<base>(b)}
+        {}
+
+        constexpr
+        DynamicArray(const DynamicArray& arr) noexcept :
+            base {static_cast<base>(arr)}
+        {}
+
+        //! Move constructor
+        constexpr
+        DynamicArray(DynamicArray&& arr) noexcept :
+            base {std::forward<base>(static_cast<base>(arr))}
+        {}
+
+        template <typename Iter>
+        constexpr
+        DynamicArray(Iter start, Iter end) :
+            base {start, end}
+        {}
+
+        constexpr
+        DynamicArray(typename base::iterator start, typename base::iterator end) :
+            base {start, end}
+        {}
+
+        constexpr
+        DynamicArray(typename base::const_iterator start, typename base::const_iterator end) :
+            base {start, end}
+        {}
+
+        constexpr
+        DynamicArray(typename base::iterator start, typename base::iterator end, size_type capacity) :
+            base {start, end, capacity}
+        {}
+
+        constexpr
+        DynamicArray(typename base::const_iterator start, typename base::const_iterator end, size_type capacity) :
+            base {start, end, capacity}
+        {}
+
+        constexpr
+        DynamicArray(std::initializer_list<value_type> list) :
+            base {list.begin(), list.end()}
+        {}
+
+        constexpr
+        DynamicArray(size_type size, value_type value) noexcept :
+            base {size, value}
+        {}
+
+        constexpr
+        DynamicArray& operator=(base other) noexcept
+        {
+            swap(*this, other);
+            return *this;
+        }
+
+        constexpr
+        DynamicArray& operator=(DynamicArray other) noexcept
+        {
+            swap(*this, other);
+            return *this;
+        }
+
+        constexpr
+        DynamicArray& operator=(DynamicArray&& arr) noexcept
+        {
+            swap(*this, arr);
+            return *this;
+        }
+
+        virtual
+        ~DynamicArray() = default;
+
+        // Member functions
+
+        [[nodiscard]] virtual constexpr
+        size_type capacity() const noexcept
+        {
+            return base::p_capacity;
+        }
+
+        [[nodiscard]] virtual constexpr
+        bool is_empty() const
+        {
+            return base::begin() == base::end();
+        }
+
+        virtual constexpr
+        iterator storage_end()
+        {
+            return iterator(base::p_start + base::p_capacity);
+        }
+
+        virtual constexpr
+        const_iterator storage_end() const
+        {
+            return const_iterator(base::p_start + base::p_capacity);
+        }
+
+        virtual constexpr
+        void resize(size_type newCapacity)
+        {
+            if (newCapacity == base::p_capacity)
+                return;
+
+            if (std::numeric_limits<size_type>::max() - base::size() < newCapacity)
+                throw hpr::LengthError("Invalid capacity value passed");
+
+            if (newCapacity > base::p_capacity)
+            {
+                DynamicArray tmp {base::begin(), base::end(), newCapacity};
+                swap(*this, tmp);
+            }
+            else
+            {
+                if (newCapacity >= base::p_size)
+                {
+                    DynamicArray tmp {base::begin(), base::end()};
+                    swap(*this, tmp);
+                }
+                else
+                {
+                    DynamicArray tmp {base::begin(), base::begin() + newCapacity};
+                    swap(*this, tmp);
+                }
+            }
+        }
+
+        virtual constexpr
+        void push(const value_type& val)
+        {
+            if (base::end() == storage_end())
+                resize(base::p_capacity * 2);
+            *base::end() = val;
+            ++base::p_size;
+        }
+
+        virtual constexpr
+        void push(value_type&& val)
+        {
+            if (base::end() == storage_end())
+                resize(base::p_capacity * 2);
+            *base::end() = std::move(val);
+            ++base::p_size;
+        }
+
+        virtual constexpr
+        void push(const DynamicArray<T>& arr)
+        {
+            for (const value_type& el : arr)
+                push(el);
+        }
+
+        virtual constexpr
+        value_type pop()
+        {
+            if (is_empty())
+                throw hpr::LengthError("Cannot pop element from empty array");
+            value_type val = base::back();
+            std::destroy_at(base::p_start + base::p_size);
+            --base::p_size;
+            return val;
+        }
+
+        virtual constexpr
+        void insert(size_type position, const value_type& val)
+        {
+            if (base::end() == storage_end())
+                resize(base::p_capacity * 2);
+            for (size_type n = base::p_size; n > position; --n)
+                *(base::p_start + n) = std::move(*(base::p_start + n - 1));
+            *(base::p_start + position) = val;
+            ++base::p_size;
+        }
+
+        virtual constexpr
+        void insert(size_type position, value_type&& val)
+        {
+            if (base::end() == storage_end())
+                resize(base::p_capacity * 2);
+            for (size_type n = base::p_size; n > position; --n)
+                *(base::p_start + n) = std::move(*(base::p_start + n - 1));
+            *(base::p_start + position) = std::move(val);
+            ++base::p_size;
+        }
+
+        virtual constexpr
+        void remove(size_type position)
+        {
+            for (size_type n = position; n < base::p_size - 1; ++n)
+                *(base::p_start + n) = std::move(*(base::p_start + n + 1));
+            std::destroy_at(base::p_start + base::p_size);
+            --base::p_size;
+        }
+
+        virtual constexpr
+        void remove(iterator position)
+        {
+            if (position + 1 != base::end())
+                std::copy(position + 1, base::end(), position);
+            std::destroy_at(base::p_start + base::p_size);
+            --base::p_size;
+        }
+
+        virtual constexpr
+        void remove(const std::function<bool(value_type)>& condition)
+        {
+            size_type newSize = base::p_size;
+            for (size_type offset = 0; offset < newSize; ++offset)
+            {
+                if (condition(*(base::p_start + offset)))
+                {
+                    for (size_type n = offset; n < newSize; ++n)
+                        *(base::p_start + n) = std::move(*(base::p_start + n + 1));
+
+                    --newSize;
+                    --offset;
+                }
+            }
+            base::p_size = newSize;
+        }
+
+        virtual constexpr
+        void clear()
+        {
+            delete[] base::p_start;
+            base::p_size = 0;
+            base::p_capacity = 1;
+            base::p_start = new value_type[base::p_capacity];
+        }
+
+        DynamicArray slice(iterator start, iterator end)
+        {
+            return DynamicArray {start, end};
+        }
+
+    };
+
+} // end namespace hpr
\ No newline at end of file
diff --git a/source/hpr/container/iterator.hpp b/source/hpr/container/iterator.hpp
new file mode 100644
index 0000000..32db78a
--- /dev/null
+++ b/source/hpr/container/iterator.hpp
@@ -0,0 +1,98 @@
+#pragma once
+
+#include <iterator>
+
+
+namespace hpr
+{
+
+    template <typename Type, typename Category = std::forward_iterator_tag>
+    class Iterator
+    {
+
+    public:
+
+        using iterator_category = Category;
+        using difference_type = std::ptrdiff_t;
+        using value_type = Type;
+        using pointer = Type*;
+        using reference = Type&;
+        using iterator = Iterator;
+        using const_pointer = Type const*;
+        using const_reference = Type const&;
+        using const_iterator = Iterator const;
+
+    protected:
+
+        pointer p_ptr;
+
+    public:
+
+        Iterator() :
+                p_ptr {nullptr}
+        {}
+
+        Iterator(pointer ptr) :
+                p_ptr {ptr}
+        {}
+
+        reference operator*()
+        {
+            return *p_ptr;
+        }
+
+        const_reference operator*() const
+        {
+            return *p_ptr;
+        }
+
+        pointer operator->()
+        {
+            return p_ptr;
+        }
+
+        const_pointer operator->() const
+        {
+            return p_ptr;
+        }
+
+        iterator& operator++()
+        {
+            p_ptr++;
+            return *this;
+        }
+
+        const_iterator& operator++() const
+        {
+            p_ptr++;
+            return *this;
+        }
+
+        iterator operator++(int)
+        {
+            iterator temp {*this};
+            ++(*this);
+            return temp;
+        }
+
+        iterator operator+(int value)
+        {
+            iterator temp {*this};
+            temp.p_ptr += value;
+            return temp;
+        }
+
+        friend
+        bool operator==(const iterator& lhs, const iterator& rhs)
+        {
+            return lhs.p_ptr == rhs.p_ptr;
+        }
+
+        friend
+        bool operator!=(const iterator& lhs, const iterator& rhs)
+        {
+            return lhs.p_ptr != rhs.p_ptr;
+        }
+    };
+
+} // end namespace hpr
diff --git a/source/hpr/container/sequence.hpp b/source/hpr/container/sequence.hpp
new file mode 100644
index 0000000..fffac31
--- /dev/null
+++ b/source/hpr/container/sequence.hpp
@@ -0,0 +1,319 @@
+#pragma once
+
+#include <hpr/container/iterator.hpp>
+#include <hpr/exception.hpp>
+#include <hpr/numeric.hpp>
+
+
+namespace hpr
+{
+    // forward declaration
+
+    template <typename T>
+    class Sequence;
+
+    // type traits
+
+    template <typename T>
+    struct is_sequence : public std::false_type
+    {};
+
+    template <typename T>
+    struct is_sequence<Sequence<T>> : public std::true_type
+    {};
+
+    // concepts
+
+    template <typename Type>
+    concept IsSequence = is_sequence<Type>::value;
+
+    // class definition
+
+    template <typename T>
+    class Sequence
+    {
+
+    public:
+
+        using difference_type = std::ptrdiff_t;
+        using value_type = T;
+        using size_type = Size;
+        using pointer = T*;
+        using reference = T&;
+        using iterator = Iterator<T>;
+        using const_pointer = T const*;
+        using const_reference = T const&;
+        using const_iterator = Iterator<const T>;
+
+    protected:
+
+        size_type p_size;
+        size_type p_capacity;
+        pointer p_start;
+
+    protected:
+
+        constexpr
+        Sequence(size_type size, size_type capacity, pointer start) :
+            p_size {size},
+            p_capacity {capacity},
+            p_start {start}
+        {}
+
+        constexpr explicit
+        Sequence(size_type size, size_type capacity) :
+            p_size {size},
+            p_capacity {capacity},
+            p_start {p_capacity ? new value_type[p_capacity] : nullptr}
+        {}
+
+    public:
+
+        friend constexpr
+        void swap(Sequence& main, Sequence& other) noexcept
+        {
+            using std::swap;
+            swap(main.p_size, other.p_size);
+            swap(main.p_capacity, other.p_capacity);
+            swap(main.p_start, other.p_start);
+        }
+
+        constexpr
+        Sequence() noexcept :
+            p_size {0},
+            p_capacity {1},
+            p_start {new value_type[p_capacity]}
+        {}
+
+        constexpr
+        Sequence(const Sequence& seq) noexcept :
+            p_size {seq.p_size},
+            p_capacity {seq.p_capacity},
+            p_start {p_capacity ? new value_type[seq.p_capacity] : nullptr}
+        {
+            std::copy(seq.p_start, seq.p_start + seq.p_size, p_start);
+        }
+
+        //! Move constructor
+        constexpr
+        Sequence(Sequence&& seq) noexcept :
+            Sequence {}
+        {
+            swap(*this, seq);
+        }
+
+        template <typename Iter>
+        constexpr
+        Sequence(Iter start, Iter end) :
+            p_size {static_cast<size_type>(std::distance(start, end))},
+            p_capacity {p_size},
+            p_start {p_capacity ? new value_type[p_capacity] : nullptr}
+        {
+            std::copy(start, end, p_start);
+        }
+
+        constexpr
+        Sequence(iterator start, iterator end) :
+            p_size {static_cast<size_type>(std::distance(start, end))},
+            p_capacity {p_size},
+            p_start {p_capacity ? new value_type[p_capacity] : nullptr}
+        {
+            std::copy(start, end, p_start);
+        }
+
+        constexpr
+        Sequence(const_iterator start, const_iterator end) :
+            p_size {static_cast<size_type>(std::distance(start, end))},
+            p_capacity {p_size},
+            p_start {p_capacity ? new value_type[p_capacity] : nullptr}
+        {
+            std::copy(start, end, p_start);
+        }
+
+        constexpr
+        Sequence(iterator start, iterator end, size_type capacity) :
+            p_size {static_cast<size_type>(std::distance(start, end))},
+            p_capacity {capacity},
+            p_start {p_capacity ? new value_type[p_capacity] : nullptr}
+        {
+            std::copy(start, end, p_start);
+        }
+
+        constexpr
+        Sequence(const_iterator start, const_iterator end, size_type capacity) :
+            p_size {static_cast<size_type>(std::distance(start, end))},
+            p_capacity {capacity},
+            p_start {p_capacity ? new value_type[p_capacity] : nullptr}
+        {
+            std::copy(start, end, p_start);
+        }
+
+        /*constexpr
+        Sequence(std::initializer_list<value_type> list) :
+            Sequence {list.begin(), list.end()}
+        {}*/
+
+        constexpr
+        Sequence(size_type size, value_type value) noexcept:
+            p_size {size},
+            p_capacity {size},
+            p_start {p_capacity ? new value_type[p_capacity] : nullptr}
+        {
+            for (auto n = 0; n < size; ++n)
+                *(p_start + n) = value;
+        }
+
+        constexpr
+        Sequence& operator=(Sequence other) noexcept
+        {
+            swap(*this, other);
+            return *this;
+        }
+
+        constexpr
+        Sequence& operator=(Sequence&& seq) noexcept
+        {
+            swap(*this, seq);
+            return *this;
+        }
+
+        virtual
+        ~Sequence()
+        {
+            delete[] p_start;
+        }
+
+
+        // Member functions
+
+        virtual constexpr
+        iterator begin()
+        {
+            return iterator(p_start);
+        }
+
+        [[nodiscard]] virtual constexpr
+        const_iterator begin() const
+        {
+            return const_iterator(p_start);
+        }
+
+        virtual constexpr
+        iterator end()
+        {
+            return iterator(p_start + p_size);
+        }
+
+        virtual constexpr
+        const_iterator end() const
+        {
+            return const_iterator(p_start + p_size);
+        }
+
+        virtual constexpr
+        iterator storage_end()
+        {
+            return iterator(p_start + p_capacity);
+        }
+
+        virtual constexpr
+        const_iterator storage_end() const
+        {
+            return const_iterator(p_start + p_capacity);
+        }
+
+        [[nodiscard]] virtual constexpr
+        size_type size() const noexcept
+        {
+            return p_size;
+        }
+
+        virtual constexpr
+        reference operator[](size_type n)
+        {
+            if (n >= size())
+                throw hpr::OutOfRange();
+            return *(p_start + n);
+        }
+
+        virtual constexpr
+        const_reference operator[](size_type n) const
+        {
+            if (n >= size())
+                throw hpr::OutOfRange();
+            return *(p_start + n);
+        }
+
+        virtual constexpr
+        reference at(size_type n)
+        {
+            return operator[](n);
+        }
+
+        virtual constexpr
+        const_reference at(size_type n) const
+        {
+            return operator[](n);
+        }
+
+        virtual constexpr
+        reference front()
+        {
+            return *begin();
+        }
+
+        virtual constexpr
+        const_reference front() const
+        {
+            return *begin();
+        }
+
+        virtual constexpr
+        reference back()
+        {
+            return *(p_start + p_size - 1);
+        }
+
+        virtual constexpr
+        const_reference back() const
+        {
+            return *(p_start + p_size - 1);
+        }
+
+        virtual constexpr
+        pointer data()
+        {
+            return p_start;
+        }
+
+        [[nodiscard]] virtual constexpr
+        const_pointer data() const
+        {
+            return p_start;
+        }
+
+        // Friend functions
+
+        friend
+        bool operator==(const Sequence& lhs, const Sequence& rhs)
+        {
+            for (auto n = 0; n < lhs.size(); ++n)
+                if (lhs[n] != rhs[n])
+                    return false;
+            return true;
+        }
+
+        friend
+        Sequence operator+(const Sequence& lhs, const Sequence& rhs)
+        {
+            Sequence seq {rhs.size() + lhs.size(), {}};
+            for (auto n = 0; n < rhs.size(); ++n)
+            {
+                seq[n] = rhs[n];
+                seq[n + rhs.size()] = lhs[n];
+            }
+            return seq;
+        }
+    };
+
+} // end namespace hpr
\ No newline at end of file
diff --git a/source/hpr/container/static_array.hpp b/source/hpr/container/static_array.hpp
new file mode 100644
index 0000000..cd58a02
--- /dev/null
+++ b/source/hpr/container/static_array.hpp
@@ -0,0 +1,146 @@
+#pragma once
+
+#include <hpr/container/sequence.hpp>
+
+
+namespace hpr
+{
+    // forward declaration
+
+    template <typename T, Size S>
+    requires (S > 0)
+    class StaticArray;
+
+    // type traits
+
+    template <typename T, Size S>
+    struct is_sequence<StaticArray<T, S>> : public std::true_type
+    {};
+
+    // aliases
+
+    template <typename T, Size S>
+    using sarray = StaticArray<T, S>;
+
+    // class definition
+
+    template <typename T, Size S>
+    requires (S > 0)
+    class StaticArray : public Sequence<T>
+    {
+
+        using base = Sequence<T>;
+
+    public:
+
+        using difference_type = std::ptrdiff_t;
+        using value_type = T;
+        using size_type = Size;
+        using pointer = T*;
+        using reference = T&;
+        using iterator = Iterator<T>;
+        using const_pointer = T const*;
+        using const_reference = T const&;
+        using const_iterator = Iterator<const T>;
+
+    public:
+
+        friend constexpr
+        void swap(StaticArray& main, StaticArray& other) noexcept
+        {
+            using std::swap;
+            swap(static_cast<base&>(main), static_cast<base&>(other));
+        }
+
+        //! Default constructor
+        constexpr
+        StaticArray() noexcept :
+            base {S, S, new value_type[S] {}}
+        {}
+
+        constexpr explicit
+        StaticArray(const base& b) noexcept :
+            base {b}
+        {}
+
+        //! Copy constructor
+        constexpr
+        StaticArray(const StaticArray& arr) noexcept :
+            base {static_cast<base>(arr)}
+        {}
+
+        //! Move constructor
+        constexpr
+        StaticArray(StaticArray&& arr) noexcept :
+            base {std::forward<base>(static_cast<base>(arr))}
+        {}
+
+        constexpr
+        StaticArray(typename base::iterator start, typename base::iterator end) :
+            base {start, end}
+        {}
+
+        constexpr
+        StaticArray(typename base::const_iterator start, typename base::const_iterator end) :
+            base {start, end}
+        {}
+
+        constexpr
+        StaticArray(typename base::iterator start, typename base::iterator end, size_type capacity) :
+            base {start, end, capacity}
+        {}
+
+        constexpr
+        StaticArray(typename base::const_iterator start, typename base::const_iterator end, size_type capacity) :
+            base {start, end, capacity}
+        {}
+
+        constexpr
+        StaticArray(std::initializer_list<value_type> list) :
+            base {list.begin(), list.end()}
+        {}
+
+        template <std::convertible_to<value_type>... Args>
+        constexpr explicit
+        StaticArray(const value_type& v, const Args& ...args)
+        requires (1 + sizeof...(args) == S) :
+            StaticArray {std::initializer_list<value_type>({v, static_cast<value_type>(args)...})}
+        {}
+
+        template <std::convertible_to<value_type>... Args>
+        constexpr explicit
+        StaticArray(value_type&& v, Args&& ...args)
+        requires (1 + sizeof...(args) == S) :
+            StaticArray {std::initializer_list<value_type>({std::forward<value_type>(v), std::forward<value_type>(static_cast<value_type>(args))...})}
+        {}
+
+        template <size_type SS, std::convertible_to<value_type>... Args>
+        constexpr
+        StaticArray(const StaticArray<value_type, SS>& subArr, const value_type& v, const Args& ...args) noexcept
+        requires (SS + 1 + sizeof...(args) == S) :
+            base {S, S, new value_type[S] {}}
+        {
+            std::copy(subArr.begin(), subArr.end(), base::begin());
+            std::initializer_list<value_type> list {v, static_cast<value_type>(args)...};
+            std::copy(list.begin(), list.end(), base::begin() + subArr.size());
+        }
+
+        constexpr
+        StaticArray& operator=(StaticArray other)
+        {
+            swap(*this, other);
+            return *this;
+        }
+
+        constexpr
+        StaticArray& operator=(StaticArray&& other) noexcept
+        {
+            swap(*this, other);
+            return *this;
+        }
+
+        virtual
+        ~StaticArray() = default;
+    };
+
+} // end namespace hpr
diff --git a/source/hpr/container/tree_node.hpp b/source/hpr/container/tree_node.hpp
new file mode 100644
index 0000000..af3d1b0
--- /dev/null
+++ b/source/hpr/container/tree_node.hpp
@@ -0,0 +1,181 @@
+#pragma once
+
+#include <hpr/container/dynamic_array.hpp>
+
+#include <memory>
+
+
+namespace hpr
+{
+
+    template <typename Type>
+    class TreeNode
+    {
+
+    public:
+
+        using value_type = Type;
+        using pointer = Type*;
+        using reference = Type&;
+        using const_pointer = Type const*;
+        using const_reference = Type const&;
+
+    protected:
+
+        pointer p_data;
+
+        darray<TreeNode*> p_descendants;
+        TreeNode* p_ancestor;
+
+    public:
+
+        friend constexpr
+        void swap(TreeNode& main, TreeNode& other) noexcept
+        {
+            using std::swap;
+            swap(main.p_data, other.p_data);
+            swap(main.p_descendants, other.p_descendants);
+            swap(main.p_ancestor, other.p_ancestor);
+        }
+
+        inline
+        TreeNode() :
+            p_data {new value_type {}},
+            p_descendants {},
+            p_ancestor {}
+        {}
+
+        inline
+        TreeNode(const TreeNode& node) :
+            p_data {new value_type {*node.p_data}},
+            p_descendants {},
+            p_ancestor {}
+        {}
+
+        inline
+        TreeNode(TreeNode&& node)  noexcept :
+            p_data {},
+            p_descendants {},
+            p_ancestor {}
+        {
+            swap(*this, node);
+            /*std::swap(p_data, node.p_data);
+            std::swap(p_descendants, node.p_descendants);
+            std::swap(p_ancestor, node.p_ancestor);*/
+        }
+
+        inline
+        TreeNode& operator=(TreeNode other)
+        {
+            swap(*this, other);
+            return *this;
+        }
+
+        inline explicit
+        TreeNode(const value_type& data) :
+            p_data {new value_type {data}},
+            p_descendants {},
+            p_ancestor {}
+        {}
+
+        inline
+        TreeNode(const value_type& data, const darray<TreeNode*>& descendants, TreeNode* ancestor) :
+            p_data {new value_type {data}},
+            p_descendants {descendants},
+            p_ancestor {ancestor}
+        {
+            for (auto descendant : p_descendants)
+                descendant->ancestor(this);
+        }
+
+        inline
+        TreeNode(const value_type& data, const darray<TreeNode*>& descendants) :
+            p_data {new value_type {data}},
+            p_descendants {descendants},
+            p_ancestor {}
+        {
+            for (auto descendant : p_descendants)
+                descendant->ancestor(this);
+        }
+
+        virtual
+        ~TreeNode()
+        {
+            delete p_data;
+            //delete p_ancestor;
+        }
+
+        // Member functions
+
+        inline
+        pointer data()
+        {
+            return p_data;
+        }
+
+        inline
+        void ancestor(TreeNode* node)
+        {
+            p_ancestor = node;
+        }
+
+        inline
+        TreeNode* ancestor()
+        {
+            return p_ancestor;
+        }
+
+        inline
+        void descendants(const darray<TreeNode*>& descendants)
+        {
+            p_descendants = descendants;
+        }
+
+        inline
+        darray<TreeNode*>& descendants()
+        {
+            return p_descendants;
+        }
+
+        darray<TreeNode*> traverse_descendants()
+        {
+            std::function<darray<TreeNode*>(TreeNode*)> collect = [&](TreeNode* node)
+            {
+                darray<TreeNode*> ds;
+
+                if (!node->descendants().is_empty())
+                {
+                    for (TreeNode* dnode: node->descendants())
+                    {
+                        ds.push(dnode);
+                        if (!dnode->descendants().is_empty())
+                            ds.push(collect(dnode));
+                    }
+                }
+
+                return ds;
+            };
+
+            return collect(this);
+        }
+
+        darray<TreeNode*> traverse_ancestors()
+        {
+            std::function<darray<TreeNode*>(TreeNode*)> collect = [&](TreeNode* node)
+            {
+                darray<TreeNode*> ds;
+
+                if (node->p_ancestor)
+                {
+                    ds.push(node);
+                    ds.push(collect(node->p_ancestor));
+                }
+
+                return ds;
+            };
+
+            return collect(this);
+        }
+    };
+
+} // end namespace hpr
\ No newline at end of file
diff --git a/source/hpr/exception.hpp b/source/hpr/exception.hpp
new file mode 100644
index 0000000..41ad7e7
--- /dev/null
+++ b/source/hpr/exception.hpp
@@ -0,0 +1,63 @@
+#pragma once
+
+#include <exception>
+#include <string>
+#include <source_location>
+#include <sstream>
+#include <utility>
+#include <vector>
+
+
+namespace hpr
+{
+    class Exception : public std::exception
+    {
+
+    protected:
+
+        std::string p_message;
+
+    public:
+
+        inline explicit
+        Exception(std::source_location location = std::source_location::current()) :
+                p_message {}
+        {
+            std::stringstream _message;
+            _message << "\t" << p_message
+                     << "\n  where:\t\t" << location.file_name() << ":" << location.line() << ":" << location.column()
+                     << "\n  function:\t" <<  location.function_name();
+            p_message = _message.str();
+        }
+
+        inline explicit
+        Exception(const std::string& message, std::source_location location = std::source_location::current()) :
+                p_message {message}
+        {
+            std::stringstream _message;
+            _message << "\t" << p_message
+                     << "\n  where:\t" << location.file_name() << ":" << location.line() << ":" << location.column()
+                     << "\n  function:\t" <<  location.function_name();
+            p_message = _message.str();
+        }
+
+        [[nodiscard]]
+        const char* what() const noexcept override {
+                std::vector<float> vv;
+                return p_message.data();
+        }
+    };
+
+    struct OutOfRange : public Exception
+    {
+        inline explicit OutOfRange(std::source_location location = std::source_location::current()) : Exception {"Out of range", location} {}
+        inline explicit OutOfRange(const std::string& message, std::source_location location = std::source_location::current()) : Exception {message, location} {}
+    };
+
+    struct LengthError : public Exception
+    {
+        inline explicit LengthError(std::source_location location = std::source_location::current()) : Exception {"Length error", location} {}
+        inline explicit LengthError(const std::string& message, std::source_location location = std::source_location::current()) : Exception {message, location} {}
+    };
+
+} // end namespace hpr
\ No newline at end of file
diff --git a/source/hpr/hpr.cpp b/source/hpr/hpr.cpp
new file mode 100644
index 0000000..1e1dd56
--- /dev/null
+++ b/source/hpr/hpr.cpp
@@ -0,0 +1,5 @@
+
+#include <hpr/exception.hpp>
+#include <hpr/numeric.hpp>
+#include <hpr/container.hpp>
+#include <hpr/math.hpp>
diff --git a/source/hpr/main.cpp b/source/hpr/main.cpp
deleted file mode 100644
index db5efd2..0000000
--- a/source/hpr/main.cpp
+++ /dev/null
@@ -1,9 +0,0 @@
-#include <iostream>
-
-using namespace std;
-
-int main(int argc, char** argv)
-{
-    cout << "hello xmake!" << endl;
-    return 0;
-}
diff --git a/source/hpr/math.hpp b/source/hpr/math.hpp
new file mode 100644
index 0000000..13f8985
--- /dev/null
+++ b/source/hpr/math.hpp
@@ -0,0 +1,5 @@
+#pragma once
+
+#include <hpr/math/vector.hpp>
+#include <hpr/math/quaternion.hpp>
+#include <hpr/math/matrix.hpp>
diff --git a/source/hpr/math/matrix.hpp b/source/hpr/math/matrix.hpp
new file mode 100644
index 0000000..5ccd103
--- /dev/null
+++ b/source/hpr/math/matrix.hpp
@@ -0,0 +1,509 @@
+#pragma once
+
+#include <hpr/numeric.hpp>
+#include <hpr/container/static_array.hpp>
+#include <hpr/math/quaternion.hpp>
+
+
+namespace hpr
+{
+    // forward declarations
+
+    template <IsReal Type, Size Rows, Size Cols> requires (Rows >= 0 && Cols >= 0)
+    class Matrix;
+
+    template <IsReal Type, Size Rows, Size Cols>
+    using SubMatrix = typename std::conditional<(Rows >= 2 && Cols >= 2), Matrix<Type, Rows - 1, Cols - 1>, Matrix<Type, 1, 1>>::type;
+
+    // type traits
+
+    template <typename T>
+    struct is_matrix : public std::false_type {};
+
+    template <typename T, Size Rows, Size Cols>
+    struct is_matrix<Matrix<T, Rows, Cols>> : public std::true_type {};
+
+    // concepts
+
+    template <typename T>
+    concept IsMatrix = is_matrix<T>::value;
+
+    // aliases
+
+    template <typename Type, size_t Row, size_t Col>
+    using mat = Matrix<Type, Row, Col>;
+
+    using mat2 = Matrix<scalar, 2, 2>;
+    using mat3 = Matrix<scalar, 3, 3>;
+    using mat4 = Matrix<scalar, 4, 4>;
+
+
+    template <IsReal Type, Size Rows, Size Cols> requires (Rows >= 0 && Cols >= 0)
+    class Matrix : public StaticArray<Type, Rows * Cols>
+    {
+
+        using base = StaticArray<Type, Rows * Cols>;
+
+    public:
+
+        using value_type = Type;
+        using size_type = Size;
+        using pointer = Type*;
+        using reference = Type&;
+        using iterator = Iterator<Type>;
+        using const_reference = Type const&;
+        using const_iterator = Iterator<const Type>;
+
+    protected:
+
+        size_type p_rows;
+        size_type p_cols;
+
+    public:
+
+        friend constexpr
+        void swap(Matrix& main, Matrix& other)
+        {
+            using std::swap;
+            swap(static_cast<base&>(main), static_cast<base&>(other));
+            swap(main.p_rows, other.p_rows);
+            swap(main.p_cols, other.p_cols);
+        }
+
+        inline
+        Matrix() :
+            base {},
+            p_rows {Rows},
+            p_cols {Cols}
+        {}
+
+        inline
+        Matrix(const Matrix& ms) :
+            base {static_cast<base>(ms)},
+            p_rows {Rows},
+            p_cols {Cols}
+        {}
+
+        inline
+        Matrix(Matrix&& ms) noexcept:
+            base {std::forward<base>(static_cast<base>(ms))},
+            p_rows {Rows},
+            p_cols {Cols}
+        {}
+
+        inline
+        Matrix& operator=(const Matrix& ms)
+        {
+            //base::operator=(ms);
+            swap(*this, ms);
+            return *this;
+        }
+
+        inline explicit
+        Matrix(const base& vs) :
+                base {vs},
+                p_rows {Rows},
+                p_cols {Cols}
+        {}
+
+        inline explicit
+        Matrix(base&& vs) noexcept:
+            base {std::forward<base>(vs)},
+            p_rows {Rows},
+            p_cols {Cols}
+        {}
+
+        inline
+        Matrix(typename base::iterator start, typename base::iterator end) :
+            base {start, end},
+            p_rows {Rows},
+            p_cols {Cols}
+        {}
+
+        inline
+        Matrix(typename base::const_iterator start, typename base::const_iterator end) :
+            base {start, end},
+            p_rows {Rows},
+            p_cols {Cols}
+        {}
+
+        inline
+        Matrix(std::initializer_list<value_type> list) :
+            base {list},
+            p_rows {Rows},
+            p_cols {Cols}
+        {}
+
+        template <IsReal... Args>
+        inline
+        Matrix(value_type&& v, Args&& ...args) requires (1 + sizeof...(args) == Rows * Cols):
+            base {v, static_cast<value_type>(std::forward<Args>(args))...},
+            p_rows {Rows},
+            p_cols {Cols}
+        {}
+
+        inline
+        Matrix(const value_type& v) :
+            base {},
+            p_rows {Rows},
+            p_cols {Cols}
+        {
+            for (Size n = 0; n < Rows * Cols; ++n)
+                (*this)[n] = v;
+        }
+
+        inline
+        Matrix& operator=(const value_type& v)
+        {
+            for (Size n = 0; n < Rows * Cols; ++n)
+                (*this)[n] = v;
+            return *this;
+        }
+
+        inline explicit
+        Matrix(const Quaternion& q) requires (Rows == 3 && Cols == 3 || Rows == 4 && Cols == 4) :
+            base {},
+            p_rows {Rows},
+            p_cols {Cols}
+        {
+            const scalar s = pow(norm(q), -2);
+
+            (*this)(0, 0) = 1 - 2 * s * (q[2] * q[2] + q[3] * q[3]);
+            (*this)(1, 0) = 2 * s * (q[1] * q[2] - q[3] * q[0]);
+            (*this)(2, 0) = 2 * s * (q[1] * q[3] - q[2] * q[0]);
+
+            (*this)(0, 1) = 2 * s * (q[1] * q[2] + q[3] * q[0]);
+            (*this)(1, 1) = 1 - 2 * s * (q[1] * q[1] + q[3] * q[3]);
+            (*this)(2, 1) = 2 * s * (q[2] * q[3] + q[1] * q[0]);
+
+            (*this)(0, 2) = 2 * s * (q[1] * q[3] + q[2] * q[0]);
+            (*this)(1, 2) = 2 * s * (q[2] * q[3] + q[1] * q[0]);
+            (*this)(2, 2) = 1 - 2 * s * (q[1] * q[1] + q[2] * q[2]);
+
+            if constexpr (Rows == 4)
+                (*this)(3, 3) = 1;
+        }
+
+
+        // access
+
+        inline
+        reference operator()(size_type row, size_type col)
+        {
+            if (row >= p_rows || std::numeric_limits<size_type>::max() - p_rows < row)
+                throw std::out_of_range("Row index is out of range");
+            if (col >= p_cols || std::numeric_limits<size_type>::max() - p_cols < col)
+                throw std::out_of_range("Column index is out of range");
+            return (*this)[col + p_rows * row];
+        }
+
+        inline
+        const_reference operator()(size_type row, size_type col) const
+        {
+            if (row >= p_rows || std::numeric_limits<size_type>::max() - p_rows < row)
+                throw std::out_of_range("Row index is out of range");
+            if (col >= p_cols || std::numeric_limits<size_type>::max() - p_cols < col)
+                throw std::out_of_range("Column index is out of range");
+            return (*this)[col + p_rows * row];
+        }
+
+        Vector<value_type, Cols> row(size_type row) const
+        {
+            Vector<value_type, Cols> vs;
+            for (auto n = 0; n < Cols; ++n)
+                vs[n] = (*this)(row, n);
+            return vs;
+        }
+
+        void row(size_type row, const Vector<value_type, Cols>& vs)
+        {
+            for (auto n = 0; n < Cols; ++n)
+                (*this)(row, n) = vs[n];
+        }
+
+        Vector<value_type, Rows> col(size_type col) const
+        {
+            Vector<value_type, Rows> vs;
+            for (auto n = 0; n < Rows; ++n)
+                vs[n] = (*this)(n, col);
+            return vs;
+        }
+
+        void col(size_type col, const Vector<value_type, Rows>& vs)
+        {
+            for (auto n = 0; n < Rows; ++n)
+                (*this)(n, col) = vs[n];
+        }
+
+        [[nodiscard]] constexpr size_type rows() const { return p_rows; }
+        [[nodiscard]] constexpr size_type cols() const { return p_cols; }
+
+        // member functions
+
+        [[nodiscard]]
+        constexpr
+        bool is_square() const
+        {
+            return p_rows == p_cols;
+        }
+
+        inline
+        Matrix& fill(value_type value)
+        {
+            for (auto n = 0; n < this->size(); ++n)
+                (*this)[n] = value;
+            return *this;
+        }
+
+        // Global functions
+
+        static inline
+        Matrix identity()
+        {
+            Matrix ms;
+            for (auto n = 0; n < Rows; ++n)
+                //for (auto k = 0; k < Cols; ++k)
+                    ms(n, n) = 1;
+            return ms;
+        }
+    };
+
+    // global operators
+
+    template <IsReal T, Size R, Size C> inline Matrix<T, R, C> operator+(const Matrix<T, R, C>& lhs) { Matrix<T, R, C> ms; for (Size n = 0; n < lhs.size(); ++n) ms[n] = lhs[n]; return ms; }
+    template <IsReal T, Size R, Size C> inline Matrix<T, R, C> operator-(const Matrix<T, R, C>& lhs) { Matrix<T, R, C> ms; for (Size n = 0; n < lhs.size(); ++n) ms[n] = -lhs[n]; return ms; }
+
+    template <IsReal T, Size R, Size C> inline Matrix<T, R, C>& operator+=(Matrix<T, R, C>& lhs, const Matrix<T, R, C>& rhs) { for (Size n = 0; n < lhs.size(); ++n) lhs[n] += rhs[n]; return lhs; }
+    template <IsReal T, Size R, Size C> inline Matrix<T, R, C>& operator-=(Matrix<T, R, C>& lhs, const Matrix<T, R, C>& rhs) { for (Size n = 0; n < lhs.size(); ++n) lhs[n] -= rhs[n]; return lhs; }
+    template <IsReal T, Size R, Size C> inline Matrix<T, R, C>& operator*=(Matrix<T, R, C>& lhs, const Matrix<T, C, R>& rhs) { Matrix<T, R, C> temp {lhs}; for (Size n = 0; n < R; ++n) for (Size k = 0; k < C; ++k) lhs(n, k) = sum(temp.col(k) * rhs.row(n)); return lhs; }
+
+    template <IsReal T, Size R, Size C> inline Matrix<T, R, C> operator+(const Matrix<T, R, C>& lhs, const Matrix<T, R, C>& rhs) { Matrix<T, R, C> ms {lhs}; for (Size n = 0; n < lhs.size(); ++n) ms[n] += rhs[n]; return ms; }
+    template <IsReal T, Size R, Size C> inline Matrix<T, R, C> operator-(const Matrix<T, R, C>& lhs, const Matrix<T, R, C>& rhs) { Matrix<T, R, C> ms {lhs}; for (Size n = 0; n < lhs.size(); ++n) ms[n] -= rhs[n]; return ms; }
+    template <IsReal T, Size R, Size C> inline Matrix<T, R, C> operator*(const Matrix<T, R, C>& lhs, const Matrix<T, C, R>& rhs) { Matrix<T, R, C> ms; for (Size n = 0; n < R; ++n) for (Size k = 0; k < C; ++k) ms(n, k) = sum(lhs.col(k) * rhs.row(n)); return ms; }
+
+    template <IsReal T, Size R, Size C> inline bool operator==(const Matrix<T, R, C>& lhs, const Matrix<T, R, C>& rhs) { for (Size n = 0; n < lhs.size(); ++n) if (lhs[n] != rhs[n]) return false; return true; }
+    template <IsReal T, Size R, Size C> inline bool operator!=(const Matrix<T, R, C>& lhs, const Matrix<T, R, C>& rhs) { for (Size n = 0; n < lhs.size(); ++n) if (lhs[n] == rhs[n]) return false; return true; }
+
+
+    template <IsReal T, Size R, Size C> inline Matrix<T, R, C>& operator+=(Matrix<T, R, C>& lhs, const T& rhs) { for (Size n = 0; n < lhs.size(); ++n) lhs[n] += rhs; return lhs; }
+    template <IsReal T, Size R, Size C> inline Matrix<T, R, C>& operator-=(Matrix<T, R, C>& lhs, const T& rhs) { for (Size n = 0; n < lhs.size(); ++n) lhs[n] -= rhs; return lhs; }
+    template <IsReal T, Size R, Size C> inline Matrix<T, R, C>& operator*=(Matrix<T, R, C>& lhs, const T& rhs) { for (Size n = 0; n < lhs.size(); ++n) lhs[n] *= rhs; return lhs; }
+    template <IsReal T, Size R, Size C> inline Matrix<T, R, C>& operator/=(Matrix<T, R, C>& lhs, const T& rhs) { for (Size n = 0; n < lhs.size(); ++n) lhs[n] /= rhs; return lhs; }
+
+    template <IsReal T, Size R, Size C> inline Matrix<T, R, C> operator+(const Matrix<T, R, C>& lhs, const T& rhs) { Matrix<T, R, C> ms {lhs}; for (Size n = 0; n < lhs.size(); ++n) ms[n] += rhs; return ms; }
+    template <IsReal T, Size R, Size C> inline Matrix<T, R, C> operator-(const Matrix<T, R, C>& lhs, const T& rhs) { Matrix<T, R, C> ms {lhs}; for (Size n = 0; n < lhs.size(); ++n) ms[n] -= rhs; return ms; }
+    template <IsReal T, Size R, Size C> inline Matrix<T, R, C> operator*(const Matrix<T, R, C>& lhs, const T& rhs) { Matrix<T, R, C> ms {lhs}; for (Size n = 0; n < lhs.size(); ++n) ms[n] *= rhs; return ms; }
+    template <IsReal T, Size R, Size C> inline Matrix<T, R, C> operator/(const Matrix<T, R, C>& lhs, const T& rhs) { Matrix<T, R, C> ms {lhs}; for (Size n = 0; n < lhs.size(); ++n) ms[n] /= rhs; return ms; }
+
+
+    template <IsReal T, Size R, Size C> inline Vector<T, R> operator*(const Matrix<T, R, C>& ms, const Vector<T, R>& vs) { Vector<T, R> res; for (Size n = 0; n < R; ++n) res[n] = sum(ms.row(n) * vs); return res; }
+    template <IsReal T, Size R, Size C> inline Vector<T, C> operator*(const Vector<T, R>& vs, const Matrix<T, R, C>& ms) { Vector<T, C> res; for (Size n = 0; n < C; ++n) res[n] = sum(ms.col(n) * vs); return res; }
+
+    template <IsReal T, Size R, Size C> inline bool operator==(const Matrix<T, R, C>& lhs, const Vector<T, R * C>& rhs) { return false; }
+    template <IsReal T, Size R, Size C> inline bool operator!=(const Matrix<T, R, C>& lhs, const Vector<T, R * C>& rhs) { return true; }
+
+    // matrix operations
+
+    //! Transpose matrix
+    template <IsReal T, Size R, Size C>
+    inline
+    Matrix<T, R, C> transpose(const Matrix<T, R, C>& ms)
+    {
+        Matrix<T, R, C> res;
+        for (Size n = 0; n < R; ++n)
+            for (Size k = 0; k < C; ++k)
+                res(k, n) = ms(n, k);
+        return res;
+    }
+
+    //! Trace of a matrix
+    template <IsReal T, Size R, Size C>
+    inline
+    T trace(const Matrix<T, R, C>& ms) requires (R == C)
+    {
+        T res;
+        for (auto n = 0; n < R; ++n)
+            res += ms(n, n);
+        return res;
+    }
+
+    //! Minor of a matrix
+    template <IsReal T, Size R, Size C>
+    inline
+    SubMatrix<T, R, C> minor(const Matrix<T, R, C>& ms, Size row, Size col)
+    {
+        if (ms.size() < 4)
+            throw std::runtime_error("Matrix should be greater 2x2");
+
+        SubMatrix<T, R, C> minor;
+        auto minor_iter = minor.begin();
+        for (auto n = 0; n < R; ++n)
+            for (auto k = 0; k < C; ++k)
+                if (k != col && n != row)
+                    *(minor_iter++) = ms[k + ms.rows() * n];
+        return minor;
+    }
+
+    //! Determinant of a matrix
+    template <IsReal T, Size R, Size C>
+    inline
+    scalar det(const Matrix<T, R, C>& ms) requires (R == C)
+    {
+        if (ms.size() == 1)
+            return ms[0];
+
+        else if (ms.size() == 4)
+            return ms(0, 0) * ms(1, 1) - ms(0, 1) * ms(1, 0);
+
+        else {
+            scalar res = 0;
+            for (auto n = 0; n < ms.cols(); ++n)
+                res += pow(-1, n) * ms(0, n) * det(minor(ms, 0, n));
+            return res;
+        }
+    }
+
+    //! Adjoint matrix
+    template <IsReal T, Size R, Size C>
+    inline
+    Matrix<T, R, C> adj(const Matrix<T, R, C>& ms)
+    {
+        Matrix<T, R, C> res;
+        for (auto n = 0; n < R; ++n)
+            for (auto k = 0; k < C; ++k)
+                res(n, k) = pow(-1, n + k) * det(minor(ms, n, k));
+        return transpose(res);
+    }
+
+    //! Inverse matrix
+    template <IsReal T, Size R, Size C>
+    inline
+    Matrix<T, R, C> inv(const Matrix<T, R, C>& ms)
+    {
+        return adj(ms) / det(ms);
+    }
+
+    // Transforms
+
+    template <IsReal T>
+    inline
+    Matrix<T, 4, 4> translate(const Matrix<T, 4, 4>& ms, const Vector<T, 3>& vs)
+    {
+        Matrix<T, 4, 4> res {ms};
+        res.col(3, ms.row(0) * vs[0] + ms.row(1) * vs[1] + ms.row(2) * vs[2] + ms.row(3));
+
+        return res;
+    }
+
+    template <IsReal T>
+    inline
+    Matrix<T, 4, 4> rotate(const Matrix<T, 4, 4>& ms, const Vector<T, 3>& vs, T angle)
+    {
+        const T cosv = cos(angle);
+        const T sinv = sin(angle);
+        Vector<T, 3> axis {normalize(vs)};
+        Vector<T, 3> temp {(static_cast<T>(1) - cosv) * axis};
+
+        Matrix<T, 4, 4> rot;
+        rot(0, 0) = cosv + temp[0] * axis[0];
+        rot(0, 1) = temp[0] * axis[1] + sinv * axis[2];
+        rot(0, 2) = temp[0] * axis[2] - sinv * axis[1];
+        rot(1, 0) = temp[1] * axis[0] - sinv * axis[2];
+        rot(1, 1) = cosv + temp[1] * axis[1];
+        rot(1, 2) = temp[1] * axis[2] + sinv * axis[0];
+        rot(2, 0) = temp[2] * axis[0] + sinv * axis[1];
+        rot(2, 1) = temp[2] * axis[1] - sinv * axis[0];
+        rot(2, 2) = cosv + temp[2] * axis[2];
+
+        Matrix<T, 4, 4> res {ms};
+        res.row(0, ms.row(0) * rot(0, 0) + ms.row(1) * rot(0, 1) + ms.row(2) * rot(0, 2));
+        res.row(1, ms.row(0) * rot(1, 0) + ms.row(1) * rot(1, 1) + ms.row(2) * rot(1, 2));
+        res.row(2, ms.row(0) * rot(2, 0) + ms.row(1) * rot(2, 1) + ms.row(2) * rot(2, 2));
+        res.row(3, ms.row(3));
+
+        return res;
+    }
+
+    template <IsReal T>
+    inline
+    Matrix<T, 4, 4> rotate(const Matrix<T, 4, 4>& ms, const Quaternion& q)
+    {
+        return ms * Matrix<T, 4, 4>(q);
+    }
+
+    template <IsReal T>
+    inline
+    Matrix<T, 4, 4> scale(const Matrix<T, 4, 4>& ms, const Vector<T, 3>& vs)
+    {
+        Matrix<T, 4, 4> res;
+        res.row(0, ms.row(0) * vs[0]);
+        res.row(1, ms.row(1) * vs[1]);
+        res.row(2, ms.row(2) * vs[2]);
+        res.row(3, ms.row(3));
+
+        return res;
+    }
+
+    template <typename T>
+    inline
+    Matrix<T, 4, 4> lookAt(const Vector<T, 3>& eye, const Vector<T, 3>& center, const Vector<T, 3>& up)
+    {
+        const Vector<T, 3> forward {normalize(center - eye)};
+        const Vector<T, 3> right {normalize(cross(forward, up))};
+        const Vector<T, 3> nup {cross(right, forward)};
+        const Vector<T, 3> translation {dot(right, eye), dot(nup, eye), -dot(forward, eye)};
+
+        Matrix<T, 4, 4> res = Matrix<T, 4, 4>::identity();
+        res.row(0, Vector<T, 4>(right, 0));
+        res.row(1, Vector<T, 4>(nup, 0));
+        res.row(2, Vector<T, 4>(-forward, 0));
+        res.col(3, Vector<T, 4>(-translation, static_cast<T>(1)));
+
+        return res;
+    }
+
+    // Clip space
+
+    template <IsReal T>
+    inline
+    Matrix<T, 4, 4> ortho(T left, T right, T bottom, T top)
+    {
+        Matrix<T, 4, 4> ms = Matrix<T, 4, 4>::identity();
+        ms(0, 0) = static_cast<T>(2) / (right - left);
+        ms(1, 1) = static_cast<T>(2) / (top - bottom);
+        ms(2, 2) = -static_cast<T>(1);
+        ms(3, 0) = -(right + left) / (right - left);
+        ms(3, 1) = -(top + bottom) / (top - bottom);
+        return ms;
+    }
+
+    template <IsReal T>
+    inline
+    Matrix<T, 4, 4> ortho(T left, T right, T bottom, T top, T zNear, T zFar)
+    {
+        Matrix<T, 4, 4> ms = Matrix<T, 4, 4>::identity();
+        ms(0, 0) = static_cast<T>(2) / (right - left);
+        ms(1, 1) = static_cast<T>(2) / (top - bottom);
+        ms(2, 2) = -static_cast<T>(2) / (zFar - zNear);
+        ms(0, 3) = -(right + left) / (right - left);
+        ms(1, 3) = -(top + bottom) / (top - bottom);
+        ms(2, 3) = -(zFar + zNear) / (zFar - zNear);
+        return ms;
+    }
+
+    template <IsReal T>
+    inline
+    Matrix<T, 4, 4> perspective(T fovy, T aspect, T zNear, T zFar)
+    {
+        assert(abs(aspect - std::numeric_limits<T>::epsilon()) > 0);
+        Matrix<T, 4, 4> ms;
+        const T halfFovyTan = tan(fovy / 2);
+        ms(0, 0) = static_cast<T>(1) / (aspect * halfFovyTan);
+        ms(1, 1) = static_cast<T>(1) / halfFovyTan;
+        ms(2, 2) = -(zFar + zNear) / (zFar - zNear);
+        ms(3, 2) = -static_cast<T>(1);
+        ms(2, 3) = -(static_cast<T>(2) * zFar * zNear) / (zFar - zNear);
+        return ms;
+    }
+
+} // end namespace hpr
\ No newline at end of file
diff --git a/source/hpr/math/quaternion.hpp b/source/hpr/math/quaternion.hpp
new file mode 100644
index 0000000..e987bc2
--- /dev/null
+++ b/source/hpr/math/quaternion.hpp
@@ -0,0 +1,290 @@
+#pragma once
+
+#include <hpr/math/vector.hpp>
+#include <hpr/exception.hpp>
+
+
+namespace hpr
+{
+    // forward declarations
+
+    class Quaternion;
+
+    inline
+    Quaternion inverse(const Quaternion& q);
+
+    inline
+    Quaternion operator*(const Quaternion& lhs, const Quaternion& rhs);
+
+    // aliases
+
+    using quat = Quaternion;
+
+    // class declaration
+    class Quaternion
+    {
+
+    public:
+
+        enum RotationSequence
+        {
+            ZYX, ZYZ, ZXY, ZXZ, YXZ, YXY, YZX, YZY, XYZ, XYX, XZY, XZX
+        };
+
+    protected:
+
+        scalar p_real;
+        vec3 p_imag;
+
+    public:
+
+        inline
+        Quaternion() :
+            p_real{},
+            p_imag{}
+        {}
+
+        inline
+        Quaternion(const scalar real, const vec3& imag) :
+            p_real {real},
+            p_imag {imag}
+        {}
+
+        inline explicit
+        Quaternion(const scalar real) :
+            p_real {real},
+            p_imag {}
+        {}
+
+        inline explicit
+        Quaternion(const vec3& imag) :
+            p_real {},
+            p_imag {imag}
+        {}
+
+        inline
+        Quaternion(const vec3& vs, const scalar& theta) :
+            p_real {cos(0.5 * theta)},
+            p_imag {sin(0.5 * theta) * vs / mag(vs)}
+        {}
+
+        static inline
+        Quaternion unit(const vec3& vs)
+        {
+            return Quaternion(sqrt(1 - norm(vs)), vs);
+        }
+
+        inline
+        Quaternion(const RotationSequence rs, const vec3& angles)
+        {
+            switch (rs)
+            {
+                case XYZ:
+                    *this = Quaternion(vec3(0, 1, 0), angles[0]) *
+                            Quaternion(vec3(0, 1, 0), angles[1]) *
+                            Quaternion(vec3(0, 0, 1), angles[2]);
+                    break;
+
+                default:
+                    throw std::runtime_error("Unknown rotation sequence");
+            }
+        }
+
+        inline
+        scalar real() const
+        {
+            return p_real;
+        }
+
+        inline
+        scalar& real()
+        {
+            return p_real;
+        }
+
+        inline
+        vec3 imag() const
+        {
+            return p_imag;
+        }
+
+        inline
+        vec3& imag()
+        {
+            return p_imag;
+        }
+
+        inline
+        scalar& operator[](Size n)
+        {
+            if (n > 3)
+                throw hpr::OutOfRange();
+            if (n == 0)
+                return p_real;
+            else
+                return p_imag[n - 1];
+        }
+
+        inline
+        scalar operator[](Size n) const
+        {
+            if (n > 3)
+                throw hpr::OutOfRange();
+            if (n == 0)
+                return p_real;
+            else
+                return p_imag[n - 1];
+        }
+
+        inline
+        void operator+=(const Quaternion& q)
+        {
+            p_real += q.p_real;
+            p_imag += q.p_imag;
+        }
+
+        inline
+        void operator-=(const Quaternion& q)
+        {
+            p_real -= q.p_real;
+            p_imag -= q.p_imag;
+        }
+
+        inline
+        void operator*=(const Quaternion& q)
+        {
+            scalar temp = p_real;
+            p_real = p_real * q.p_real - dot(p_imag, q.p_imag);
+            p_imag = temp * q.p_imag + q.p_real * p_imag + cross(p_imag, q.p_imag);
+        }
+
+        inline
+        void operator/=(const Quaternion& q)
+        {
+            operator*=(inverse(q));
+        }
+
+        inline
+        void operator*=(const scalar s)
+        {
+            p_real *= s;
+            p_imag *= s;
+        }
+
+        inline
+        void operator/=(const scalar s)
+        {
+            p_real /= s;
+            p_imag /= s;
+        }
+    };
+
+    inline
+    bool equal(const Quaternion& lhs, const Quaternion& rhs)
+    {
+        return lhs.real() == rhs.real() && lhs.imag() == rhs.imag();
+    }
+
+    inline
+    bool operator==(const Quaternion& lhs, const Quaternion& rhs)
+    {
+        return equal(lhs, rhs);
+    }
+
+    inline
+    bool operator!=(const Quaternion& lhs, const Quaternion& rhs)
+    {
+        return !equal(lhs, rhs);
+    }
+
+    inline
+    Quaternion operator+(const Quaternion& lhs, const Quaternion& rhs)
+    {
+        return {lhs.real() + rhs.real(), lhs.imag() + rhs.imag()};
+    }
+
+    inline
+    Quaternion operator-(const Quaternion& q)
+    {
+        return {q.real(), q.imag()};
+    }
+
+    inline
+    Quaternion operator-(const Quaternion& lhs, const Quaternion& rhs)
+    {
+        return {lhs.real() - rhs.real(), lhs.imag() - rhs.imag()};
+    }
+
+    inline
+    Quaternion operator*(const Quaternion& lhs, const Quaternion& rhs)
+    {
+        return {
+            lhs.real() * rhs.real() - dot(lhs.imag(), rhs.imag()),
+            lhs.real() * rhs.imag() + rhs.real() * lhs.imag() + cross(lhs.imag(), rhs.imag())
+        };
+    }
+
+    inline
+    Quaternion operator/(const Quaternion& lhs, const Quaternion& rhs)
+    {
+        return lhs * inverse(rhs);
+    }
+
+    inline
+    Quaternion operator*(const scalar s, const Quaternion& q)
+    {
+        return {s * q.real(), s * q.imag()};
+    }
+
+    inline
+    Quaternion operator*(const Quaternion& q, const scalar s)
+    {
+        return {q.real() * s, q.imag() * s};
+    }
+
+    inline
+    Quaternion operator/(const Quaternion& q, const scalar s)
+    {
+        return {q.real() / s, q.imag() / s};
+    }
+
+    inline
+    scalar norm(const Quaternion& q)
+    {
+        return sqrt(pow(q.real(), 2) + dot(q.imag(), q.imag()));
+    }
+
+    inline
+    Quaternion conjugate(const Quaternion& q)
+    {
+        return {q.real(), -q.imag()};
+    }
+
+    inline
+    Quaternion inverse(const Quaternion& q)
+    {
+        return conjugate(q) / pow(norm(q), 2);
+    }
+
+    inline
+    Quaternion normalize(const Quaternion& q)
+    {
+        return q / norm(q);
+    }
+
+    inline
+    vec3 rotate(const vec3& point, const vec3& axis, const scalar& angle)
+    {
+        Quaternion p {point};
+        Quaternion q {normalize(axis), angle};
+        return (q * p * inverse(q)).imag();
+    }
+
+    void decompose(const Quaternion& q, vec3& axis, scalar& angle)
+    {
+        const scalar qnorm = norm(q.imag());
+        axis = q.imag() / qnorm;
+        angle = 2 * atan2(qnorm, q.real());
+    }
+
+} // end namespace hpr
\ No newline at end of file
diff --git a/source/hpr/math/vector.hpp b/source/hpr/math/vector.hpp
new file mode 100644
index 0000000..c2cf8d1
--- /dev/null
+++ b/source/hpr/math/vector.hpp
@@ -0,0 +1,321 @@
+#pragma once
+
+#include <hpr/numeric.hpp>
+#include <hpr/container/static_array.hpp>
+
+
+namespace hpr
+{
+    // forward declarations
+
+    template <IsReal T, Size S> requires (S >= 0)
+    class Vector;
+
+    template <IsReal T, Size S>
+    using SubVector = typename std::conditional<S >= 2, Vector<T, S - 1>, Vector<T, 1>>::type;
+
+    // type traits
+
+    template <typename T>
+    struct is_vector : public std::false_type {};
+
+    template <typename T, Size S>
+    struct is_vector<Vector<T, S>> : public std::true_type {};
+
+    // concepts
+
+    template <typename T>
+    concept IsVector = is_vector<T>::value;
+
+    // aliases
+
+    template <typename Type, Size S>
+    using vec = Vector<Type, S>;
+
+    using vec2 = Vector<scalar, 2>;
+    using vec3 = Vector<scalar, 3>;
+    using vec4 = Vector<scalar, 4>;
+
+
+    template <IsReal Type, Size S> requires (S >= 0)
+    class Vector : public StaticArray<Type, S>
+    {
+
+        using base = StaticArray<Type, S>;
+
+    public:
+
+        using value_type = Type;
+        using size_type = Size;
+        using pointer = Type*;
+        using reference = Type&;
+        using iterator = Iterator<Type>;
+        using const_iterator = Iterator<const Type>;
+
+    public:
+
+        //! null constructor
+        constexpr
+        Vector() :
+            base {}
+        {}
+
+        //! copy constructor
+        constexpr
+        Vector(const Vector& vs) :
+            base {static_cast<base>(vs)}
+        {}
+
+        //! move constructor
+        constexpr
+        Vector(Vector&& vs) noexcept :
+            base {std::forward<base>(static_cast<base>(vs))}
+        {}
+
+        //! copy assignment operator
+        constexpr
+        Vector& operator=(const Vector& vs)
+        {
+            base::operator=(vs);
+            return *this;
+        }
+
+        //! move assignment operator
+        constexpr
+        Vector& operator=(Vector&& vs) noexcept
+        {
+            swap(*this, vs);
+            return *this;
+        }
+
+        //! destructor
+        virtual
+        ~Vector() = default;
+
+        //! copy constructor from base
+        constexpr
+        Vector(const base& arr) :
+            base {arr}
+        {}
+
+        //! move constructor from base
+        constexpr
+        Vector(base&& arr) :
+            base {std::forward<base>(arr)}
+        {}
+
+        //! construct from iterators
+        constexpr
+        Vector(typename base::iterator start, typename base::iterator end) :
+            base {start, end}
+        {}
+
+        //! construct from constant iterators
+        constexpr
+        Vector(typename base::const_iterator start, typename base::const_iterator end) :
+            base {start, end}
+        {}
+
+        //! construct from initializer list
+        constexpr
+        Vector(std::initializer_list<value_type> list) :
+            base {list}
+        {}
+
+        //! copy constructor with variadic args
+        template <IsReal... Args>
+        constexpr
+        Vector(const value_type& v, const Args& ...args) requires (S == 1 + sizeof...(args)):
+            base {v, static_cast<value_type>(args)...}
+        {}
+
+        //! move constructor with variadic args
+        template <IsReal... Args>
+        constexpr
+        Vector(value_type&& v, Args&& ...args) requires (S == 1 + sizeof...(args)):
+            base {v, static_cast<value_type>(std::forward<Args>(args))...}
+        {}
+
+        //! copy constructor with sub vector and value
+        constexpr
+        Vector(const SubVector<Type, S>& svs, const value_type& v) requires (S >= 1):
+            base {}
+        {
+            for (auto n = 0; n < svs.size(); ++n)
+                (*this)[n] = svs[n];
+            (*this)[svs.size()] = v;
+        }
+
+        //! copy constructor from greater vector
+        template <Size GS> requires (GS > S)
+        constexpr explicit
+        Vector(const Vector<Type, GS>& vs) :
+            base {vs.begin(), vs.begin() + S}
+        {}
+
+    };
+
+    // global operators
+
+    template <IsReal T, Size S> inline Vector<T, S> operator+(const Vector<T, S>& lhs) { Vector<T, S> vs; for (Size n = 0; n < S; ++n) vs[n] = lhs[n]; return vs; }
+    template <IsReal T, Size S> inline Vector<T, S> operator-(const Vector<T, S>& lhs) { Vector<T, S> vs; for (Size n = 0; n < S; ++n) vs[n] = -lhs[n]; return vs; }
+
+    template <IsReal T, Size S> inline Vector<T, S>& operator+=(Vector<T, S>& lhs, const Vector<T, S>& rhs) { for (Size n = 0; n < S; ++n) lhs[n] += rhs[n]; return lhs; }
+    template <IsReal T, Size S> inline Vector<T, S>& operator-=(Vector<T, S>& lhs, const Vector<T, S>& rhs) { for (Size n = 0; n < S; ++n) lhs[n] -= rhs[n]; return lhs; }
+    template <IsReal T, Size S> inline Vector<T, S>& operator*=(Vector<T, S>& lhs, const Vector<T, S>& rhs) { for (Size n = 0; n < S; ++n) lhs[n] *= rhs[n]; return lhs; }
+    template <IsReal T, Size S> inline Vector<T, S>& operator/=(Vector<T, S>& lhs, const Vector<T, S>& rhs) { for (Size n = 0; n < S; ++n) lhs[n] /= rhs[n]; return lhs; }
+
+    template <IsReal T, Size S> inline Vector<T, S> operator+(const Vector<T, S>& lhs, const Vector<T, S>& rhs) { Vector<T, S> vs {lhs}; for (Size n = 0; n < S; ++n) vs[n] += rhs[n]; return vs; }
+    template <IsReal T, Size S> inline Vector<T, S> operator-(const Vector<T, S>& lhs, const Vector<T, S>& rhs) { Vector<T, S> vs {lhs}; for (Size n = 0; n < S; ++n) vs[n] -= rhs[n]; return vs; }
+    template <IsReal T, Size S> inline Vector<T, S> operator*(const Vector<T, S>& lhs, const Vector<T, S>& rhs) { Vector<T, S> vs {lhs}; for (Size n = 0; n < S; ++n) vs[n] *= rhs[n]; return vs; }
+    template <IsReal T, Size S> inline Vector<T, S> operator/(const Vector<T, S>& lhs, const Vector<T, S>& rhs) { Vector<T, S> vs {lhs}; for (Size n = 0; n < S; ++n) vs[n] /= rhs[n]; return vs; }
+
+    template <IsReal T, Size S> inline bool operator==(const Vector<T, S>& lhs, const Vector<T, S>& rhs) { for (Size n = 0; n < S; ++n) if (lhs[n] != rhs[n]) return false; return true; }
+    template <IsReal T, Size S> inline bool operator!=(const Vector<T, S>& lhs, const Vector<T, S>& rhs) { for (Size n = 0; n < S; ++n) if (lhs[n] == rhs[n]) return false; return true; }
+
+
+    template <IsReal T, Size S> inline Vector<T, S>& operator+=(Vector<T, S>& lhs, const T& rhs) { for (Size n = 0; n < S; ++n) lhs[n] += rhs; return lhs; }
+    template <IsReal T, Size S> inline Vector<T, S>& operator-=(Vector<T, S>& lhs, const T& rhs) { for (Size n = 0; n < S; ++n) lhs[n] -= rhs; return lhs; }
+    template <IsReal T, Size S> inline Vector<T, S>& operator*=(Vector<T, S>& lhs, const T& rhs) { for (Size n = 0; n < S; ++n) lhs[n] *= rhs; return lhs; }
+    template <IsReal T, Size S> inline Vector<T, S>& operator/=(Vector<T, S>& lhs, const T& rhs) { for (Size n = 0; n < S; ++n) lhs[n] /= rhs; return lhs; }
+
+    template <IsReal T, Size S> inline Vector<T, S> operator+(const Vector<T, S>& lhs, const T& rhs) { Vector<T, S> vs {lhs}; for (Size n = 0; n < S; ++n) vs[n] += rhs; return vs; }
+    template <IsReal T, Size S> inline Vector<T, S> operator-(const Vector<T, S>& lhs, const T& rhs) { Vector<T, S> vs {lhs}; for (Size n = 0; n < S; ++n) vs[n] -= rhs; return vs; }
+    template <IsReal T, IsReal T2, Size S> inline Vector<T, S> operator*(const Vector<T, S>& lhs, const T2& rhs) { Vector<T, S> vs {lhs}; for (Size n = 0; n < S; ++n) vs[n] *= static_cast<T>(rhs); return vs; }
+    template <IsReal T, IsReal T2, Size S> inline Vector<T, S> operator/(const Vector<T, S>& lhs, const T2& rhs) { Vector<T, S> vs {lhs}; for (Size n = 0; n < S; ++n) vs[n] /= static_cast<T>(rhs); return vs; }
+
+    template <IsReal T, Size S> inline Vector<T, S> operator+(const T& lhs, const Vector<T, S>& rhs) { Vector<T, S> vs {rhs}; for (Size n = 0; n < S; ++n) vs[n] += lhs; return vs; }
+    template <IsReal T, Size S> inline Vector<T, S> operator-(const T& lhs, const Vector<T, S>& rhs) { Vector<T, S> vs {rhs}; for (Size n = 0; n < S; ++n) vs[n] -= lhs; return vs; }
+    template <IsReal T, Size S> inline Vector<T, S> operator*(const T& lhs, const Vector<T, S>& rhs) { Vector<T, S> vs {rhs}; for (Size n = 0; n < S; ++n) vs[n] *= lhs; return vs; }
+    template <IsReal T, Size S> inline Vector<T, S> operator/(const T& lhs, const Vector<T, S>& rhs) { Vector<T, S> vs {rhs}; for (Size n = 0; n < S; ++n) vs[n] /= lhs; return vs; }
+
+    // boolean operations
+
+    template <typename Type, Size S>
+    inline
+    Vector<bool, S> equal(const Vector<Type, S>& lhs, const Vector<Type, S>& rhs, const scalar& precision = hpr::precision())
+    {
+        Vector<bool, S> vs;
+        for (auto n = 0; n < S; ++n)
+            vs[n] = equal(lhs[n], rhs[n], precision);
+        return vs;
+    }
+
+    template <Size S>
+    inline
+    bool any(const Vector<bool, S>& vs)
+    {
+        bool res = false;
+        for (auto e : vs)
+            res = res || e;
+        return res;
+    }
+
+    template <Size S>
+    inline
+    bool all(const Vector<bool, S>& vs)
+    {
+        bool res = true;
+        for (auto e : vs)
+            res = res && e;
+        return res;
+    }
+
+    // per element operations
+
+    template <typename Type, Size S>
+    inline
+    Vector<Type, S> abs(const Vector<Type, S>& vs)
+    {
+        Vector<Type, S> res;
+        for (auto n = 0; n < S; ++n)
+            res[n] = abs(vs[n]);
+        return res;
+    }
+
+    template <typename Type, Size S>
+    inline
+    Type sum(const Vector<Type, S>& vs)
+    {
+        Type sum {};
+        for (const Type& v : vs)
+            sum += v;
+        return sum;
+    }
+
+    template <typename Type, Size S>
+    inline
+    Vector<Type, S> pow(const Vector<Type, S>& vs, scalar degree)
+    {
+        Vector<Type, S> res;
+        for (auto n = 0; n < S; ++n)
+            res[n] = pow(vs[n], degree);
+        return res;
+    }
+
+    // vector operations
+
+    template <typename Type, Size S>
+    inline
+    Type norm(const Vector<Type, S>& vs)
+    {
+        return sqrt(sum(pow(abs(vs), 2)));
+    }
+
+    template <typename Type, Size S>
+    inline
+    Type dot(const Vector<Type, S>& lhs, const Vector<Type, S>& rhs)
+    {
+        return sum(lhs * rhs);
+    }
+
+    template <typename Type, Size S>
+    inline
+    Type length(const Vector<Type, S>& vs)
+    {
+        return sqrt(dot(vs, vs));
+    }
+
+    template <typename Type, Size S>
+    inline
+    Type mag(const Vector<Type, S>& vs)
+    {
+        return length(vs);
+    }
+
+    template <typename Type, Size S>
+    inline
+    Type distance(const Vector<Type, S>& vs1, const Vector<Type, S>& vs2)
+    {
+        return length(vs1 - vs2);
+    }
+
+    template <typename Type, Size S>
+    inline
+    Vector<Type, S> normalize(const Vector<Type, S>& vs)
+    {
+        return vs * isqrt(dot(vs, vs));
+    }
+
+    template <typename Type>
+    inline
+    Type angle(const Vector<Type, 3>& lhs, const Vector<Type, 3>& rhs)
+    {
+        scalar cos = dot(lhs, rhs) / (norm(lhs) * norm(rhs));
+        return acos(cos); //clip(cos, -1., 1.));
+    }
+
+    // vector 3 operations
+
+    template <typename Type>
+    inline
+    Vector<Type, 3> cross(const Vector<Type, 3>& lhs, const Vector<Type, 3>& rhs)
+    {
+        return Vector<Type, 3>(
+            lhs[1] * rhs[2] - lhs[2] * rhs[1],
+            lhs[2] * rhs[0] - lhs[0] * rhs[2],
+            lhs[0] * rhs[1] - lhs[1] * rhs[0]
+        );
+    }
+
+} // end namespace hpr
\ No newline at end of file
diff --git a/source/hpr/numeric.hpp b/source/hpr/numeric.hpp
new file mode 100644
index 0000000..84c0c21
--- /dev/null
+++ b/source/hpr/numeric.hpp
@@ -0,0 +1,305 @@
+#pragma once
+
+#include <type_traits>
+#include <cmath>
+#include <limits>
+#include <numbers>
+
+#ifdef __MSC_VER
+#define HPR_CONSTEXPR
+#else
+#define HPR_CONSTEXPR constexpr
+#endif
+
+
+namespace hpr
+{
+    using Size = std::size_t;
+
+    template <typename Type>
+    struct is_size : public std::integral_constant<bool, std::is_integral<Type>::value && std::is_unsigned<Type>::value> {};
+
+    template <typename Type>
+    concept IsSize = is_size<Type>::value || std::convertible_to<Type, Size>;
+
+
+    template <typename Type>
+    struct is_integer : public std::is_integral<Type> {};
+
+    template <typename Type>
+    concept IsInteger = is_integer<Type>::value;
+
+
+    template <typename Type>
+    struct is_scalar : public std::is_floating_point<Type> {};
+
+    template <typename Type>
+    concept IsScalar = is_scalar<Type>::value;
+
+    template <typename Type>
+    concept IsReal = is_integer<Type>::value || is_scalar<Type>::value;
+
+
+#ifndef HPR_SCALAR
+#define HPR_SCALAR float
+#endif
+
+#ifdef HPR_SCALAR_IMPLEMENTATION
+    template <IsScalar T>
+    class Scalar
+    {
+
+    public:
+
+        using type = Scalar<T>;
+        using value_type = T;
+
+    protected:
+
+        value_type p_value;
+
+    public:
+
+        // constructors
+
+        constexpr Scalar() : p_value {} {}
+
+        template <IsScalar X> constexpr Scalar(const Scalar<X>& value) : p_value {static_cast<value_type>(value.p_value)} {}
+
+        template <IsReal X> constexpr Scalar(const X& value) : p_value {static_cast<value_type>(value)} {}
+
+        template <IsScalar X> constexpr type& operator=(const Scalar<X>& value) { p_value = static_cast<value_type>(value.p_value); return *this; }
+
+        template <IsReal X> constexpr type& operator=(const X& value) { p_value = static_cast<value_type>(value); return *this; }
+
+        virtual constexpr ~Scalar() = default;
+
+        // conversion
+
+        constexpr operator double() const { return static_cast<double>(p_value); }
+
+        constexpr operator float() const { return static_cast<float>(p_value); }
+
+        constexpr operator long double() const { return static_cast<long double>(p_value); }
+
+        constexpr operator bool() const { return static_cast<bool>(p_value); }
+
+        // access
+
+        [[nodiscard]] constexpr value_type value() const { return p_value; }
+
+        constexpr value_type& value() { return p_value; }
+    };
+
+
+    using scalar_type = HPR_SCALAR;
+    using scalar = Scalar<HPR_SCALAR>;
+
+
+    // template<> scalar::value_type scalar::s_precision = static_cast<scalar::value_type>(1e-15);
+
+    /// scalar vs scalar
+
+    constexpr scalar operator+(const scalar& s) { return s; }
+    constexpr scalar operator-(const scalar& s) { return -s.value(); }
+    constexpr bool operator!(const scalar& s) { return !static_cast<bool>(s.value()); }
+
+    constexpr scalar& operator+=(scalar& lhs, const scalar& rhs) { lhs.value() += static_cast<scalar::value_type>(rhs.value()); return lhs; }
+    constexpr scalar& operator-=(scalar& lhs, const scalar& rhs) { lhs.value() -= static_cast<scalar::value_type>(rhs.value()); return lhs; }
+    constexpr scalar& operator*=(scalar& lhs, const scalar& rhs) { lhs.value() *= static_cast<scalar::value_type>(rhs.value()); return lhs; }
+    constexpr scalar& operator/=(scalar& lhs, const scalar& rhs) { lhs.value() /= static_cast<scalar::value_type>(rhs.value()); return lhs; }
+
+    constexpr scalar operator+(const scalar& lhs, const scalar& rhs) { return lhs.value() + rhs.value(); }
+    constexpr scalar operator-(const scalar& lhs, const scalar& rhs) { return lhs.value() - rhs.value(); }
+    constexpr scalar operator*(const scalar& lhs, const scalar& rhs) { return lhs.value() * rhs.value(); }
+    constexpr scalar operator/(const scalar& lhs, const scalar& rhs) { return lhs.value() / rhs.value(); }
+
+    constexpr bool operator==(const scalar& lhs, const scalar& rhs) { return lhs.value() == rhs.value(); }
+    constexpr bool operator!=(const scalar& lhs, const scalar& rhs) { return lhs.value() != rhs.value(); }
+    constexpr bool operator&&(const scalar& lhs, const scalar& rhs) { return static_cast<bool>(lhs) && static_cast<bool>(rhs); }
+    constexpr bool operator||(const scalar& lhs, const scalar& rhs) { return static_cast<bool>(lhs) || static_cast<bool>(rhs); }
+    constexpr bool operator>(const scalar& lhs, const scalar& rhs) { return lhs.value() > rhs.value(); }
+    constexpr bool operator<(const scalar& lhs, const scalar& rhs) { return lhs.value() < rhs.value(); }
+    constexpr bool operator>=(const scalar& lhs, const scalar& rhs) { return lhs.value() >= rhs.value(); }
+    constexpr bool operator<=(const scalar& lhs, const scalar& rhs) { return lhs.value() <= rhs.value(); }
+
+    // std::ostream& operator<<(std::ostream& stream, const scalar& s) { return stream << s.value(); }
+
+    /// scalar vs Scalar<T>
+
+    template <IsScalar T> constexpr scalar& operator+=(scalar& lhs, const Scalar<T>& rhs) { lhs.value() += static_cast<scalar::value_type>(rhs.value()); return lhs; }
+    template <IsScalar T> constexpr scalar& operator-=(scalar& lhs, const Scalar<T>& rhs) { lhs.value() -= static_cast<scalar::value_type>(rhs.value()); return lhs; }
+    template <IsScalar T> constexpr scalar& operator*=(scalar& lhs, const Scalar<T>& rhs) { lhs.value() *= static_cast<scalar::value_type>(rhs.value()); return lhs; }
+    template <IsScalar T> constexpr scalar& operator/=(scalar& lhs, const Scalar<T>& rhs) { lhs.value() /= static_cast<scalar::value_type>(rhs.value()); return lhs; }
+
+    template <IsScalar T> constexpr scalar operator+(const scalar& lhs, const Scalar<T>& rhs) { return lhs.value() + static_cast<scalar::value_type>(rhs.value()); }
+    template <IsScalar T> constexpr scalar operator+(const Scalar<T>& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs.value()) + rhs.value(); }
+    template <IsScalar T> constexpr scalar operator-(const scalar& lhs, const Scalar<T>& rhs) { return lhs.value() - static_cast<scalar::value_type>(rhs.value()); }
+    template <IsScalar T> constexpr scalar operator-(const Scalar<T>& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs.value()) - rhs.value(); }
+    template <IsScalar T> constexpr scalar operator*(const scalar& lhs, const Scalar<T>& rhs) { return lhs.value() * static_cast<scalar::value_type>(rhs.value()); }
+    template <IsScalar T> constexpr scalar operator*(const Scalar<T>& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs.value()) * rhs.value(); }
+    template <IsScalar T> constexpr scalar operator/(const scalar& lhs, const Scalar<T>& rhs) { return lhs.value() / static_cast<scalar::value_type>(rhs.value()); }
+    template <IsScalar T> constexpr scalar operator/(const Scalar<T>& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs.value()) / rhs.value(); }
+
+    template <IsScalar T> constexpr bool operator==(const scalar& lhs, const Scalar<T>& rhs) { return lhs.value() == static_cast<scalar::value_type>(rhs.value()); }
+    template <IsScalar T> constexpr bool operator==(const Scalar<T>& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs.value()) == rhs.value(); }
+    template <IsScalar T> constexpr bool operator!=(const scalar& lhs, const Scalar<T>& rhs) { return lhs.value() != static_cast<scalar::value_type>(rhs.value()); }
+    template <IsScalar T> constexpr bool operator!=(const Scalar<T>& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs.value()) != rhs.value(); }
+    template <IsScalar T> constexpr bool operator&&(const scalar& lhs, const Scalar<T>& rhs) { return static_cast<bool>(lhs) && static_cast<bool>(rhs); }
+    template <IsScalar T> constexpr bool operator&&(const Scalar<T>& lhs, const scalar& rhs) { return static_cast<bool>(lhs) && static_cast<bool>(rhs); }
+    template <IsScalar T> constexpr bool operator||(const scalar& lhs, const Scalar<T>& rhs) { return static_cast<bool>(lhs) || static_cast<bool>(rhs); }
+    template <IsScalar T> constexpr bool operator||(const Scalar<T>& lhs, const scalar& rhs) { return static_cast<bool>(lhs) || static_cast<bool>(rhs); }
+    template <IsScalar T> constexpr bool operator>(const scalar& lhs, const Scalar<T>& rhs) { return lhs.value() > static_cast<scalar::value_type>(rhs.value()); }
+    template <IsScalar T> constexpr bool operator>(const Scalar<T>& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs.value()) > rhs.value(); }
+    template <IsScalar T> constexpr bool operator<(const scalar& lhs, const Scalar<T>& rhs) { return lhs.value() < static_cast<scalar::value_type>(rhs.value()); }
+    template <IsScalar T> constexpr bool operator<(const Scalar<T>& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs.value()) < rhs.value(); }
+    template <IsScalar T> constexpr bool operator>=(const scalar& lhs, const Scalar<T>& rhs) { return lhs.value() >= static_cast<scalar::value_type>(rhs.value()); }
+    template <IsScalar T> constexpr bool operator>=(const Scalar<T>& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs.value()) >= rhs.value(); }
+    template <IsScalar T> constexpr bool operator<=(const scalar& lhs, const Scalar<T>& rhs) { return lhs.value() <= static_cast<scalar::value_type>(rhs.value()); }
+    template <IsScalar T> constexpr bool operator<=(const Scalar<T>& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs.value()) <= rhs.value(); }
+
+    template <IsScalar T> std::ostream& operator<<(std::ostream& stream, const Scalar<T>& s) { return stream << s.value(); }
+    template <IsScalar T> std::istream& operator>>(std::istream& stream, Scalar<T>& s) { return stream >> s.value(); }
+
+    /// scalar vs real
+
+    template <IsReal T> constexpr scalar& operator+=(scalar& lhs, const T& rhs) { lhs.value() += static_cast<scalar::value_type>(rhs); return lhs; }
+    template <IsReal T> constexpr scalar& operator-=(scalar& lhs, const T& rhs) { lhs.value() -= static_cast<scalar::value_type>(rhs); return lhs; }
+    template <IsReal T> constexpr scalar& operator*=(scalar& lhs, const T& rhs) { lhs.value() *= static_cast<scalar::value_type>(rhs); return lhs; }
+    template <IsReal T> constexpr scalar& operator/=(scalar& lhs, const T& rhs) { lhs.value() /= static_cast<scalar::value_type>(rhs); return lhs; }
+
+    template <IsReal T> constexpr T& operator+=(T& lhs, const scalar& rhs) { lhs += static_cast<T>(rhs); return lhs; }
+    template <IsReal T> constexpr T& operator-=(T& lhs, const scalar& rhs) { lhs -= static_cast<T>(rhs); return lhs; }
+    template <IsReal T> constexpr T& operator*=(T& lhs, const scalar& rhs) { lhs *= static_cast<T>(rhs); return lhs; }
+    template <IsReal T> constexpr T& operator/=(T& lhs, const scalar& rhs) { lhs /= static_cast<T>(rhs); return lhs; }
+
+    template <IsReal T> constexpr scalar operator+(const scalar& lhs, const T& rhs) { return lhs.value() + static_cast<scalar::value_type>(rhs); }
+    template <IsReal T> constexpr scalar operator+(const T& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs) + rhs.value(); }
+    template <IsReal T> constexpr scalar operator-(const scalar& lhs, const T& rhs) { return lhs.value() - static_cast<scalar::value_type>(rhs); }
+    template <IsReal T> constexpr scalar operator-(const T& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs) - rhs.value(); }
+    template <IsReal T> constexpr scalar operator*(const scalar& lhs, const T& rhs) { return lhs.value() * static_cast<scalar::value_type>(rhs); }
+    template <IsReal T> constexpr scalar operator*(const T& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs) * rhs.value(); }
+    template <IsReal T> constexpr scalar operator/(const scalar& lhs, const T& rhs) { return lhs.value() / static_cast<scalar::value_type>(rhs); }
+    template <IsReal T> constexpr scalar operator/(const T& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs) / rhs.value(); }
+
+    template <IsReal T> constexpr bool operator==(const scalar& lhs, const T& rhs) { return lhs.value() == static_cast<scalar::value_type>(rhs); }
+    template <IsReal T> constexpr bool operator==(const T& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs) == rhs.value(); }
+    template <IsReal T> constexpr bool operator!=(const scalar& lhs, const T& rhs) { return lhs.value() != static_cast<scalar::value_type>(rhs); }
+    template <IsReal T> constexpr bool operator!=(const T& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs) != rhs.value(); }
+    template <IsReal T> constexpr bool operator&&(const scalar& lhs, const T& rhs) { return static_cast<bool>(lhs) && static_cast<bool>(rhs); }
+    template <IsReal T> constexpr bool operator&&(const T& lhs, const scalar& rhs) { return static_cast<bool>(lhs) && static_cast<bool>(rhs); }
+    template <IsReal T> constexpr bool operator||(const scalar& lhs, const T& rhs) { return static_cast<bool>(lhs) || static_cast<bool>(rhs); }
+    template <IsReal T> constexpr bool operator||(const T& lhs, const scalar& rhs) { return static_cast<bool>(lhs) || static_cast<bool>(rhs); }
+    template <IsReal T> constexpr bool operator>(const scalar& lhs, const T& rhs) { return lhs.value() > static_cast<scalar::value_type>(rhs); }
+    template <IsReal T> constexpr bool operator>(const T& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs) > rhs.value(); }
+    template <IsReal T> constexpr bool operator<(const scalar& lhs, const T& rhs) { return lhs.value() < static_cast<scalar::value_type>(rhs); }
+    template <IsReal T> constexpr bool operator<(const T& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs) < rhs.value(); }
+    template <IsReal T> constexpr bool operator>=(const scalar& lhs, const T& rhs) { return lhs.value() >= static_cast<scalar::value_type>(rhs); }
+    template <IsReal T> constexpr bool operator>=(const T& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs) >= rhs.value(); }
+    template <IsReal T> constexpr bool operator<=(const scalar& lhs, const T& rhs) { return lhs.value() <= static_cast<scalar::value_type>(rhs); }
+    template <IsReal T> constexpr bool operator<=(const T& lhs, const scalar& rhs) { return static_cast<scalar::value_type>(lhs) <= rhs.value(); }
+
+#else
+
+    using scalar_type = HPR_SCALAR;
+    using scalar = HPR_SCALAR;
+
+#endif
+
+
+    // constants
+
+    HPR_CONSTEXPR scalar pi() { return std::numbers::pi_v<scalar_type>; }
+
+    HPR_CONSTEXPR scalar e() { return std::numbers::e_v<scalar_type>; }
+
+    // transcendentals
+
+    template <IsReal T> HPR_CONSTEXPR scalar cos(const T& s) { return std::cos(static_cast<scalar_type>(s));}
+
+    template <IsReal T> HPR_CONSTEXPR scalar acos(const T& s) { return std::acos(static_cast<scalar_type>(s)); }
+
+    template <IsReal T> HPR_CONSTEXPR scalar cosh(const T& s) { return std::cosh(static_cast<scalar_type>(s)); }
+
+    template <IsReal T> HPR_CONSTEXPR scalar acosh(const T& s) { return std::acosh(static_cast<scalar_type>(s)); }
+
+    template <IsReal T> HPR_CONSTEXPR scalar sin(const T& s) { return std::sin(static_cast<scalar_type>(s)); }
+
+    template <IsReal T> HPR_CONSTEXPR scalar asin(const T& s) { return std::asin(static_cast<scalar_type>(s)); }
+
+    template <IsReal T> HPR_CONSTEXPR scalar sinh(const T& s) { return std::sinh(static_cast<scalar_type>(s)); }
+
+    template <IsReal T> HPR_CONSTEXPR scalar asinh(const T& s) { return std::asinh(static_cast<scalar_type>(s)); }
+
+    template <IsReal T> HPR_CONSTEXPR scalar tan(const T& s) { return std::tan(static_cast<scalar_type>(s)); }
+
+    template <IsReal T> HPR_CONSTEXPR scalar atan(const T& s) { return std::atan(static_cast<scalar_type>(s)); }
+
+    template <IsReal T, IsReal X> HPR_CONSTEXPR scalar atan2(const T& s, const X& s2) { return std::atan2(static_cast<scalar_type>(s), static_cast<scalar_type>(s2)); }
+
+    template <IsReal T> HPR_CONSTEXPR scalar tanh(const T& s) { return std::tanh(static_cast<scalar_type>(s)); }
+
+    template <IsReal T> HPR_CONSTEXPR scalar atanh(const T& s) { return std::atanh(static_cast<scalar_type>(s)); }
+
+    template <IsReal T> HPR_CONSTEXPR scalar exp(const T& s) { return std::exp(static_cast<scalar_type>(s)); }
+
+    template <IsReal T> HPR_CONSTEXPR scalar log(const T& s) { return std::log(static_cast<scalar_type>(s)); }
+
+    template <IsReal T> HPR_CONSTEXPR scalar log10(const T& s) { return std::log10(static_cast<scalar_type>(s)); }
+
+    template <IsReal T, IsReal X> HPR_CONSTEXPR scalar pow(const T& s, const X& d) { return std::pow(static_cast<scalar_type>(s), static_cast<scalar_type>(d)); }
+
+    template <IsReal T> HPR_CONSTEXPR scalar sqrt(const T& s) { return std::sqrt(static_cast<scalar_type>(s)); }
+
+    template <IsReal T> HPR_CONSTEXPR scalar isqrt(const T& s) { return static_cast<T>(1) / sqrt(static_cast<scalar_type>(s)); }
+
+    //
+
+    static scalar_type scalar_precision = static_cast<scalar_type>(1e-15);
+
+    inline scalar_type precision()
+    {
+        return scalar_precision;
+    }
+
+    template <IsReal T> inline void precision(const T& precision)
+    {
+        scalar_precision = static_cast<scalar_type>(precision);
+    }
+
+
+    HPR_CONSTEXPR scalar inf() { return std::numeric_limits<scalar_type>::infinity(); }
+
+    HPR_CONSTEXPR scalar epsilon() { return std::numeric_limits<scalar_type>::epsilon(); }
+
+    HPR_CONSTEXPR scalar isnan(const scalar& s) { return std::isnan(static_cast<scalar_type>(s)); }
+
+    HPR_CONSTEXPR scalar abs(const scalar& s) { return std::abs(static_cast<scalar_type>(s)); }
+
+    HPR_CONSTEXPR scalar mag(const scalar& s) { return std::abs(static_cast<scalar_type>(s)); }
+
+    HPR_CONSTEXPR bool equal(const scalar& lhs, const scalar& rhs, const scalar& precision = hpr::precision()) { return abs(lhs - rhs) < precision; }
+
+    //! Convert degrees to radians
+    HPR_CONSTEXPR scalar rad(const scalar& s) { return s * pi() / static_cast<scalar_type>(180); }
+
+    //! Convert radians to degrees
+    HPR_CONSTEXPR scalar deg(const scalar& s) { return s * static_cast<scalar_type>(180) / pi(); }
+
+    HPR_CONSTEXPR scalar min(const scalar& s1, const scalar& s2) { return std::min(static_cast<scalar_type>(s1), static_cast<scalar_type>(s2));}
+
+    HPR_CONSTEXPR scalar max(const scalar& s1, const scalar& s2) { return std::max(static_cast<scalar_type>(s1), static_cast<scalar_type>(s2)); }
+
+    HPR_CONSTEXPR scalar clip(const scalar& s, const scalar& sMin, const scalar& sMax) { return min(sMax, max(s, sMin)); }
+
+
+} // end namespace hpr
+
+#ifdef HPR_SCALAR_IMPLEMENTATION
+
+namespace std
+{
+    template <class T> struct is_floating_point<hpr::scalar> : true_type {};
+} // end namespace std
+
+#endif
\ No newline at end of file
diff --git a/source/hpr/tests/test_container.cpp b/source/hpr/tests/test_container.cpp
new file mode 100644
index 0000000..092791d
--- /dev/null
+++ b/source/hpr/tests/test_container.cpp
@@ -0,0 +1,100 @@
+#include <gtest/gtest.h>
+
+
+#include <hpr/container/static_array.hpp>
+
+TEST(containers, StaticArray)
+{
+    hpr::StaticArray<float, 3> arr {7, 3, 2};
+    hpr::StaticArray<float, 4> sarr {arr, 5};
+    hpr::StaticArray<float, 4> sarr2 {7, 3, 2, 5};
+    //hpr::StaticArray<float, 0> sarr4;
+
+    EXPECT_EQ(sarr, sarr2);
+    //EXPECT_EQ(sarr4.data(), nullptr);
+    //EXPECT_EQ(sarr4.is_empty(), true);
+    EXPECT_EQ(sarr.size(), 4);
+    EXPECT_EQ(sarr2.size(), 4);
+    EXPECT_EQ(sarr[0], 7);
+    EXPECT_EQ(sarr[1], 3);
+    EXPECT_EQ(sarr[2], 2);
+    EXPECT_EQ(sarr[3], 5);
+}
+
+#include <hpr/container/dynamic_array.hpp>
+
+TEST(containers, DynamicArray)
+{
+    hpr::DynamicArray<float> arr {1, 3, 2};
+    hpr::DynamicArray<float> arr2 {1, 3, 2};
+    EXPECT_EQ(arr, arr2);
+    EXPECT_TRUE(arr == arr2);
+    arr.remove(1);
+    EXPECT_EQ(arr, hpr::darray<float>({1, 2}));
+    auto iter = arr2.begin();
+    ++iter;
+    arr2.remove(iter);
+    EXPECT_EQ(arr2, hpr::darray<float>({1, 2}));
+
+    hpr::DynamicArray<float> arr3 {1, 3, 0, 2, 9, 0, 5};
+    arr3.remove([](float num) { return num == 0; });
+    EXPECT_EQ(arr3, hpr::darray<float>({1, 3, 2, 9, 5}));
+    EXPECT_EQ(arr3.size(), 5);
+    arr3.insert(3, 19);
+    EXPECT_EQ(arr3.size(), 6);
+    EXPECT_EQ(arr3[3], 19);
+    EXPECT_EQ(arr3.pop(), 5);
+    EXPECT_EQ(arr3.size(), 5);
+    arr3.resize(3);
+    EXPECT_EQ(arr3.size(), 3);
+    EXPECT_EQ(arr3.capacity(), 3);
+    EXPECT_EQ(arr3.back(), 2);
+    arr3.resize(4);
+    EXPECT_EQ(arr3.back(), 2);
+    arr3.push(17);
+    arr3.push(14);
+    EXPECT_EQ(arr3.back(), 14);
+    EXPECT_EQ(arr3.size(), 5);
+
+    arr3.clear();
+    EXPECT_EQ(arr3.is_empty(), true);
+
+    hpr::DynamicArray<float*> arr4;
+    arr4.push(new float(5));
+    arr4.push(new float(7));
+    arr4.push(new float(9));
+    EXPECT_EQ(*arr4[0], 5.f);
+    EXPECT_EQ(*arr4[2], 9.f);
+    for (auto& v : arr4)
+        delete v;
+
+
+}
+
+#include <hpr/container/tree_node.hpp>
+
+TEST(containers, TreeNode)
+{
+    hpr::TreeNode<float> node1 (5);
+    //std::shared_ptr<hpr::TreeNode<float>> ptr {node1};
+    hpr::TreeNode<float> node2 {7, {&node1}};
+    hpr::TreeNode<float> node3 {9, {&node2, &node1}};
+    hpr::TreeNode<float> node4 {11, {&node3}};
+
+    EXPECT_EQ(*node1.data(), 5);
+    //EXPECT_EQ(*node1->data(), *node2.ancestor()->data());
+
+    hpr::darray<hpr::TreeNode<float>*> tr = node1.traverse_descendants();
+    EXPECT_EQ(tr.size(), 0);
+    hpr::darray<hpr::TreeNode<float>*> tr2 = node4.traverse_descendants();
+    EXPECT_EQ(tr2.size(), 4);
+    hpr::darray<hpr::TreeNode<float>*> tr3 = node1.traverse_ancestors();
+    EXPECT_EQ(tr3.size(), 2); // node1 has changed ancestor
+    //
+}
+
+int main(int argc, char **argv)
+{
+    testing::InitGoogleTest(&argc, argv);
+    return RUN_ALL_TESTS();
+}
\ No newline at end of file
diff --git a/source/hpr/tests/test_math.cpp b/source/hpr/tests/test_math.cpp
new file mode 100644
index 0000000..038be15
--- /dev/null
+++ b/source/hpr/tests/test_math.cpp
@@ -0,0 +1,66 @@
+#include <gtest/gtest.h>
+
+#include <hpr/math/vector.hpp>
+#include <hpr/math/matrix.hpp>
+#include <hpr/math/quaternion.hpp>
+
+
+TEST(math, Vector)
+{
+    hpr::vec3 v1 {1, 3, 2};
+    hpr::vec3 v2 {5, 7, -1};
+    hpr::vec2 v31 {13, -2};
+    hpr::vec3 v32 (v31, 9);
+    EXPECT_EQ((hpr::vec3(v2.begin(), v2.end())), v2);
+    EXPECT_EQ(v32, hpr::vec3(13, -2, 9));
+    EXPECT_EQ(-v1, hpr::vec3(-1, -3, -2));
+    EXPECT_EQ(v1 + v2, hpr::vec3(6, 10, 1));
+    EXPECT_EQ(v1 - v2, hpr::vec3(-4, -4, 3));
+    EXPECT_EQ((hpr::dot(v1, v2) ), 24);
+    EXPECT_EQ((hpr::cross(hpr::vec3(1, 0, 0), hpr::vec3(0, 1, 0))), hpr::vec3(0, 0, 1));
+    EXPECT_EQ((hpr::angle(hpr::vec3(1, 0, 0), hpr::vec3(0, 0, 1))), hpr::pi() * 0.5);
+    EXPECT_EQ((hpr::normalize(hpr::vec3(1, 1, 1))), hpr::vec3(0.5773502691896258, 0.5773502691896258, 0.5773502691896258));
+}
+
+TEST(math, Matrix)
+{
+    hpr::mat2 m1;
+    hpr::vec4 v1;
+    EXPECT_FALSE(v1 == m1);
+    EXPECT_FALSE(m1 == v1);
+    hpr::mat2 m2 {3, 2, 7, 4};
+    hpr::vec2 v2 {2, 4};
+    EXPECT_EQ(m2.col(1), v2);
+    hpr::vec2 v3 {13, 51};
+    m2.col(1, v3);
+    EXPECT_EQ(m2.col(1), v3);
+    hpr::mat3 m4 {1, 2, 3, 4, 5, 6, 7, 8, 9};
+    hpr::mat2 m41 {5, 6, 8, 9};
+    //EXPECT_EQ(minor(m41, 0, 0), 9);
+    hpr::mat2 m5 {1, 2, 3, 4};
+
+    EXPECT_EQ(det(m4), 0);
+    EXPECT_EQ(hpr::det(m4), 0);
+    EXPECT_EQ(det(hpr::mat3(-9, 23, 3, 5, 5, 6, 7, -3, 9)), -786);
+
+    hpr::mat2 m6 {2, 1, 7, 4};
+    EXPECT_EQ(det(m6), 1);
+    EXPECT_EQ(adj(m6), hpr::mat2(4, -1, -7, 2));
+    EXPECT_EQ(inv(m6), hpr::mat2(4, -1, -7, 2));
+    EXPECT_EQ(det(hpr::mat3(1, 0, 0, 0, 1, 0, 0, 0, 1)), 1.);
+    //EXPECT_EQ(m4.det(), 0);
+}
+
+TEST(math, Quaternion)
+{
+    hpr::quat q;
+    hpr::vec3 np = hpr::rotate(hpr::vec3(0, 1, 0), {1, 0, 0}, hpr::pi() * 0.5);
+    EXPECT_TRUE(hpr::all(hpr::equal(np, hpr::vec3(0, 0, 1))));
+    EXPECT_EQ(hpr::norm(hpr::quat(np)), 1);
+}
+
+int main(int argc, char **argv)
+{
+    testing::InitGoogleTest(&argc, argv);
+    return RUN_ALL_TESTS();
+}
\ No newline at end of file
diff --git a/source/hpr/tests/test_numeric.cpp b/source/hpr/tests/test_numeric.cpp
new file mode 100644
index 0000000..4c154d4
--- /dev/null
+++ b/source/hpr/tests/test_numeric.cpp
@@ -0,0 +1,45 @@
+#include <gtest/gtest.h>
+
+#include <hpr/numeric.hpp>
+#include <complex>
+
+
+TEST(math, Scalar)
+{
+    EXPECT_EQ(hpr::scalar(5) + hpr::scalar(7), hpr::scalar(12));
+#ifdef HPR_SCALAR_IMPLEMENTATION
+    EXPECT_TRUE(std::is_floating_point_v<hpr::Scalar<double>>);
+    EXPECT_TRUE(std::is_arithmetic_v<hpr::Scalar<double>>);
+#endif
+    //EXPECT_EQ(5.f, hpr::Scalar<double>(5));
+    EXPECT_EQ(hpr::rad(180), hpr::pi());
+    EXPECT_EQ(hpr::deg(hpr::pi()), 180);
+    EXPECT_EQ(hpr::cos(0), 1);
+    EXPECT_EQ(hpr::sin(0), 0);
+    EXPECT_EQ(hpr::abs(hpr::scalar(-1)), 1);
+    EXPECT_EQ(hpr::pow(2, 2), 4);
+    EXPECT_EQ(hpr::precision(), static_cast<hpr::scalar>(1e-15));
+
+    EXPECT_TRUE(typeid(static_cast<float>(hpr::scalar(5))) == typeid(float));
+    EXPECT_FALSE(!hpr::scalar(-1.));
+
+    std::stringstream oss;
+    oss << hpr::cos(0);
+    EXPECT_EQ(oss.str(), "1");
+    hpr::scalar s;
+    oss >> s;
+    EXPECT_EQ(s, 1);
+
+    EXPECT_TRUE(hpr::equal(5.5453535353535395818593, 5.5453535353535395817592, 1e-18));
+    EXPECT_EQ(hpr::min(7., 5), 5);
+    EXPECT_EQ(hpr::max(7., 5), 7);
+    EXPECT_EQ(hpr::clip(7., 5, 10), 7);
+    EXPECT_EQ(hpr::clip(1., 5, 10), 5);
+    EXPECT_EQ(hpr::clip(72., 5, 10), 10);
+}
+
+int main(int argc, char **argv)
+{
+    testing::InitGoogleTest(&argc, argv);
+    return RUN_ALL_TESTS();
+}
\ No newline at end of file
diff --git a/source/hpr/xmake.lua b/source/hpr/xmake.lua
index c3e2cab..3e040de 100644
--- a/source/hpr/xmake.lua
+++ b/source/hpr/xmake.lua
@@ -4,12 +4,14 @@ target("hpr")
 
     set_policy("build.c++.modules", false)
 
+    add_options("scalar", "scalar_implementation")
 
     add_includedirs("..", { public = true })
 
     add_headerfiles("../(hpr/*.hpp)")
+    add_headerfiles("../(hpr/*/*.hpp)")
 
-    add_files("container.cpp")
+    add_files("hpr.cpp")
 
     -- export all symbols for windows/dll
     if is_plat("windows") and is_kind("shared") then
@@ -30,4 +32,26 @@ target("hpr")
             set_configvar("HPR_HAVE_PACKAGE_" .. name:upper(), 1)
         end
     end
-target_end()
\ No newline at end of file
+target_end()
+
+
+for _, file in ipairs(os.files("tests/test_*.cpp")) do
+    local name = path.basename(file)
+    target(name)
+        set_kind("binary")
+        set_default(false)
+        set_languages("c++20")
+
+        set_policy("build.c++.modules", false)
+
+        add_options("scalar", "scalar_implementation")
+
+        add_packages("gtest")
+
+        add_includedirs("..", { public = true })
+
+        add_files(file)
+
+        add_tests("default")
+    target_end()
+end
diff --git a/source/xmake.lua b/source/xmake.lua
index 4285366..704f1eb 100644
--- a/source/xmake.lua
+++ b/source/xmake.lua
@@ -1,16 +1,23 @@
 
 -- options: general
 option("creator")
-    set_category("option")
+    set_category("General")
     set_default(false)
     set_showmenu(true)
-    set_description("Enable creator application")
+    set_description("Build creator application")
+option_end()
+
+option("tests")
+    set_category("General")
+    set_default(false)
+    set_showmenu(true)
+    set_description("Build tests")
 option_end()
 
 -- options: modules
 for _, name in ipairs({"csg", "mesh"}) do
     option(name)
-        set_category("module")
+        set_category("Modules")
         set_default(true)
         set_showmenu(true)
         set_description(format("Module %s", name))
@@ -27,6 +34,30 @@ if has_config("csg") then
     add_requires("conan::opencascade/7.6.2")
 end
 
+if has_config("tests") then
+    add_requires("gtest")
+end
+
+--
+
+option("scalar_implementation")
+    set_category("Internal")
+    set_default(true)
+    set_description("Use internal scalar implementation")
+    add_defines("HPR_SCALAR_IMPLEMENTATION")
+option_end()
+
+option("scalar")
+    set_category("Internal")
+    set_values("float", "double", "long double")
+    set_default("float")
+    set_description("Internal scalar value type")
+    after_check(function(option)
+        option:add("defines", format("HPR_SCALAR=%s", option:value()))
+    end)
+option_end()
+
+
 -- includes project dirs
 includes("hpr")
 if has_config("creator") then