diff --git a/CMakeLists.txt b/CMakeLists.txt
index 59f422d..a8a307b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,14 +1,14 @@
-##======================================================================================================================
+##==============================================================================
 ##  ROTGEN - Runtime Overlay for Eigen
 ##  Copyright : CODE RECKONS
 ##  SPDX-License-Identifier: BSL-1.0
-##======================================================================================================================
+##==============================================================================
 cmake_minimum_required(VERSION 3.22)
 enable_testing()
 
-##======================================================================================================================
+##==============================================================================
 ## Setup project
-##======================================================================================================================
+##==============================================================================
 set(ROTGEN_MAJOR_VERSION 0)
 set(ROTGEN_MINOR_VERSION 0)
 set(ROTGEN_PATCH_VERSION 1)
@@ -18,16 +18,16 @@ set(PROJECT_VERSION   ${ROTGEN_VERSION})
 project(ROTGEN VERSION ${PROJECT_VERSION} DESCRIPTION "Runtime Overlay for Eigen" LANGUAGES CXX)
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${ROTGEN_SOURCE_DIR}/cmake )
 
-##======================================================================================================================
+##==============================================================================
 ## Prevent in-source build
-##======================================================================================================================
+##==============================================================================
 if (${PROJECT_SOURCE_DIR} STREQUAL ${PROJECT_BINARY_DIR})
   message(FATAL_ERROR "[${PROJECT_NAME}]: In-source build is not supported")
 endif()
 
-##======================================================================================================================
+##==============================================================================
 ## Handle options
-##======================================================================================================================
+##==============================================================================
 option(ROTGEN_FORCE_DYNAMIC                 "Force dynamic configuration for rotgen"  OFF)
 option(ROTGEN_ENABLE_EXPRESSION_TEMPLATES    "Enable expression templates"            OFF)
 set(ROTGEN_MAX_SIZE "" CACHE STRING "Max matrix number of elements for rotgen (integer)")
@@ -35,13 +35,14 @@ set(ROTGEN_MAX_SIZE "" CACHE STRING "Max matrix number of elements for rotgen (i
 include(${ROTGEN_SOURCE_DIR}/cmake/options.cmake)
 handle_option(ROTGEN_COMPILE_DEFS)
 
-##======================================================================================================================
+##==============================================================================
 ## Sources & Public Headers lists
-##======================================================================================================================
+##==============================================================================
 if(ROTGEN_FORCE_DYNAMIC)
   set ( SOURCES
         src/map/impl.cpp
         src/matrix/impl.cpp
+        src/quaternion/impl.cpp
         src/block/impl.cpp
         src/svd/impl.cpp
         src/info.cpp
@@ -53,14 +54,14 @@ else()
       )
 endif()
 
-##======================================================================================================================
+##==============================================================================
 ## Setup the library's dependencies
-##======================================================================================================================
+##==============================================================================
 find_package (Eigen3 3.4 REQUIRED NO_MODULE)
 
-##======================================================================================================================
+##==============================================================================
 ## Setup the library's build
-##======================================================================================================================
+##==============================================================================
 if(NOT MSVC)
   set(CMAKE_CXX_VISIBILITY_PRESET hidden)
 endif()
@@ -71,6 +72,11 @@ set_target_properties(rotgen PROPERTIES VERSION       ${PROJECT_VERSION})
 set_target_properties(rotgen PROPERTIES SOVERSION     ${PROJECT_VERSION_MAJOR})
 
 target_compile_options(rotgen PUBLIC -std=c++20 -Werror -Wall -Wextra -Wshadow -Wunused-variable)
+
+if(CMAKE_CXX_COMPILER_ID STREQUAL GNU)
+  target_compile_options(rotgen PUBLIC -Wno-array-bounds)
+endif()
+
 target_compile_definitions(rotgen PUBLIC ${ROTGEN_COMPILE_DEFS})
 add_library(rotgen::rotgen ALIAS rotgen)
 
@@ -81,12 +87,12 @@ target_include_directories(rotgen PUBLIC
 
 target_link_libraries (rotgen PRIVATE Eigen3::Eigen)
 
-##======================================================================================================================
+##==============================================================================
 ## Setup the library's installation target
-##======================================================================================================================
+##==============================================================================
 include(${ROTGEN_SOURCE_DIR}/cmake/config/rotgen-install.cmake)
 
-##======================================================================================================================
+##==============================================================================
 ## Setup the library's Tests
-##======================================================================================================================
+##==============================================================================
 add_subdirectory(test)
diff --git a/README.md b/README.md
index bc0a964..d19fcd2 100644
--- a/README.md
+++ b/README.md
@@ -7,6 +7,7 @@
   + Free function versions of some member function from Eigen to promote a more generic programming approach.
 
 **ROTGEN** depends on [Eigen v3.4.0](https://eigen.tuxfamily.org/index.php?title=Main_Page).
+Documentation for Eigen v3.4.0 is available here: https://libeigen.gitlab.io/eigen/docs-3.4/.
 
 ## Installation
 After cloning this repository, you can setup **ROTGEN** via CMake using the following options:
@@ -51,6 +52,7 @@ The following containers are provided
   + `rotgen::matrix`, which wraps `Eigen::Matrix`.
   + `rotgen::block`, which wraps `Eigen::Block`.
   + `rotgen::map`, which wraps `Eigen::Map`.
+  + `rotgen::quaternion`, which wraps `Eigen::Quaternion`.
   + `rotgen::ref`, which wraps `Eigen::Ref`.
 
 In dynamic mode, `rotgen::matrix`,`rotgen::block` and `rotgen::map` use a PIMPL based implementation to garantee ABI stability.
diff --git a/include/rotgen/algebra/qr.hpp b/include/rotgen/algebra/qr.hpp
index a54d669..bcd45e4 100644
--- a/include/rotgen/algebra/qr.hpp
+++ b/include/rotgen/algebra/qr.hpp
@@ -7,6 +7,8 @@
 //==================================================================================================
 #pragma once
 
+#include <rotgen/container/ref.hpp>
+
 namespace rotgen::solver
 {
   template<typename X, typename M, typename RHS>
diff --git a/include/rotgen/algebra/svd/fixed.hpp b/include/rotgen/algebra/svd/fixed.hpp
index 82b9369..5adcbaf 100644
--- a/include/rotgen/algebra/svd/fixed.hpp
+++ b/include/rotgen/algebra/svd/fixed.hpp
@@ -7,6 +7,11 @@
 //==================================================================================================
 #pragma once
 
+#include <rotgen/detail/helpers.hpp>
+
+#include <rotgen/config.hpp>
+#include <rotgen/container/matrix.hpp>
+
 #include <Eigen/Dense>
 #include <Eigen/SVD>
 
diff --git a/include/rotgen/alias.hpp b/include/rotgen/alias.hpp
index 297f8d8..abd1efd 100644
--- a/include/rotgen/alias.hpp
+++ b/include/rotgen/alias.hpp
@@ -5,6 +5,11 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
+
+#pragma once
+
+#include <rotgen/container/matrix.hpp>
+
 namespace rotgen
 {
   // matrix aliases (float and double, fixed and dynamic)
diff --git a/include/rotgen/container.hpp b/include/rotgen/container.hpp
index 9f9adff..b6e4f8d 100644
--- a/include/rotgen/container.hpp
+++ b/include/rotgen/container.hpp
@@ -7,7 +7,8 @@
 //==================================================================================================
 #pragma once
 
-#include <rotgen/container/matrix.hpp>
 #include <rotgen/container/block.hpp>
 #include <rotgen/container/map.hpp>
+#include <rotgen/container/matrix.hpp>
+#include <rotgen/container/quaternion.hpp>
 #include <rotgen/container/ref.hpp>
diff --git a/include/rotgen/container/block/fixed.hpp b/include/rotgen/container/block/fixed.hpp
index 64016a5..9263b2d 100644
--- a/include/rotgen/container/block/fixed.hpp
+++ b/include/rotgen/container/block/fixed.hpp
@@ -7,10 +7,12 @@
 //==================================================================================================
 #pragma once
 
+#include <rotgen/detail/helpers.hpp>
+
+#include <rotgen/container/matrix.hpp>
 #include <rotgen/container/strides.hpp>
 
 #include <Eigen/Dense>
-#include <iostream>
 
 namespace rotgen
 {
diff --git a/include/rotgen/container/map/dynamic.hpp b/include/rotgen/container/map/dynamic.hpp
index ffb2223..e9c7b55 100644
--- a/include/rotgen/container/map/dynamic.hpp
+++ b/include/rotgen/container/map/dynamic.hpp
@@ -11,6 +11,7 @@
 
 #include <rotgen/concepts.hpp>
 #include <rotgen/container/map/dynamic/impl.hpp>
+#include <rotgen/container/matrix.hpp>
 
 #include <cassert>
 
diff --git a/include/rotgen/container/matrix/dynamic.hpp b/include/rotgen/container/matrix/dynamic.hpp
index 8c96261..ae0d1a6 100644
--- a/include/rotgen/container/matrix/dynamic.hpp
+++ b/include/rotgen/container/matrix/dynamic.hpp
@@ -8,7 +8,6 @@
 #pragma once
 
 #include <rotgen/detail/helpers.hpp>
-
 #include <rotgen/concepts.hpp>
 #include <rotgen/container/matrix/dynamic/impl.hpp>
 
@@ -445,11 +444,7 @@ namespace rotgen
 
     parent& base() { return static_cast<parent&>(*this); }
 
-    parent const& base() const
-    {
-      return static_cast<parent const&>(*this);
-      ;
-    }
+    parent const& base() const { return static_cast<parent const&>(*this); }
   };
 
   template<typename S, int R, int C, int O, int MR, int MC>
diff --git a/include/rotgen/container/matrix/fixed.hpp b/include/rotgen/container/matrix/fixed.hpp
index a74fc24..6c11b67 100644
--- a/include/rotgen/container/matrix/fixed.hpp
+++ b/include/rotgen/container/matrix/fixed.hpp
@@ -9,8 +9,9 @@
 
 #include <rotgen/detail/helpers.hpp>
 
+#include <rotgen/concepts.hpp>
+
 #include <Eigen/Dense>
-#include <iostream>
 
 namespace rotgen
 {
diff --git a/include/rotgen/container/quaternion.hpp b/include/rotgen/container/quaternion.hpp
new file mode 100644
index 0000000..4625b77
--- /dev/null
+++ b/include/rotgen/container/quaternion.hpp
@@ -0,0 +1,19 @@
+//==============================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==============================================================================
+
+#pragma once
+
+// https://libeigen.gitlab.io/eigen/docs-3.4/classEigen_1_1Quaternion.html
+
+#include <rotgen/config.hpp>
+
+#if defined(ROTGEN_FORCE_DYNAMIC)
+#include <rotgen/container/quaternion/dynamic.hpp>
+#else
+#include <rotgen/container/quaternion/fixed.hpp>
+#endif
diff --git a/include/rotgen/container/quaternion/dynamic.hpp b/include/rotgen/container/quaternion/dynamic.hpp
new file mode 100644
index 0000000..30446d1
--- /dev/null
+++ b/include/rotgen/container/quaternion/dynamic.hpp
@@ -0,0 +1,84 @@
+//==============================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==============================================================================
+
+#pragma once
+
+#include <rotgen/detail/helpers.hpp>
+
+#include <rotgen/config.hpp>
+#include <rotgen/container/matrix.hpp>
+#include <rotgen/container/quaternion/dynamic/impl.hpp>
+#include <rotgen/container/ref.hpp>
+
+namespace rotgen
+{
+  template<typename Scalar, int Options>
+  class quaternion : public find_quaternion<Scalar>
+  {
+  public:
+    using parent = find_quaternion<Scalar>;
+
+    using CoeffReturnType = Scalar;
+    using NonConstCoeffReturnType = Scalar&;
+
+    // Constructors
+
+    quaternion() : parent() {}
+
+    quaternion(parent const& other) : parent(other) {}
+
+    quaternion(parent&& other) : parent(other) {}
+
+    quaternion(Scalar w, Scalar x, Scalar y, Scalar z) : parent(w, x, y, z) {}
+
+    // Coefficient accessors
+
+    CoeffReturnType w() const { return parent::w(); }
+
+    NonConstCoeffReturnType w() { return parent::w(); }
+
+    CoeffReturnType x() const { return parent::x(); }
+
+    NonConstCoeffReturnType x() { return parent::x(); }
+
+    CoeffReturnType y() const { return parent::y(); }
+
+    NonConstCoeffReturnType y() { return parent::y(); }
+
+    CoeffReturnType z() const { return parent::z(); }
+
+    NonConstCoeffReturnType z() { return parent::z(); }
+
+    parent& base() { return static_cast<parent&>(*this); }
+
+    parent const& base() const { return static_cast<parent const&>(*this); }
+
+    /// Multiplication method, implemented not as operator* to bypass common
+    /// operator* implementation and avoid Eigen's quirky behavior when
+    /// multiplying quaternions and matrices.
+    template<typename T> auto mul(ref<matrix<T, 3, 1> const> mat) const
+    {
+      return parent::mul(mat.base());
+    }
+  };
+
+  /// Wraps `Eigen::AngleAxis`
+  template<typename Scalar>
+  quaternion<Scalar> angleAxis(Scalar theta,
+                               rotgen::matrix<Scalar, 3, 1> const& axis)
+  {
+    return quaternion<Scalar>::angleAxis(theta, axis);
+  }
+
+  template<typename Scalar>
+  auto operator*(quaternion<Scalar> const& q,
+                 ref<matrix<Scalar, 3, 1> const> const& mat)
+  {
+    return q.mul(mat);
+  }
+}
diff --git a/include/rotgen/container/quaternion/dynamic/impl.hpp b/include/rotgen/container/quaternion/dynamic/impl.hpp
new file mode 100644
index 0000000..1e45209
--- /dev/null
+++ b/include/rotgen/container/quaternion/dynamic/impl.hpp
@@ -0,0 +1,163 @@
+//==============================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==============================================================================
+
+#pragma once
+
+#include <rotgen/detail/generators.hpp>
+
+#include <rotgen/config.hpp>
+#include <rotgen/container/map.hpp>
+#include <rotgen/container/matrix.hpp>
+
+namespace rotgen
+{
+  // Only example to support:
+
+  // ConstColVectorRef<3>& T_;
+  // Eigen::Quaternion<Real> Q;
+  // Q = Eigen::AngleAxis<Real>(theta_, N_);
+  // T_ = Q * T_;
+
+  /// Hard-wired implementation of quaternion for 32-bit floats
+  class ROTGEN_EXPORT quaternion_impl_32
+  {
+  private:
+    struct payload;
+    std::unique_ptr<payload> storage_;
+
+  public:
+    // Construction, destruction, assignment
+
+    quaternion_impl_32();
+    quaternion_impl_32(quaternion_impl_32 const&);
+    quaternion_impl_32(quaternion_impl_32&&) noexcept;
+    quaternion_impl_32& operator=(quaternion_impl_32 const&);
+    quaternion_impl_32& operator=(quaternion_impl_32&&) noexcept;
+    ~quaternion_impl_32();
+
+    quaternion_impl_32(float w, float x, float y, float z);
+
+    // Comparison
+
+    friend ROTGEN_EXPORT bool operator==(quaternion_impl_32 const&,
+                                         quaternion_impl_32 const&);
+    friend ROTGEN_EXPORT bool operator!=(quaternion_impl_32 const&,
+                                         quaternion_impl_32 const&);
+
+    // Data and payload access
+
+    std::unique_ptr<payload>& storage() { return storage_; }
+
+    std::unique_ptr<payload> const& storage() const { return storage_; }
+
+    float const* data() const;
+    float* data();
+
+    // Coefficient accessors
+
+    float w() const;
+    float& w();
+    float x() const;
+    float& x();
+    float y() const;
+    float& y();
+    float z() const;
+    float& z();
+
+    // Arithmetic operators
+
+    /// Multiplication method, implemented not as operator* to bypass common
+    /// operator* implementation and avoid Eigen's quirky behavior when
+    /// multiplying quaternions and matrices.
+    matrix<float, 3, 1> mul(map<matrix<float, 3, 1> const>) const;
+
+    // Generators
+
+    static quaternion_impl_32 angleAxis(float theta,
+                                        matrix<float, 3, 1> const& axis);
+  };
+
+  /// Hard-wired implementation of quaternion for 64-bit floats
+  class ROTGEN_EXPORT quaternion_impl_64
+  {
+  private:
+    struct payload;
+    std::unique_ptr<payload> storage_;
+
+  public:
+    // Construction, destruction, assignment
+
+    quaternion_impl_64();
+    quaternion_impl_64(quaternion_impl_64 const&);
+    quaternion_impl_64(quaternion_impl_64&&) noexcept;
+    quaternion_impl_64& operator=(quaternion_impl_64 const&);
+    quaternion_impl_64& operator=(quaternion_impl_64&&) noexcept;
+    ~quaternion_impl_64();
+
+    quaternion_impl_64(double w, double x, double y, double z);
+
+    // Comparison
+
+    friend ROTGEN_EXPORT bool operator==(quaternion_impl_64 const&,
+                                         quaternion_impl_64 const&);
+    friend ROTGEN_EXPORT bool operator!=(quaternion_impl_64 const&,
+                                         quaternion_impl_64 const&);
+
+    // Data and payload access
+
+    std::unique_ptr<payload>& storage() { return storage_; }
+
+    std::unique_ptr<payload> const& storage() const { return storage_; }
+
+    double const* data() const;
+    double* data();
+
+    // Coefficient accessors
+
+    double w() const;
+    double& w();
+    double x() const;
+    double& x();
+    double y() const;
+    double& y();
+    double z() const;
+    double& z();
+
+    // Arithmetic operators
+
+    /// Multiplication method, implemented not as operator* to bypass common
+    /// operator* implementation and avoid Eigen's quirky behavior when
+    /// multiplying quaternions and matrices.
+    matrix<double, 3, 1> mul(map<matrix<double, 3, 1> const>) const;
+
+    // Generators
+
+    static quaternion_impl_64 angleAxis(double theta,
+                                        matrix<double, 3, 1> const& axis);
+  };
+
+  // find_quaternion
+
+  template<typename Scalar> struct _find_quaternion_impl;
+
+  template<> struct _find_quaternion_impl<float>
+  {
+    using type = quaternion_impl_32;
+  };
+
+  template<> struct _find_quaternion_impl<double>
+  {
+    using type = quaternion_impl_64;
+  };
+
+  /// Resolves to the right type of quaternion depending on the input type
+  /// (float or double).
+  template<typename Scalar>
+  using find_quaternion = typename _find_quaternion_impl<Scalar>::type;
+
+}
diff --git a/include/rotgen/container/quaternion/fixed.hpp b/include/rotgen/container/quaternion/fixed.hpp
new file mode 100644
index 0000000..1705c3e
--- /dev/null
+++ b/include/rotgen/container/quaternion/fixed.hpp
@@ -0,0 +1,98 @@
+//==============================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==============================================================================
+
+#pragma once
+
+#include <rotgen/detail/helpers.hpp>
+
+#include <rotgen/concepts.hpp>
+#include <rotgen/config.hpp>
+#include <rotgen/container/matrix.hpp>
+#include <rotgen/container/ref.hpp>
+
+#include <Eigen/Dense>
+
+namespace rotgen
+{
+  template<typename Scalar, int Options>
+  class quaternion : private Eigen::Quaternion<Scalar, Options>
+  {
+  public:
+    using rotgen_tag = void;
+    using parent = Eigen::Quaternion<Scalar, Options>;
+    using value_type = Scalar;
+
+    quaternion() : parent() {}
+
+    quaternion(Scalar w_, Scalar x_, Scalar y_, Scalar z_)
+      : parent(w_, x_, y_, z_)
+    {
+    }
+
+    quaternion(concepts::entity auto const& other) : parent(other.base()) {}
+
+    template<typename OtherDerived>
+    quaternion(Eigen::QuaternionBase<OtherDerived> const& other) : parent(other)
+    {
+    }
+
+    template<typename OtherDerived>
+    quaternion(Eigen::EigenBase<OtherDerived> const& other) : parent(other)
+    {
+    }
+
+    template<typename OtherDerived>
+    quaternion& operator=(Eigen::QuaternionBase<OtherDerived> const& other)
+    {
+      parent::operator=(other);
+      return *this;
+    }
+
+    template<typename OtherDerived>
+    quaternion& operator=(Eigen::EigenBase<OtherDerived> const& other)
+    {
+      parent::operator=(other);
+      return *this;
+    }
+
+    quaternion& operator=(concepts::entity auto const& other)
+    {
+      parent::operator=(other.base());
+      return *this;
+    }
+
+    quaternion(quaternion const& other) = default;
+    quaternion(quaternion&& other) = default;
+    quaternion& operator=(quaternion const&) = default;
+    quaternion& operator=(quaternion&&) = default;
+
+    parent& base() { return static_cast<parent&>(*this); }
+
+    parent const& base() const { return static_cast<parent const&>(*this); }
+
+    using parent::w;
+    using parent::x;
+    using parent::y;
+    using parent::z;
+  };
+
+  template<typename Scalar>
+  auto operator*(quaternion<Scalar> const& q,
+                 ref<matrix<Scalar, 3, 1> const> const& mat)
+  {
+    return matrix<Scalar, 3, 1>{q.base() * mat.base()};
+  }
+
+  template<typename Scalar>
+  quaternion<Scalar> angleAxis(Scalar theta, rotgen::concepts::entity auto mat)
+  {
+    Eigen::Quaternion<Scalar> result;
+    result = Eigen::AngleAxis<Scalar>(theta, mat.base());
+    return {result};
+  }
+}
diff --git a/include/rotgen/container/ref.hpp b/include/rotgen/container/ref.hpp
index f59fa47..2d99197 100644
--- a/include/rotgen/container/ref.hpp
+++ b/include/rotgen/container/ref.hpp
@@ -8,6 +8,13 @@
 #pragma once
 
 #include <rotgen/config.hpp>
+#include <rotgen/container/block.hpp>
+#include <rotgen/container/map.hpp>
+#include <rotgen/format.hpp>
+
+#include <cassert>
+#include <iostream>
+#include <type_traits>
 
 #if defined(ROTGEN_FORCE_DYNAMIC)
 #include <rotgen/container/ref/dynamic.hpp>
diff --git a/include/rotgen/container/ref/generalize.hpp b/include/rotgen/container/ref/generalize.hpp
index 021646b..b5ad520 100644
--- a/include/rotgen/container/ref/generalize.hpp
+++ b/include/rotgen/container/ref/generalize.hpp
@@ -7,10 +7,14 @@
 //==================================================================================================
 #pragma once
 
+#include <rotgen/config.hpp>
+
 #include <type_traits>
 
 namespace rotgen
 {
+  template<typename Scalar, int Options = AutoAlign> class quaternion;
+
   //-------------------------------------------------------------------------------------------
   // Convert entity/eigen types to a proper ref so we can write less function
   // overloads
@@ -46,6 +50,16 @@ namespace rotgen
     using type = ref<T, O, S>;
   };
 
+  template<typename T> struct generalize<quaternion<T>>
+  {
+    using type = quaternion<T>;
+  };
+
+  template<typename T> struct generalize<quaternion<T> const>
+  {
+    using type = quaternion<T>;
+  };
+
   template<concepts::entity T> typename T::parent& base_of(T& a)
   {
     return a.base();
diff --git a/include/rotgen/container/strides.hpp b/include/rotgen/container/strides.hpp
index f4c8545..303b309 100644
--- a/include/rotgen/container/strides.hpp
+++ b/include/rotgen/container/strides.hpp
@@ -7,6 +7,7 @@
 //==================================================================================================
 #pragma once
 
+#include <rotgen/concepts.hpp>
 #include <rotgen/config.hpp>
 
 #if defined(ROTGEN_FORCE_DYNAMIC)
diff --git a/include/rotgen/detail/payload.hpp b/include/rotgen/detail/payload.hpp
index 96998e0..c590289 100644
--- a/include/rotgen/detail/payload.hpp
+++ b/include/rotgen/detail/payload.hpp
@@ -1,19 +1,49 @@
-//==================================================================================================
+//==============================================================================
 /*
   ROTGEN - Runtime Overlay for Eigen
   Copyright : CODE RECKONS
   SPDX-License-Identifier: BSL-1.0
 */
-//==================================================================================================
+//==============================================================================
 #pragma once
 
 #include <rotgen/detail/payload/block.hpp>
 #include <rotgen/detail/payload/map.hpp>
 #include <rotgen/detail/payload/matrix.hpp>
 
+#include <rotgen/container/quaternion/dynamic/impl.hpp>
+
+#include <Eigen/Core>
+#include <Eigen/Dense>
+#include <cstddef>
+
 namespace rotgen
 {
-  // IO Formatting paylod
+  struct quaternion_impl_32::payload
+  {
+    using data_type = Eigen::Quaternion<float>;
+    data_type data;
+
+    payload() : data() {}
+
+    payload(data_type const& o) : data(o) {}
+
+    payload(float* ptr) : data(ptr) {}
+  };
+
+  struct quaternion_impl_64::payload
+  {
+    using data_type = Eigen::Quaternion<double>;
+    data_type data;
+
+    payload() : data() {}
+
+    payload(data_type const& o) : data(o) {}
+
+    payload(double* ptr) : data(ptr) {}
+  };
+
+  // IO Payload
   struct ioformat::payload
   {
     Eigen::IOFormat instance;
diff --git a/include/rotgen/rotgen.hpp b/include/rotgen/rotgen.hpp
index 2b14aa7..982f41c 100644
--- a/include/rotgen/rotgen.hpp
+++ b/include/rotgen/rotgen.hpp
@@ -7,11 +7,9 @@
 //==================================================================================================
 #pragma once
 
-//clang-format off
-#include <rotgen/config.hpp>
-#include <rotgen/format.hpp>
-#include <rotgen/container.hpp>
-#include <rotgen/alias.hpp>
 #include <rotgen/algebra.hpp>
+#include <rotgen/alias.hpp>
+#include <rotgen/config.hpp>
+#include <rotgen/container.hpp>
+#include <rotgen/format.hpp>
 #include <rotgen/functions.hpp>
-//clang-format on
diff --git a/src/quaternion/impl.cpp b/src/quaternion/impl.cpp
new file mode 100644
index 0000000..7cc9bb9
--- /dev/null
+++ b/src/quaternion/impl.cpp
@@ -0,0 +1,261 @@
+//==============================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==============================================================================
+#include <rotgen/container/quaternion/dynamic/impl.hpp>
+
+#include <rotgen/detail/generators.hpp>
+#include <rotgen/detail/payload.hpp>
+
+#include <rotgen/config.hpp>
+#include <rotgen/container/map.hpp>
+
+#include <Eigen/Dense>
+#include <Eigen/src/Geometry/AngleAxis.h>
+
+namespace rotgen
+{
+
+  // Payload constructor
+
+  quaternion_impl_32::quaternion_impl_32()
+    : storage_(std::make_unique<payload>())
+  {
+  }
+
+  quaternion_impl_64::quaternion_impl_64()
+    : storage_(std::make_unique<payload>())
+  {
+  }
+
+  // Regular copy/move constructors/assigns
+
+  quaternion_impl_32::quaternion_impl_32(quaternion_impl_32&&) noexcept =
+    default;
+  quaternion_impl_64::quaternion_impl_64(quaternion_impl_64&&) noexcept =
+    default;
+
+  quaternion_impl_32::quaternion_impl_32(quaternion_impl_32 const& other)
+    : quaternion_impl_32()
+  {
+    storage_->data = other.storage_->data;
+  }
+
+  quaternion_impl_64::quaternion_impl_64(quaternion_impl_64 const& other)
+    : quaternion_impl_64()
+  {
+    storage_->data = other.storage_->data;
+  }
+
+  quaternion_impl_32& quaternion_impl_32::operator=(
+    quaternion_impl_32 const& other)
+  {
+    if (this != &other) storage_->data = other.storage_->data;
+    return *this;
+  }
+
+  quaternion_impl_64& quaternion_impl_64::operator=(
+    quaternion_impl_64 const& other)
+  {
+    if (this != &other) storage_->data = other.storage_->data;
+    return *this;
+  }
+
+  quaternion_impl_32& quaternion_impl_32::operator=(
+    quaternion_impl_32&&) noexcept = default;
+  quaternion_impl_64& quaternion_impl_64::operator=(
+    quaternion_impl_64&&) noexcept = default;
+
+  quaternion_impl_32::quaternion_impl_32(float w, float x, float y, float z)
+    : storage_(std::make_unique<payload>(payload::data_type{w, x, y, z}))
+  {
+  }
+
+  quaternion_impl_64::quaternion_impl_64(double w, double x, double y, double z)
+    : storage_(std::make_unique<payload>(payload::data_type{w, x, y, z}))
+  {
+  }
+
+  // Destructors
+
+  quaternion_impl_32::~quaternion_impl_32() = default;
+  quaternion_impl_64::~quaternion_impl_64() = default;
+
+  // .data()
+
+  float const* quaternion_impl_32::data() const
+  {
+    return storage_->data.coeffs().data();
+  }
+
+  double const* quaternion_impl_64::data() const
+  {
+    return storage_->data.coeffs().data();
+  }
+
+  float* quaternion_impl_32::data()
+  {
+    return storage_->data.coeffs().data();
+  }
+
+  double* quaternion_impl_64::data()
+  {
+    return storage_->data.coeffs().data();
+  }
+
+  // Coefficient accessors
+
+  float quaternion_impl_32::w() const
+  {
+    return storage_->data.w();
+  }
+
+  float& quaternion_impl_32::w()
+  {
+    return storage_->data.w();
+  }
+
+  float quaternion_impl_32::x() const
+  {
+    return storage_->data.x();
+  }
+
+  float& quaternion_impl_32::x()
+  {
+    return storage_->data.x();
+  }
+
+  float quaternion_impl_32::y() const
+  {
+    return storage_->data.y();
+  }
+
+  float& quaternion_impl_32::y()
+  {
+    return storage_->data.y();
+  }
+
+  float quaternion_impl_32::z() const
+  {
+    return storage_->data.z();
+  }
+
+  float& quaternion_impl_32::z()
+  {
+    return storage_->data.z();
+  }
+
+  double quaternion_impl_64::w() const
+  {
+    return storage_->data.w();
+  }
+
+  double& quaternion_impl_64::w()
+  {
+    return storage_->data.w();
+  }
+
+  double quaternion_impl_64::x() const
+  {
+    return storage_->data.x();
+  }
+
+  double& quaternion_impl_64::x()
+  {
+    return storage_->data.x();
+  }
+
+  double quaternion_impl_64::y() const
+  {
+    return storage_->data.y();
+  }
+
+  double& quaternion_impl_64::y()
+  {
+    return storage_->data.y();
+  }
+
+  double quaternion_impl_64::z() const
+  {
+    return storage_->data.z();
+  }
+
+  double& quaternion_impl_64::z()
+  {
+    return storage_->data.z();
+  }
+
+  // operator==
+
+  ROTGEN_EXPORT bool operator==(quaternion_impl_32 const& lhs,
+                                quaternion_impl_32 const& rhs)
+  {
+    return lhs.storage_->data == rhs.storage_->data;
+  }
+
+  ROTGEN_EXPORT bool operator==(quaternion_impl_64 const& lhs,
+                                quaternion_impl_64 const& rhs)
+  {
+    return lhs.storage_->data == rhs.storage_->data;
+  }
+
+  // operator!=
+
+  ROTGEN_EXPORT bool operator!=(quaternion_impl_32 const& lhs,
+                                quaternion_impl_32 const& rhs)
+  {
+    return lhs.storage_->data != rhs.storage_->data;
+  }
+
+  ROTGEN_EXPORT bool operator!=(quaternion_impl_64 const& lhs,
+                                quaternion_impl_64 const& rhs)
+  {
+    return lhs.storage_->data != rhs.storage_->data;
+  }
+
+  // Arithmetic operators
+
+  matrix<float, 3, 1> quaternion_impl_32::mul(
+    map<matrix<float, 3, 1> const> mat) const
+  {
+    matrix<float, 3, 1> result;
+
+    // TODO: Document Eigen quirk
+    Eigen::Map<Eigen::Matrix<float, 3, 1> const> local(mat.data());
+    result.storage()->data = this->storage()->data * local;
+    return result;
+  }
+
+  matrix<double, 3, 1> quaternion_impl_64::mul(
+    map<matrix<double, 3, 1> const> mat) const
+  {
+    matrix<double, 3, 1> result;
+
+    // TODO: Document Eigen quirk
+    Eigen::Map<Eigen::Matrix<double, 3, 1> const> local(mat.data());
+    result.storage()->data = this->storage()->data * local;
+    return result;
+  }
+
+  // rotgen::quaternion::angleAxis
+
+  quaternion_impl_32 quaternion_impl_32::angleAxis(
+    float theta, matrix<float, 3, 1> const& axis)
+  {
+    quaternion_impl_32 m;
+    m.storage_->data = Eigen::AngleAxis<float>(theta, axis.storage()->data);
+    return m;
+  }
+
+  quaternion_impl_64 quaternion_impl_64::angleAxis(
+    double theta, matrix<double, 3, 1> const& axis)
+  {
+    quaternion_impl_64 m;
+    m.storage_->data = Eigen::AngleAxis<double>(theta, axis.storage()->data);
+    return m;
+  }
+
+}
diff --git a/src/quaternion/model.inl b/src/quaternion/model.inl
new file mode 100644
index 0000000..e69de29
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index ba0a2de..396360d 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -27,6 +27,7 @@ rotgen_glob_unit(QUIET PATTERN "unit/matrix/*.cpp"    INTERFACE rotgen_test)
 rotgen_glob_unit(QUIET PATTERN "unit/block/*.cpp"     INTERFACE rotgen_test)
 rotgen_glob_unit(QUIET PATTERN "unit/map/*.cpp"       INTERFACE rotgen_test)
 rotgen_glob_unit(QUIET PATTERN "unit/meta/*.cpp"      INTERFACE rotgen_test)
+rotgen_glob_unit(QUIET PATTERN "unit/quaternion/*.cpp" INTERFACE rotgen_test)
 rotgen_glob_unit(QUIET PATTERN "unit/functions/*.cpp" INTERFACE rotgen_test)
 
 ##======================================================================================================================
diff --git a/test/unit/block/arithmetic_functions.cpp b/test/unit/block/arithmetic_functions.cpp
index 8dd3128..29b83c3 100644
--- a/test/unit/block/arithmetic_functions.cpp
+++ b/test/unit/block/arithmetic_functions.cpp
@@ -5,10 +5,12 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
-#include "unit/common/arithmetic.hpp"
 #include <rotgen/rotgen.hpp>
 
+#include "unit/common/arithmetic.hpp"
+#include "unit/common/references.hpp"
+#include "unit/tests.hpp"
+
 TTS_CASE_TPL("Test dynamic block transposition-like operations",
              rotgen::tests::types)<typename T, typename O>(
   tts::type<tts::types<T, O>>)
diff --git a/test/unit/block/basic_api.cpp b/test/unit/block/basic_api.cpp
index a0ac871..c74515c 100644
--- a/test/unit/block/basic_api.cpp
+++ b/test/unit/block/basic_api.cpp
@@ -5,9 +5,10 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
 #include <rotgen/rotgen.hpp>
 
+#include "unit/tests.hpp"
+
 // Helper: fill matrix from raw data
 // NB: This function must not be turned into a lambda, otherwise
 // `test-ubuntu-gcc-release` will not pass and the CI will fail.
diff --git a/test/unit/block/cwise.cpp b/test/unit/block/cwise.cpp
index dbfbb9c..c6b521c 100644
--- a/test/unit/block/cwise.cpp
+++ b/test/unit/block/cwise.cpp
@@ -5,10 +5,12 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
-#include "unit/common/cwise.hpp"
 #include <rotgen/rotgen.hpp>
 
+#include "unit/common/cwise.hpp"
+#include "unit/common/references.hpp"
+#include "unit/tests.hpp"
+
 TTS_CASE_TPL("Test dynamic block cwise operations",
              rotgen::tests::types)<typename T, typename O>(
   tts::type<tts::types<T, O>>)
diff --git a/test/unit/block/extract.cpp b/test/unit/block/extract.cpp
index 3b61c5c..236aa01 100644
--- a/test/unit/block/extract.cpp
+++ b/test/unit/block/extract.cpp
@@ -7,6 +7,7 @@
 //==================================================================================================
 #include <rotgen/rotgen.hpp>
 
+#include "unit/common/references.hpp"
 #include "unit/tests.hpp"
 
 template<typename EigenType, typename F>
diff --git a/test/unit/block/generators.cpp b/test/unit/block/generators.cpp
index 786d55d..e6f8fc3 100644
--- a/test/unit/block/generators.cpp
+++ b/test/unit/block/generators.cpp
@@ -8,6 +8,7 @@
 #include <rotgen/rotgen.hpp>
 
 #include "rotgen/functions/functions.hpp"
+#include "unit/common/references.hpp"
 #include "unit/tests.hpp"
 
 void test_value(auto const& matrix,
diff --git a/test/unit/block/norms.cpp b/test/unit/block/norms.cpp
index 6ac2947..a1efed9 100644
--- a/test/unit/block/norms.cpp
+++ b/test/unit/block/norms.cpp
@@ -5,10 +5,12 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
-#include "unit/common/norms.hpp"
 #include <rotgen/rotgen.hpp>
 
+#include "unit/common/norms.hpp"
+#include "unit/common/references.hpp"
+#include "unit/tests.hpp"
+
 TTS_CASE_TPL("Test dynamic block norm operations",
              rotgen::tests::types)<typename T, typename O>(
   tts::type<tts::types<T, O>>)
diff --git a/test/unit/block/operators.cpp b/test/unit/block/operators.cpp
index dffd7d6..e2002c3 100644
--- a/test/unit/block/operators.cpp
+++ b/test/unit/block/operators.cpp
@@ -5,8 +5,10 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
 #include <rotgen/rotgen.hpp>
+
+#include "unit/common/references.hpp"
+#include "unit/tests.hpp"
 #include <Eigen/Dense>
 
 template<typename MatrixType, typename T>
diff --git a/test/unit/common/arithmetic.hpp b/test/unit/common/arithmetic.hpp
index c102294..4899850 100644
--- a/test/unit/common/arithmetic.hpp
+++ b/test/unit/common/arithmetic.hpp
@@ -10,6 +10,8 @@
 #include <rotgen/rotgen.hpp>
 
 #include "tts.hpp"
+#include "unit/common/references.hpp"
+#include "unit/tests.hpp"
 #include <Eigen/Dense>
 
 namespace rotgen::tests
diff --git a/test/unit/common/cwise.hpp b/test/unit/common/cwise.hpp
index 9e903b1..313361c 100644
--- a/test/unit/common/cwise.hpp
+++ b/test/unit/common/cwise.hpp
@@ -10,6 +10,7 @@
 #include <rotgen/rotgen.hpp>
 
 #include "tts.hpp"
+#include "unit/common/references.hpp"
 #include <Eigen/Dense>
 
 namespace rotgen::tests
diff --git a/test/unit/common/norms.hpp b/test/unit/common/norms.hpp
index 9c6fca8..711c190 100644
--- a/test/unit/common/norms.hpp
+++ b/test/unit/common/norms.hpp
@@ -10,6 +10,7 @@
 #include <rotgen/rotgen.hpp>
 
 #include "tts.hpp"
+#include "unit/common/references.hpp"
 #include <Eigen/Dense>
 
 namespace rotgen::tests
diff --git a/test/unit/common/references.hpp b/test/unit/common/references.hpp
index 97c2bd3..08a07d6 100644
--- a/test/unit/common/references.hpp
+++ b/test/unit/common/references.hpp
@@ -1,10 +1,10 @@
-//==================================================================================================
+//==============================================================================
 /*
   ROTGEN - Runtime Overlay for Eigen
   Copyright : CODE RECKONS
   SPDX-License-Identifier: BSL-1.0
 */
-//==================================================================================================
+//==============================================================================
 #pragma once
 
 #include <rotgen/rotgen.hpp>
@@ -59,30 +59,33 @@ namespace rotgen::tests
       for (rotgen::Index c = 0; c < output.cols(); ++c) output(r, c) = fn(r, c);
   }
 
+  auto default_init_function = [](auto row, auto col) {
+    return row + 3 * col - 2.5;
+  };
+
   auto generate_matrix_references()
   {
-    auto fn = [](auto r, auto c) { return r + 3 * c - 2.5; };
     std::vector<rotgen::tests::matrix_descriptor> cases = {
       // Singular matrix
-      {1, 1, fn},
+      {1, 1, default_init_function},
       // Square matrix below MAX_SIZE
-      {3, 3, fn},
+      {3, 3, default_init_function},
       // Square matrix at MAX_SIZE
-      {4, 4, fn},
+      {4, 4, default_init_function},
       // Square matrix above MAX_SIZE
-      {7, 7, fn},
+      {7, 7, default_init_function},
       // Tall matrix below MAX_SIZE
-      {5, 2, fn},
+      {5, 2, default_init_function},
       // Tall matrix at MAX_SIZE
-      {8, 2, fn},
+      {8, 2, default_init_function},
       // Tall matrix above MAX_SIZE
-      {10, 3, fn},
+      {10, 3, default_init_function},
       // Thick matrix below MAX_SIZE
-      {2, 5, fn},
+      {2, 5, default_init_function},
       // Thick matrix at MAX_SIZE
-      {2, 8, fn},
+      {2, 8, default_init_function},
       // Thick matrix above MAX_SIZE
-      {3, 10, fn}};
+      {3, 10, default_init_function}};
 
     return cases;
   }
diff --git a/test/unit/functions/qr.cpp b/test/unit/functions/qr.cpp
index 3b0cd01..9d4f299 100644
--- a/test/unit/functions/qr.cpp
+++ b/test/unit/functions/qr.cpp
@@ -5,9 +5,10 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
 #include <rotgen/rotgen.hpp>
 
+#include "unit/tests.hpp"
+
 TTS_CASE_TPL("System solver using QR",
              rotgen::tests::types)<typename T, typename O>(
   tts::type<tts::types<T, O>>)
diff --git a/test/unit/functions/rowwise.cpp b/test/unit/functions/rowwise.cpp
index 8afdce1..9f67e41 100644
--- a/test/unit/functions/rowwise.cpp
+++ b/test/unit/functions/rowwise.cpp
@@ -7,6 +7,7 @@
 //==================================================================================================
 #include <rotgen/rotgen.hpp>
 
+#include "unit/common/references.hpp"
 #include "unit/tests.hpp"
 
 TTS_CASE_TPL("rowwise API", rotgen::tests::types)<typename T, typename O>(
diff --git a/test/unit/functions/svd.cpp b/test/unit/functions/svd.cpp
index 0e836ee..23671a1 100644
--- a/test/unit/functions/svd.cpp
+++ b/test/unit/functions/svd.cpp
@@ -5,9 +5,10 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
 #include <rotgen/rotgen.hpp>
 
+#include "unit/tests.hpp"
+
 TTS_CASE_TPL("SVD decomposition - Dynamic case",
              rotgen::tests::types)<typename T, typename O>(
   tts::type<tts::types<T, O>>)
@@ -19,7 +20,8 @@ TTS_CASE_TPL("SVD decomposition - Dynamic case",
     rotgen::matrix<T, rotgen::Dynamic, rotgen::Dynamic, O::value>::Random(5, 5);
   auto decomp = rotgen::svd(m);
 
-  do {
+  do
+  {
     rank = decomp.rank();
 
     auto u = decomp.U(rank);
@@ -56,7 +58,8 @@ TTS_CASE_TPL("SVD decomposition - Static case",
   auto m = rotgen::matrix<T, 5, 5, O::value>::Random();
   auto decomp = rotgen::svd(m);
 
-  do {
+  do
+  {
     rank = decomp.rank();
 
     auto u = decomp.U(rank);
diff --git a/test/unit/map/arithmetic_functions.cpp b/test/unit/map/arithmetic_functions.cpp
index 21346d9..7544014 100644
--- a/test/unit/map/arithmetic_functions.cpp
+++ b/test/unit/map/arithmetic_functions.cpp
@@ -5,10 +5,12 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
-#include "unit/common/arithmetic.hpp"
 #include <rotgen/rotgen.hpp>
 
+#include "unit/common/arithmetic.hpp"
+#include "unit/common/references.hpp"
+#include "unit/tests.hpp"
+
 TTS_CASE_TPL("Test dynamic map transposition-like operations",
              rotgen::tests::types)<typename T, typename O>(
   tts::type<tts::types<T, O>>)
diff --git a/test/unit/map/basic_api.cpp b/test/unit/map/basic_api.cpp
index a7bedd8..773f6ad 100644
--- a/test/unit/map/basic_api.cpp
+++ b/test/unit/map/basic_api.cpp
@@ -5,9 +5,10 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
 #include <rotgen/rotgen.hpp>
 
+#include "unit/tests.hpp"
+
 TTS_CASE_TPL("Function size", rotgen::tests::types)<typename T, typename O>(
   tts::type<tts::types<T, O>>)
 {
diff --git a/test/unit/map/constructors.cpp b/test/unit/map/constructors.cpp
index d034087..57d4f56 100644
--- a/test/unit/map/constructors.cpp
+++ b/test/unit/map/constructors.cpp
@@ -5,9 +5,10 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
 #include <rotgen/rotgen.hpp>
 
+#include "unit/tests.hpp"
+
 TTS_CASE_TPL("map constructor from pointer and single static size",
              rotgen::tests::types)<typename T, typename O>(
   tts::type<tts::types<T, O>>)
diff --git a/test/unit/map/cwise.cpp b/test/unit/map/cwise.cpp
index d4b5980..9dda215 100644
--- a/test/unit/map/cwise.cpp
+++ b/test/unit/map/cwise.cpp
@@ -5,10 +5,12 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
-#include "unit/common/cwise.hpp"
 #include <rotgen/rotgen.hpp>
 
+#include "unit/common/cwise.hpp"
+#include "unit/common/references.hpp"
+#include "unit/tests.hpp"
+
 TTS_CASE_TPL("Test dynamic map cwise operations",
              rotgen::tests::types)<typename T, typename O>(
   tts::type<tts::types<T, O>>)
diff --git a/test/unit/map/norms.cpp b/test/unit/map/norms.cpp
index 8802d33..2b93205 100644
--- a/test/unit/map/norms.cpp
+++ b/test/unit/map/norms.cpp
@@ -5,10 +5,12 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
-#include "unit/common/norms.hpp"
 #include <rotgen/rotgen.hpp>
 
+#include "unit/common/norms.hpp"
+#include "unit/common/references.hpp"
+#include "unit/tests.hpp"
+
 TTS_CASE_TPL("Test dynamic map norm operations",
              rotgen::tests::types)<typename T, typename O>(
   tts::type<tts::types<T, O>>)
diff --git a/test/unit/map/operators.cpp b/test/unit/map/operators.cpp
index fc65c21..ec0ed2a 100644
--- a/test/unit/map/operators.cpp
+++ b/test/unit/map/operators.cpp
@@ -5,9 +5,9 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
 #include <rotgen/rotgen.hpp>
-#include <vector>
+
+#include "unit/tests.hpp"
 
 template<typename MatrixType>
 void test_map_operations(rotgen::Index rows,
diff --git a/test/unit/matrix/arithmetic_functions.cpp b/test/unit/matrix/arithmetic_functions.cpp
index 75474fb..2930dd5 100644
--- a/test/unit/matrix/arithmetic_functions.cpp
+++ b/test/unit/matrix/arithmetic_functions.cpp
@@ -5,10 +5,13 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
-#include "unit/common/arithmetic.hpp"
+
 #include <rotgen/rotgen.hpp>
 
+#include "unit/common/arithmetic.hpp"
+#include "unit/common/references.hpp"
+#include "unit/tests.hpp"
+
 TTS_CASE_TPL("Test dynamic matrix transposition-like operations",
              rotgen::tests::types)<typename T, typename O>(
   tts::type<tts::types<T, O>>)
diff --git a/test/unit/matrix/basic_api.cpp b/test/unit/matrix/basic_api.cpp
index 1b8656a..e54a59b 100644
--- a/test/unit/matrix/basic_api.cpp
+++ b/test/unit/matrix/basic_api.cpp
@@ -5,9 +5,10 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
 #include <rotgen/rotgen.hpp>
 
+#include "unit/tests.hpp"
+
 TTS_CASE_TPL("Function size", rotgen::tests::types)<typename T, typename O>(
   tts::type<tts::types<T, O>>)
 {
diff --git a/test/unit/matrix/cwise.cpp b/test/unit/matrix/cwise.cpp
index dba7dde..ede5967 100644
--- a/test/unit/matrix/cwise.cpp
+++ b/test/unit/matrix/cwise.cpp
@@ -5,10 +5,12 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
-#include "unit/common/cwise.hpp"
 #include <rotgen/rotgen.hpp>
 
+#include "unit/common/cwise.hpp"
+#include "unit/common/references.hpp"
+#include "unit/tests.hpp"
+
 TTS_CASE_TPL("Test dynamic matrix cwise operations",
              rotgen::tests::types)<typename T, typename O>(
   tts::type<tts::types<T, O>>)
diff --git a/test/unit/matrix/generators.cpp b/test/unit/matrix/generators.cpp
index e03744f..3c12ca3 100644
--- a/test/unit/matrix/generators.cpp
+++ b/test/unit/matrix/generators.cpp
@@ -5,9 +5,10 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
 #include <rotgen/rotgen.hpp>
 
+#include "unit/tests.hpp"
+
 void test_value(auto const& matrix,
                 std::size_t rows,
                 std::size_t cols,
diff --git a/test/unit/matrix/inverse.cpp b/test/unit/matrix/inverse.cpp
index 83bcf34..dd9dacb 100644
--- a/test/unit/matrix/inverse.cpp
+++ b/test/unit/matrix/inverse.cpp
@@ -7,6 +7,7 @@
 //==================================================================================================
 #include <rotgen/rotgen.hpp>
 
+#include "unit/common/references.hpp"
 #include "unit/tests.hpp"
 
 TTS_CASE_TPL("Test dynamic matrix inverse",
diff --git a/test/unit/matrix/norms.cpp b/test/unit/matrix/norms.cpp
index da09570..cf9f812 100644
--- a/test/unit/matrix/norms.cpp
+++ b/test/unit/matrix/norms.cpp
@@ -5,10 +5,12 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
-#include "unit/common/norms.hpp"
 #include <rotgen/rotgen.hpp>
 
+#include "unit/common/norms.hpp"
+#include "unit/common/references.hpp"
+#include "unit/tests.hpp"
+
 TTS_CASE_TPL("Test dynamic matrix norm operations",
              rotgen::tests::types)<typename T, typename O>(
   tts::type<tts::types<T, O>>)
diff --git a/test/unit/matrix/operators.cpp b/test/unit/matrix/operators.cpp
index 2e43eba..e0b8c93 100644
--- a/test/unit/matrix/operators.cpp
+++ b/test/unit/matrix/operators.cpp
@@ -5,9 +5,10 @@
   SPDX-License-Identifier: BSL-1.0
 */
 //==================================================================================================
-#include "unit/tests.hpp"
 #include <rotgen/rotgen.hpp>
 
+#include "unit/tests.hpp"
+
 template<typename MatrixType>
 void test_matrix_operations(rotgen::Index rows,
                             rotgen::Index cols,
diff --git a/test/unit/quaternion/basic_api.cpp b/test/unit/quaternion/basic_api.cpp
new file mode 100644
index 0000000..e2f729b
--- /dev/null
+++ b/test/unit/quaternion/basic_api.cpp
@@ -0,0 +1,28 @@
+//==============================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==============================================================================
+
+#include <rotgen/rotgen.hpp>
+
+#include "unit/tests.hpp"
+
+TTS_CASE_TPL("Test coefficient accessors",
+             float,
+             double)<typename Scalar>(tts::type<Scalar>)
+{
+  rotgen::quaternion<Scalar> quaterion;
+
+  quaterion.w() = 1;
+  quaterion.x() = 2;
+  quaterion.y() = 3;
+  quaterion.z() = 4;
+
+  TTS_EQUAL(quaterion.w(), rotgen::Index{1});
+  TTS_EQUAL(quaterion.x(), rotgen::Index{2});
+  TTS_EQUAL(quaterion.y(), rotgen::Index{3});
+  TTS_EQUAL(quaterion.z(), rotgen::Index{4});
+};
diff --git a/test/unit/quaternion/constructors.cpp b/test/unit/quaternion/constructors.cpp
new file mode 100644
index 0000000..24200a0
--- /dev/null
+++ b/test/unit/quaternion/constructors.cpp
@@ -0,0 +1,32 @@
+//==================================================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==================================================================================================
+
+#include <rotgen/rotgen.hpp>
+
+#include "unit/tests.hpp"
+
+TTS_CASE_TPL("Default matrix dynamic constructor",
+             float,
+             double)<typename Scalar>(tts::type<Scalar>)
+{
+  rotgen::quaternion<Scalar> rotgen_quaterion;
+
+  TTS_PASS();
+};
+
+TTS_CASE_TPL("Default matrix dynamic constructor",
+             float,
+             double)<typename Scalar>(tts::type<Scalar>)
+{
+  rotgen::quaternion<Scalar> quaterion{1, 2, 3, 4};
+
+  TTS_EQUAL(quaterion.w(), rotgen::Index{1});
+  TTS_EQUAL(quaterion.x(), rotgen::Index{2});
+  TTS_EQUAL(quaterion.y(), rotgen::Index{3});
+  TTS_EQUAL(quaterion.z(), rotgen::Index{4});
+};
diff --git a/test/unit/quaternion/generators.cpp b/test/unit/quaternion/generators.cpp
new file mode 100644
index 0000000..6c439fe
--- /dev/null
+++ b/test/unit/quaternion/generators.cpp
@@ -0,0 +1,114 @@
+//==============================================================================
+/*
+  ROTGEN - Runtime Overlay for Eigen
+  Copyright : CODE RECKONS
+  SPDX-License-Identifier: BSL-1.0
+*/
+//==============================================================================
+
+#include <rotgen/rotgen.hpp>
+
+#include "unit/common/references.hpp"
+#include "unit/tests.hpp"
+#include <Eigen/Eigen>
+
+TTS_CASE_TPL("Test quaternion::angleAxis",
+             float,
+             double)<typename Scalar>(tts::type<Scalar>)
+{
+  // Eigen
+  Eigen::Quaternion<Scalar> eigen_result;
+  eigen_result =
+    Eigen::AngleAxis<Scalar>(13, Eigen::Matrix<Scalar, 3, 1>{1, 2, 3});
+
+  // Rotgen
+  rotgen::quaternion<Scalar> rotgen_result;
+  rotgen_result =
+    rotgen::angleAxis<Scalar>(13, rotgen::matrix<Scalar, 3, 1>{1, 2, 3});
+
+  TTS_EQUAL(rotgen_result.w(), eigen_result.w());
+  TTS_EQUAL(rotgen_result.x(), eigen_result.x());
+  TTS_EQUAL(rotgen_result.y(), eigen_result.y());
+  TTS_EQUAL(rotgen_result.z(), eigen_result.z());
+};
+
+TTS_CASE_TPL("Test quaternion::operator* with rotgen::matrix",
+             float,
+             double)<typename Scalar>(tts::type<Scalar>)
+{
+  // Eigen
+  Eigen::Quaternion<Scalar> eigen_quat;
+  eigen_quat =
+    Eigen::AngleAxis<Scalar>(13, Eigen::Matrix<Scalar, 3, 1>{1, 2, 3});
+
+  Eigen::Matrix<Scalar, 3, 1> eigen_other{3, 2, 1};
+  auto eigen_result = eigen_quat * eigen_other;
+
+  // Rotgen
+  rotgen::quaternion<Scalar> rotgen_quat;
+  rotgen_quat =
+    rotgen::angleAxis<Scalar>(13, rotgen::matrix<Scalar, 3, 1>{1, 2, 3});
+
+  rotgen::matrix<Scalar, 3, 1> rotgen_other{3, 2, 1};
+  auto rotgen_result = rotgen_quat * rotgen_other;
+
+  for (auto i = 0; i < eigen_result.size(); i++)
+    TTS_ULP_EQUAL(rotgen_result(i), eigen_result(i), 2);
+};
+
+TTS_CASE_TPL("Test quaternion::operator* with rotgen::map",
+             float,
+             double)<typename Scalar>(tts::type<Scalar>)
+{
+  Scalar data[] = {3, 2, 1};
+
+  // Eigen
+  Eigen::Quaternion<Scalar> eigen_quat;
+  eigen_quat =
+    Eigen::AngleAxis<Scalar>(13, Eigen::Matrix<Scalar, 3, 1>{1, 2, 3});
+
+  Eigen::Map<Eigen::Matrix<Scalar, 3, 1>> eigen_other{data};
+  auto eigen_result = eigen_quat * eigen_other;
+
+  // Rotgen
+  rotgen::quaternion<Scalar> rotgen_quat;
+  rotgen_quat =
+    rotgen::angleAxis<Scalar>(13, rotgen::matrix<Scalar, 3, 1>{1, 2, 3});
+
+  rotgen::map<rotgen::matrix<Scalar, 3, 1>> rotgen_other{data};
+  auto rotgen_result = rotgen_quat * rotgen_other;
+
+  for (auto i = 0; i < eigen_result.size(); i++)
+    TTS_ULP_EQUAL(rotgen_result(i), eigen_result(i), 2);
+};
+
+TTS_CASE_TPL("Test quaternion::operator* with rotgen::block",
+             float,
+             double)<typename Scalar>(tts::type<Scalar>)
+{
+
+  // Eigen
+  Eigen::Quaternion<Scalar> eigen_quat;
+  eigen_quat =
+    Eigen::AngleAxis<Scalar>(13, Eigen::Matrix<Scalar, 3, 1>{1, 2, 3});
+
+  Eigen::Matrix<Scalar, 16, 16> eigen_original_mat;
+  rotgen::tests::prepare(16, 16, rotgen::tests::default_init_function,
+                         eigen_original_mat);
+  auto eigen_other = eigen_original_mat.template block<3, 1>(0, 0);
+  auto eigen_result = eigen_quat * eigen_other;
+
+  // Rotgen
+  rotgen::quaternion<Scalar> rotgen_quat;
+  rotgen_quat =
+    rotgen::angleAxis<Scalar>(13, rotgen::matrix<Scalar, 3, 1>{1, 2, 3});
+
+  rotgen::matrix<Scalar, 16, 16> rotgen_original_mat;
+  rotgen::tests::prepare(16, 16, rotgen::tests::default_init_function,
+                         rotgen_original_mat);
+  auto rotgen_other = rotgen::extract<3, 1>(rotgen_original_mat, 0, 0);
+  auto rotgen_result = rotgen_quat * rotgen_other;
+
+  for (auto i = 0; i < eigen_result.size(); i++)
+    TTS_ULP_EQUAL(rotgen_result(i), eigen_result(i), 2);
+};
diff --git a/test/unit/tests.hpp b/test/unit/tests.hpp
index 4e45407..c3fbe3d 100644
--- a/test/unit/tests.hpp
+++ b/test/unit/tests.hpp
@@ -7,11 +7,10 @@
 //==================================================================================================
 #pragma once
 
-#include "tts.hpp"
-#include "unit/common/references.hpp"
-#include <rotgen/config.hpp>
 #include <rotgen/concepts.hpp>
-#include <functional>
+#include <rotgen/config.hpp>
+
+#include <tts.hpp>
 
 namespace rotgen::tests
 {