math: Use templates, fix alignment and 'iterator' (pointer) r/w access; enforce 'constexpr'; mat4f dropped Stack for constexpr; Dropped math.cpp having clear dependencies

Dropped math.cpp having clear dependencies

Adding ctor w/ initializer_list

Enforcing 'constexpr' where possible also enforces strict aliasing etc,
hence resulting in a more tight definition.

Template types for floating-point or integral types are used,
added the generic inner class typedef variations used in methods for clarity
and proper meta-programming.

alignas() on base type enforces utilization with arrays.

operator[] uses address of first element w/ strict alignment
and proper type (for constexpr).

14 files changed, 1630 insertions, 1429 deletions

diff --git a/JaulibSetup.cmake b/JaulibSetup.cmake
index 2848594..14c8ef2 100644
--- a/JaulibSetup.cmake
+++ b/JaulibSetup.cmake
@@ -80,6 +80,8 @@ message(STATUS "${PROJECT_NAME} USE_STRIP = ${USE_STRIP} (pre-set)")
 if(DEBUG)
     if(CMAKE_COMPILER_IS_GNUCC)
         set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -ggdb -DDEBUG -fno-omit-frame-pointer ${GCC_FLAGS_STACK} -no-pie")
+        set(CMAKE_SHARED_LINKER_FLAGS  "${CMAKE_SHARED_LINKER_FLAGS} -no-pie")
+        set(CMAKE_EXE_LINKER_FLAGS  "${CMAKE_EXE_LINKER_FLAGS} -no-pie")
     else()
         set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -DDEBUG")
     endif(CMAKE_COMPILER_IS_GNUCC)
@@ -96,33 +98,26 @@ if(DEBUG)
             set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${GCC_FLAGS_SANITIZE_THREAD}")
         endif(CMAKE_COMPILER_IS_GNUCC)
     endif(INSTRUMENTATION)
-elseif(GPROF)
-    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -g -ggdb -pg")
-elseif(PERF_ANALYSIS)
-    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -g -ggdb")
 else()
-    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3")
-    find_program(STRIP strip)
-    if (STRIP STREQUAL "STRIP-NOTFOUND")
-        set(USE_STRIP OFF)
-        message(STATUS "${PROJECT_NAME} USE_STRIP:=false, strip not found")
-    elseif(NOT DEFINED USE_STRIP)
-        set(USE_STRIP ON)
-        message(STATUS "${PROJECT_NAME} USE_STRIP:=true, !DEBUG and not set")
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -DNDEBUG")
+    if(GPROF)
+        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -ggdb -pg")
+    elseif(PERF_ANALYSIS)
+        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -ggdb")
+    else()
+        find_program(STRIP strip)
+        if (STRIP STREQUAL "STRIP-NOTFOUND")
+            set(USE_STRIP OFF)
+            message(STATUS "${PROJECT_NAME} USE_STRIP:=false, strip not found")
+        elseif(NOT DEFINED USE_STRIP)
+            set(USE_STRIP ON)
+            message(STATUS "${PROJECT_NAME} USE_STRIP:=true, !DEBUG and not set")
+        endif()
     endif()
 endif(DEBUG)
 
 message(STATUS "${PROJECT_NAME} USE_STRIP = ${USE_STRIP} (final)")
 
-if(DEBUG)
-    if(CMAKE_COMPILER_IS_GNUCC)
-        set(CMAKE_SHARED_LINKER_FLAGS  "${CMAKE_SHARED_LINKER_FLAGS} -no-pie")
-        set(CMAKE_EXE_LINKER_FLAGS  "${CMAKE_EXE_LINKER_FLAGS} -no-pie")
-    endif(CMAKE_COMPILER_IS_GNUCC)
-else()
-    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DNDEBUG")
-endif(DEBUG)
-
 if(DONT_USE_RTTI)
     message(STATUS "${PROJECT_NAME} RTTI disabled")
     set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fno-rtti")
diff --git a/include/jau/math/aabbox3f.hpp b/include/jau/math/aabbox3f.hpp
index aa0e1db..bf0ee90 100644
--- a/include/jau/math/aabbox3f.hpp
+++ b/include/jau/math/aabbox3f.hpp
@@ -24,12 +24,9 @@
 #ifndef JAU_AABBOX3F_HPP_
 #define JAU_AABBOX3F_HPP_
 
+#include <jau/functional.hpp>
 #include <jau/math/vec3f.hpp>
 
-namespace jau::math::geom::plane {
-    class AffineTransform; // forward
-}
-
 namespace jau::math {
 
     /** \addtogroup Math
@@ -260,13 +257,48 @@ namespace jau::math {
             }
 
             /**
+             * General purpose Vec3f transform function
+             */
+            typedef jau::function<jau::math::Vec3f(const jau::math::Vec3f&)> transform_vec3f_func;
+
+            /**
              * Resize the aabbox3f to encapsulate another AABox, which will be <i>transformed</i> on the fly first.
              * @param newBox aabbox3f to be encapsulated in
-             * @param t the {@link AffineTransform} applied on <i>newBox</i> on the fly
+             * @param transform the transform function, applied on <i>newBox</i> on the fly
              * @param tmpV3 temporary storage
              * @return this aabbox3f for chaining
              */
-            AABBox3f& resize(const AABBox3f& newBox, const geom::plane::AffineTransform& t, Vec3f& tmpV3) noexcept;
+            AABBox3f& resize(const AABBox3f& newBox, transform_vec3f_func& transform) noexcept {
+                /** test low */
+                {
+                    const jau::math::Vec3f newBL = transform(newBox.low());
+                    if (newBL.x < m_lo.x) {
+                        m_lo.x = newBL.x;
+                    }
+                    if (newBL.y < m_lo.y) {
+                        m_lo.y = newBL.y;
+                    }
+                    if (newBL.z < m_lo.z) {
+                        m_lo.z = newBL.z;
+                    }
+                }
+
+                /** test high */
+                {
+                    const jau::math::Vec3f newTR = transform(newBox.high());
+                    if (newTR.x > m_hi.x) {
+                        m_hi.x = newTR.x;
+                    }
+                    if (newTR.y > m_hi.y) {
+                        m_hi.y = newTR.y;
+                    }
+                    if (newTR.z > m_hi.z) {
+                        m_hi.z = newTR.z;
+                    }
+                }
+                computeCenter();
+                return *this;
+            }
 
             /**
              * Resize the aabbox3f to encapsulate the passed
diff --git a/include/jau/math/geom/plane/affine_transform.hpp b/include/jau/math/geom/plane/affine_transform.hpp
index 87fa485..60be2a0 100644
--- a/include/jau/math/geom/plane/affine_transform.hpp
+++ b/include/jau/math/geom/plane/affine_transform.hpp
@@ -507,6 +507,17 @@ namespace jau::math::geom::plane {
 
         /**
          * @param src source of transformation
+         * @return resulting Vec2f
+         */
+        Vec2f transform(const Vec2f& src) const noexcept {
+            const float x = src.x;
+            const float y = src.y;
+            return Vec2f( x * m00 + y * m01 + m02,
+                          x * m10 + y * m11 + m12 );
+        }
+
+        /**
+         * @param src source of transformation
          * @param dst destination of transformation, maybe be equal to <code>src</code>
          * @return dst for chaining
          */
@@ -519,6 +530,39 @@ namespace jau::math::geom::plane {
             return dst;
         }
 
+      private:
+        Vec3f transformVec3f(const Vec3f& src) const noexcept {
+            const float x = src.x;
+            const float y = src.y;
+            return Vec3f( x * m00 + y * m01 + m02,
+                          x * m10 + y * m11 + m12,
+                          src.z
+                        ); // just copy z
+        }
+
+      public:
+        /**
+         * @param src source of transformation
+         * @param dst destination of transformation, maybe be equal to <code>src</code>
+         * @return resulting Vec3f
+         */
+        Vec3f transform(const Vec3f& src) const noexcept {
+            return transformVec3f(src);
+        }
+
+        /**
+         * Resize the aabbox3f to encapsulate another AABox, which will be <i>transformed</i> on the fly first.
+         * @param newBox aabbox3f to be encapsulated in
+         * @param t the {@link AffineTransform} applied on <i>newBox</i> on the fly
+         * @param tmpV3 temporary storage
+         * @return this aabbox3f for chaining
+         */
+        AABBox3f& resizeBox(AABBox3f& box, const AABBox3f& newBox) noexcept {
+            AABBox3f::transform_vec3f_func transform(this, &AffineTransform::transformVec3f);
+            return box.resize(newBox, transform);
+
+        }
+
         /**
          *
          * @param src
diff --git a/include/jau/math/mat4f.hpp b/include/jau/math/mat4f.hpp
index e80e8d3..f3c51f1 100644
--- a/include/jau/math/mat4f.hpp
+++ b/include/jau/math/mat4f.hpp
@@ -27,9 +27,11 @@
 #include <cmath>
 #include <cstdarg>
 #include <cstdint>
+#include <cassert>
 #include <limits>
 #include <string>
 #include <vector>
+#include <initializer_list>
 #include <iostream>
 
 #include <jau/float_math.hpp>
@@ -50,10 +52,12 @@ namespace jau::math {
      *  @{
      */
 
+    template<typename Value_type,
+             std::enable_if_t<std::is_floating_point_v<Value_type>, bool>>
     class Quaternion; // forward
 
 /**
- * Basic 4x4 float matrix implementation using fields for intensive use-cases (host operations).
+ * Basic 4x4 value_type matrix implementation using fields for intensive use-cases (host operations).
  * <p>
  * Implementation covers {@link FloatUtil} matrix functionality, exposed in an object oriented manner.
  * </p>
@@ -102,111 +106,86 @@ namespace jau::math {
  * </ul>
  * </p>
  */
-class Mat4f {
-    private:
-        float m00, m10, m20, m30;
-        float m01, m11, m21, m31;
-        float m02, m12, m22, m32;
-        float m03, m13, m23, m33;
-
-        friend geom::Frustum;
-        friend Quaternion;
-
-        class Stack {
-            private:
-                int growSize;
-                std::vector<float> buffer;
-
-            public:
-                /**
-                 * @param initialSize initial size
-                 * @param growSize grow size if {@link #position()} is reached, maybe <code>0</code>
-                 */
-                Stack(int initialSize, int growSize_)
-                : growSize(growSize_), buffer(initialSize) {}
-
-                size_t growIfNecessary(int length) {
-                    const size_t p = buffer.size();
-                    const size_t nsz = buffer.size() + length;
-                    if( nsz > buffer.capacity() ) {
-                        buffer.reserve(buffer.size() + std::max(length, growSize));
-                    }
-                    buffer.resize(nsz);
-                    return p;
-                }
-
-                Mat4f& push(Mat4f& src) {
-                    size_t p = growIfNecessary(16);
-                    src.get(buffer, p);
-                    return src;
-                }
-
-                Mat4f& pop(Mat4f& dest) {
-                    size_t sz = buffer.size();
-                    if( sz < 16 ) {
-                        throw jau::IndexOutOfBoundsException(0, sz, E_FILE_LINE);
-                    }
-                    size_t p = sz - 16;
-                    dest.load(buffer, p);
-                    buffer.resize(p);
-                    return dest;
-                }
-        };
-
-        Stack stack; // start w/ zero size, growSize is half GL-min size (32)
+
+template<typename Value_type,
+         std::enable_if_t<std::is_floating_point_v<Value_type>, bool> = true>
+class alignas(Value_type) Matrix4 {
+  public:
+    typedef Value_type               value_type;
+    typedef value_type*              pointer;
+    typedef const value_type*        const_pointer;
+    typedef value_type&              reference;
+    typedef const value_type&        const_reference;
+    typedef value_type*              iterator;
+    typedef const value_type*        const_iterator;
+
+    typedef Vector4F<value_type, std::is_floating_point_v<Value_type>> Vec4;
+
+    constexpr static const value_type zero = value_type(0);
+    constexpr static const value_type one  = value_type(1);
+    constexpr static const value_type two  = value_type(2);
+    constexpr static const value_type half = one/two;
+
+  private:
+    //     RC
+    value_type m00, m10, m20, m30; // column 0
+    value_type m01, m11, m21, m31; // column 1
+    value_type m02, m12, m22, m32; // column 2
+    value_type m03, m13, m23, m33; // column 3
+
+    friend geom::Frustum;
+    friend Quaternion<value_type, std::is_floating_point_v<Value_type>>;
 
   public:
 
     /**
      * Creates a new identity matrix.
      */
-    Mat4f() noexcept
-    : m00(1.0f), m10(0.0f), m20(0.0f), m30(0.0f),
-      m01(0.0f), m11(1.0f), m21(0.0f), m31(0.0f),
-      m02(0.0f), m12(0.0f), m22(1.0f), m32(0.0f),
-      m03(0.0f), m13(0.0f), m23(0.0f), m33(1.0f),
-      stack(0, 16*16)
+    constexpr Matrix4() noexcept
+    : m00(one), m10(zero), m20(zero), m30(zero),
+      m01(zero), m11(one), m21(zero), m31(zero),
+      m02(zero), m12(zero), m22(one), m32(zero),
+      m03(zero), m13(zero), m23(zero), m33(one)
     { }
 
     /**
-     * Creates a new matrix based on given float[4*4] column major order.
+     * Creates a new matrix based on given value_type[4*4] column major order.
      * @param m 4x4 matrix in column-major order
      */
-    Mat4f(const float m[]) noexcept
-    : stack(0, 16*16)
-    {
-        load(m);
-    }
+    constexpr Matrix4(const_iterator m) noexcept
+    : m00(m[0+0*4]), m10(m[1+0*4]), m20(m[2+0*4]), m30(m[3+0*4]), // column 0
+      m01(m[0+1*4]), m11(m[1+1*4]), m21(m[2+1*4]), m31(m[3+1*4]), // column 1
+      m02(m[0+2*4]), m12(m[1+2*4]), m22(m[2+2*4]), m32(m[3+2*4]), // column 2
+      m03(m[0+3*4]), m13(m[1+3*4]), m23(m[2+3*4]), m33(m[3+3*4])  // column 3
+    {}
 
     /**
-     * Creates a new matrix based on given std::vector 4x4 column major order.
-     * @param m 4x4 matrix std::vector in column-major order starting at {@code src_off}
-     * @param m_off offset for matrix {@code src}
+     * Creates a new matrix based on given value_type initializer list in column major order.
+     * @param m source initializer list value_type data to be copied into this new instance, implied size must be >= 16
      */
-    Mat4f(const std::vector<float>& m, const size_t m_off) noexcept
-    : stack(0, 16*16)
+    constexpr Matrix4(std::initializer_list<value_type> m) noexcept
+    : m00(m.begin()[0+0*4]), m10(m.begin()[1+0*4]), m20(m.begin()[2+0*4]), m30(m.begin()[3+0*4]), // column 0
+      m01(m.begin()[0+1*4]), m11(m.begin()[1+1*4]), m21(m.begin()[2+1*4]), m31(m.begin()[3+1*4]), // column 1
+      m02(m.begin()[0+2*4]), m12(m.begin()[1+2*4]), m22(m.begin()[2+2*4]), m32(m.begin()[3+2*4]), // column 2
+      m03(m.begin()[0+3*4]), m13(m.begin()[1+3*4]), m23(m.begin()[2+3*4]), m33(m.begin()[3+3*4])  // column 3
     {
-        load(m, m_off);
+        assert(m.size() >= 16 );
     }
 
     /**
      * Creates a new matrix copying the values of the given {@code src} matrix.
-     *
-     * The stack is not copied.
      */
-    Mat4f(const Mat4f& o) noexcept
+    constexpr Matrix4(const Matrix4& o) noexcept
     : m00(o.m00), m10(o.m10), m20(o.m20), m30(o.m30),
       m01(o.m01), m11(o.m11), m21(o.m21), m31(o.m31),
       m02(o.m02), m12(o.m12), m22(o.m22), m32(o.m32),
-      m03(o.m03), m13(o.m13), m23(o.m23), m33(o.m33),
-      stack(0, 16*16) { }
+      m03(o.m03), m13(o.m13), m23(o.m23), m33(o.m33)
+    { }
 
     /**
      * Copy assignment using the the values of the given {@code src} matrix.
-     *
-     * The stack is not copied.
      */
-    Mat4f& operator=(const Mat4f& o) noexcept {
+    constexpr Matrix4& operator=(const Matrix4& o) noexcept {
         m00 = o.m00; m10 = o.m10; m20 = o.m20; m30 = o.m30;
         m01 = o.m01; m11 = o.m11; m21 = o.m21; m31 = o.m31;
         m02 = o.m02; m12 = o.m12; m22 = o.m22; m32 = o.m32;
@@ -214,7 +193,7 @@ class Mat4f {
         return *this;
     }
 
-    constexpr bool equals(const Mat4f& o, const float epsilon=std::numeric_limits<float>::epsilon()) const noexcept {
+    constexpr bool equals(const Matrix4& o, const value_type epsilon=std::numeric_limits<value_type>::epsilon()) const noexcept {
         if( this == &o ) {
             return true;
         } else {
@@ -236,7 +215,7 @@ class Mat4f {
                    jau::equals(m33, o.m33, epsilon);
         }
     }
-    constexpr bool operator==(const Mat4f& rhs) const noexcept {
+    constexpr bool operator==(const Matrix4& rhs) const noexcept {
         return equals(rhs);
     }
 
@@ -245,38 +224,22 @@ class Mat4f {
     //
 
     /**
-     * Returns writable reference to the {@code i}th component of the given column-major order matrix, 0 <= i < 16 w/o boundary check
+     * Returns writable reference to the {@code i}th component of this column-major order matrix, 0 <= i < 16 w/o boundary check
      */
-    float& operator[](size_t i) noexcept {
-        return reinterpret_cast<float*>(this)[i];
+    constexpr reference operator[](size_t i) noexcept {
+        assert( i < 16 );
+        return (&m00)[i];
     }
 
-    /** Sets the {@code i}th component with float {@code v} 0 <= i < 16 w/ boundary check*/
-    void set(const jau::nsize_t i, const float v) {
-        switch (i) {
-            case 0+4*0: m00 = v; break;
-            case 1+4*0: m10 = v; break;
-            case 2+4*0: m20 = v; break;
-            case 3+4*0: m30 = v; break;
-
-            case 0+4*1: m01 = v; break;
-            case 1+4*1: m11 = v; break;
-            case 2+4*1: m21 = v; break;
-            case 3+4*1: m31 = v; break;
-
-            case 0+4*2: m02 = v; break;
-            case 1+4*2: m12 = v; break;
-            case 2+4*2: m22 = v; break;
-            case 3+4*2: m32 = v; break;
-
-            case 0+4*3: m03 = v; break;
-            case 1+4*3: m13 = v; break;
-            case 2+4*3: m23 = v; break;
-            case 3+4*3: m33 = v; break;
-            default: throw jau::IndexOutOfBoundsException(i, 16, E_FILE_LINE);
-        }
+    /** Sets the {@code i}th component of this column-major order matrix with value_type {@code v}, 0 <= i < 16 w/o boundary check*/
+    constexpr void set(const jau::nsize_t i, const value_type v) noexcept {
+        assert( i < 16 );
+        (&m00)[i] = v;
     }
 
+    explicit operator pointer() noexcept { return &m00; }
+    constexpr iterator begin() noexcept { return &m00; }
+
     /**
      * Set this matrix to identity.
      * <pre>
@@ -288,12 +251,12 @@ class Mat4f {
      * </pre>
      * @return this matrix for chaining
      */
-    Mat4f& loadIdentity() noexcept {
-       m00 = m11 = m22 = m33 = 1.0f;
+    constexpr Matrix4& loadIdentity() noexcept {
+       m00 = m11 = m22 = m33 = one;
        m01 = m02 = m03 =
        m10 = m12 = m13 =
        m20 = m21 = m23 =
-       m30 = m31 = m32 = 0.0f;
+       m30 = m31 = m32 = zero;
        return *this;
     }
 
@@ -302,7 +265,7 @@ class Mat4f {
      * @param src the source values
      * @return this matrix for chaining
      */
-    Mat4f& load(const Mat4f& src) noexcept {
+    constexpr Matrix4& load(const Matrix4& src) noexcept {
         m00 = src.m00; m10 = src.m10; m20 = src.m20; m30 = src.m30;
         m01 = src.m01; m11 = src.m11; m21 = src.m21; m31 = src.m31;
         m02 = src.m02; m12 = src.m12; m22 = src.m22; m32 = src.m32;
@@ -312,10 +275,11 @@ class Mat4f {
 
     /**
      * Load the values of the given matrix {@code src} to this matrix w/o boundary check.
-     * @param src 4x4 matrix float[16] in column-major order
+     * @param src 4x4 matrix value_type[16] in column-major order
      * @return this matrix for chaining
      */
-    Mat4f& load(const float src[]) noexcept {
+    constexpr Matrix4& load(const_iterator src) noexcept {
+      // RC
         m00 = src[0+0*4]; // column 0
         m10 = src[1+0*4];
         m20 = src[2+0*4];
@@ -335,35 +299,6 @@ class Mat4f {
         return *this;
     }
 
-    /**
-     * Load the values of the given matrix {@code src} to this matrix w/ boundary check.
-     * @param src 4x4 matrix std::vector in column-major order starting at {@code src_off}
-     * @param src_off offset for matrix {@code src}
-     * @return this matrix for chaining
-     */
-    Mat4f& load(const std::vector<float>& src, const size_t src_off) {
-        if( src.size() < src_off+16 || src_off > std::numeric_limits<size_t>::max() - 15) {
-            throw jau::IndexOutOfBoundsException(src_off, 16, E_FILE_LINE);
-        }
-        m00 = src[src_off+0+0*4];
-        m10 = src[src_off+1+0*4];
-        m20 = src[src_off+2+0*4];
-        m30 = src[src_off+3+0*4];
-        m01 = src[src_off+0+1*4];
-        m11 = src[src_off+1+1*4];
-        m21 = src[src_off+2+1*4];
-        m31 = src[src_off+3+1*4];
-        m02 = src[src_off+0+2*4];
-        m12 = src[src_off+1+2*4];
-        m22 = src[src_off+2+2*4];
-        m32 = src[src_off+3+2*4];
-        m03 = src[src_off+0+3*4];
-        m13 = src[src_off+1+3*4];
-        m23 = src[src_off+2+3*4];
-        m33 = src[src_off+3+3*4];
-        return *this;
-    }
-
     //
     // Read out Matrix via get(..)
     //
@@ -371,82 +306,57 @@ class Mat4f {
     /**
      * Returns read-only {@code i}th component of the given column-major order matrix, 0 <= i < 16 w/o boundary check
      */
-    float operator[](size_t i) const noexcept {
-        return reinterpret_cast<const float*>(this)[i];
+    constexpr value_type operator[](size_t i) const noexcept {
+        assert( i < 16 );
+        return (&m00)[i];
     }
 
-    /** Returns the {@code i}th component of the given column-major order matrix, 0 <= i < 16, w/ boundary check */
-    float get(const jau::nsize_t i) const {
-        switch (i) {
-            case 0+4*0: return m00;
-            case 1+4*0: return m10;
-            case 2+4*0: return m20;
-            case 3+4*0: return m30;
-
-            case 0+4*1: return m01;
-            case 1+4*1: return m11;
-            case 2+4*1: return m21;
-            case 3+4*1: return m31;
-
-            case 0+4*2: return m02;
-            case 1+4*2: return m12;
-            case 2+4*2: return m22;
-            case 3+4*2: return m32;
-
-            case 0+4*3: return m03;
-            case 1+4*3: return m13;
-            case 2+4*3: return m23;
-            case 3+4*3: return m33;
-
-            default: throw jau::IndexOutOfBoundsException(i, 16, E_FILE_LINE);
-        }
+    /** Returns the {@code i}th component of the given column-major order matrix, 0 <= i < 16, w/o boundary check */
+    constexpr value_type get(const jau::nsize_t i) const noexcept {
+        assert( i < 16 );
+        return (&m00)[i];
     }
 
+    explicit operator const_pointer() const noexcept { return &m00; }
+    constexpr const_iterator cbegin() const noexcept { return &m00; }
+
     /**
-     * Get the named column of the given column-major matrix to v_out w/ boundary check.
+     * Get the named column of the given column-major matrix to v_out w/o boundary check.
      * @param column named column to copy
      * @param v_out the column-vector storage
      * @return given result vector <i>v_out</i> for chaining
      */
-    Vec4f& getColumn(const jau::nsize_t column, Vec4f& v_out) const {
-        if( column > 3 ) {
-            throw jau::IndexOutOfBoundsException(3+column*4, 16, E_FILE_LINE);
-        }
-        v_out.set( get(0+column*4),
-                   get(1+column*4),
-                   get(2+column*4),
-                   get(3+column*4) );
-        return v_out;
+    constexpr Vec4& getColumn(const jau::nsize_t column, Vec4& v_out) const noexcept {
+        assert( column < 4 );
+        return v_out.set( get(0+column*4),
+                          get(1+column*4),
+                          get(2+column*4),
+                          get(3+column*4) );
     }
+
     /**
-     * Get the named column of the given column-major matrix to v_out w/ boundary check.
+     * Get the named column of the given column-major matrix to v_out w/o boundary check.
      * @param column named column to copy
      * @return result vector holding the requested column
      */
-    Vec4f getColumn(const jau::nsize_t column) const {
-        if( column > 3 ) {
-            throw jau::IndexOutOfBoundsException(2+column*4, 16, E_FILE_LINE);
-        }
-        return Vec4f( get(0+column*4),
-                      get(1+column*4),
-                      get(2+column*4),
-                      get(3+column*4) );
+    constexpr Vec4 getColumn(const jau::nsize_t column) const noexcept {
+        assert( column < 4 );
+        return Vec4( get(0+column*4),
+                     get(1+column*4),
+                     get(2+column*4),
+                     get(3+column*4) );
     }
 
     /**
-     * Get the named column of the given column-major matrix to v_out w/ boundary check.
+     * Get the named column of the given column-major matrix to v_out w/o boundary check.
      * @param column named column to copy
      * @param v_out the column-vector storage
      * @return given result vector <i>v_out</i> for chaining
      */
-    Vec3f& getColumn(const jau::nsize_t column, Vec3f& v_out) const {
-        if( column > 3 ) {
-            throw jau::IndexOutOfBoundsException(2+column*4, 16, E_FILE_LINE);
-        }
-        v_out.set( get(0+column*4),
-                   get(1+column*4),
-                   get(2+column*4) );
-        return v_out;
+    constexpr Vec3f& getColumn(const jau::nsize_t column, Vec3f& v_out) const noexcept {
+        return v_out.set( get(0+column*4),
+                          get(1+column*4),
+                          get(2+column*4) );
     }
 
     /**
@@ -455,84 +365,44 @@ class Mat4f {
      * @param v_out the row-vector storage
      * @return given result vector <i>v_out</i> for chaining
      */
-    Vec4f& getRow(const jau::nsize_t row, Vec4f& v_out) const {
-        if( row > 3 ) {
-            throw jau::IndexOutOfBoundsException(row+3*4, 16, E_FILE_LINE);
-        }
-        v_out.set( get(row+0*4),
-                   get(row+1*4),
-                   get(row+2*4),
-                   get(row+3*4) );
-        return v_out;
+    constexpr Vec4& getRow(const jau::nsize_t row, Vec4& v_out) const noexcept {
+        return v_out.set( get(row+0*4),
+                          get(row+1*4),
+                          get(row+2*4),
+                          get(row+3*4) );
     }
     /**
-     * Get the named column of the given column-major matrix to v_out w/ boundary check.
+     * Get the named column of the given column-major matrix to v_out w/o boundary check.
      * @param row named row to copy
      * @return result vector holding the requested row
      */
-    Vec4f getRow(const jau::nsize_t row) const {
-        if( row > 3 ) {
-            throw jau::IndexOutOfBoundsException(row+3*4, 16, E_FILE_LINE);
-        }
-        return Vec4f( get(row+0*4),
-                      get(row+1*4),
-                      get(row+2*4),
-                      get(row+3*4) );
+    constexpr Vec4 getRow(const jau::nsize_t row) const noexcept {
+        return Vec4( get(row+0*4),
+                     get(row+1*4),
+                     get(row+2*4),
+                     get(row+3*4) );
     }
 
     /**
-     * Get the named row of the given column-major matrix to v_out w/ boundary check.
+     * Get the named row of the given column-major matrix to v_out w/o boundary check.
      * @param row named row to copy
-     * @param v_out the row-vector storage
+     * @param v_out the row-vector assert( i < 16 )e
      * @return given result vector <i>v_out</i> for chaining
      */
-    Vec3f& getRow(const jau::nsize_t row, Vec3f& v_out) const {
-        if( row > 3 ) {
-            throw jau::IndexOutOfBoundsException(row+2*4, 16, E_FILE_LINE);
-        }
-        v_out.set( get(row+0*4),
-                   get(row+1*4),
-                   get(row+2*4) );
-        return v_out;
+    constexpr Vec3f& getRow(const jau::nsize_t row, Vec3f& v_out) const noexcept {
+        assert( row <= 2 );
+        return v_out.set( get(row+0*4),
+                          get(row+1*4),
+                          get(row+2*4) );
     }
 
     /**
-     * Get this matrix into the given float[16] array at {@code dst_off} in column major order w/ boundary check.
+     * Get this matrix into the given value_type[16] array in column major order w/o boundary check.
      *
-     * @param dst float[16] array storage in column major order
-     * @param dst_off offset
+     * @param dst value_type[16] array storage in column major order
      * @return {@code dst} for chaining
      */
-    float* get(float* dst, size_t dst_off) const {
-        if( dst_off > std::numeric_limits<size_t>::max() - 15) {
-            throw jau::IndexOutOfBoundsException(dst_off, 16, E_FILE_LINE);
-        }
-        dst[dst_off++] = m00;
-        dst[dst_off++] = m10;
-        dst[dst_off++] = m20;
-        dst[dst_off++] = m30;
-        dst[dst_off++] = m01;
-        dst[dst_off++] = m11;
-        dst[dst_off++] = m21;
-        dst[dst_off++] = m31;
-        dst[dst_off++] = m02;
-        dst[dst_off++] = m12;
-        dst[dst_off++] = m22;
-        dst[dst_off++] = m32;
-        dst[dst_off++] = m03;
-        dst[dst_off++] = m13;
-        dst[dst_off++] = m23;
-        dst[dst_off++] = m33;
-        return dst;
-    }
-
-    /**
-     * Get this matrix into the given float[16] array in column major order w/o boundary check.
-     *
-     * @param dst float[16] array storage in column major order
-     * @return {@code dst} for chaining
-     */
-    float* get(float* dst) const noexcept {
+    constexpr iterator get(iterator dst) const noexcept {
         dst[0+0*4] = m00;
         dst[1+0*4] = m10;
         dst[2+0*4] = m20;
@@ -553,16 +423,14 @@ class Mat4f {
     }
 
     /**
-     * Get this matrix into the given {@link FloatBuffer} in column major order w/ boundary check.
+     * Get this matrix into the given {@link FloatBuffer} in column major order.
      *
      * @param dst 4x4 matrix std::vector in column-major order starting at {@code dst_off}
      * @param dst_off offset for matrix {@code dst}
      * @return {@code dst} for chaining
      */
-    std::vector<float>& get(std::vector<float>& dst, size_t dst_off) const {
-        if( dst.size() < dst_off+16 || dst_off > std::numeric_limits<size_t>::max() - 15) {
-            throw jau::IndexOutOfBoundsException(dst_off, 16, E_FILE_LINE);
-        }
+    constexpr std::vector<value_type>& get(std::vector<value_type>& dst, size_t dst_off) const noexcept {
+        assert( dst.size() >= dst_off+16 && dst_off <= std::numeric_limits<size_t>::max() - 15 );
         dst[dst_off++] = m00;
         dst[dst_off++] = m10;
         dst[dst_off++] = m20;
@@ -590,8 +458,8 @@ class Mat4f {
      * Returns the determinant of this matrix
      * @return the matrix determinant
      */
-    float determinant() const noexcept {
-        float ret = 0;
+    value_type determinant() const noexcept {
+        value_type ret = 0;
         ret += m00 * ( + m11*(m22*m33 - m23*m32) - m12*(m21*m33 - m23*m31) + m13*(m21*m32 - m22*m31));
         ret -= m01 * ( + m10*(m22*m33 - m23*m32) - m12*(m20*m33 - m23*m30) + m13*(m20*m32 - m22*m30));
         ret += m02 * ( + m10*(m21*m33 - m23*m31) - m11*(m20*m33 - m23*m30) + m13*(m20*m31 - m21*m30));
@@ -604,8 +472,8 @@ class Mat4f {
      *
      * @return this matrix for chaining
      */
-    Mat4f& transpose() noexcept {
-        float tmp;
+    Matrix4& transpose() noexcept {
+        value_type tmp;
 
         tmp = m10;
         m10 = m01;
@@ -640,7 +508,7 @@ class Mat4f {
      * @param src source 4x4 matrix
      * @return this matrix (result) for chaining
      */
-    Mat4f& transpose(const Mat4f& src) noexcept {
+    Matrix4& transpose(const Matrix4& src) noexcept {
         if( &src == this ) {
             return transpose();
         }
@@ -671,56 +539,56 @@ class Mat4f {
      * @return false if this matrix is singular and inversion not possible, otherwise true
      */
     bool invert() noexcept {
-        const float amax = absMax();
-        if( 0.0f == amax ) {
+        const value_type amax = absMax();
+        if( zero == amax ) {
             return false;
         }
-        const float scale = 1.0f/amax;
-        const float a00 = m00*scale;
-        const float a10 = m10*scale;
-        const float a20 = m20*scale;
-        const float a30 = m30*scale;
-
-        const float a01 = m01*scale;
-        const float a11 = m11*scale;
-        const float a21 = m21*scale;
-        const float a31 = m31*scale;
-
-        const float a02 = m02*scale;
-        const float a12 = m12*scale;
-        const float a22 = m22*scale;
-        const float a32 = m32*scale;
-
-        const float a03 = m03*scale;
-        const float a13 = m13*scale;
-        const float a23 = m23*scale;
-        const float a33 = m33*scale;
-
-        const float b00 = + a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31);
-        const float b01 = -( + a10*(a22*a33 - a23*a32) - a12*(a20*a33 - a23*a30) + a13*(a20*a32 - a22*a30));
-        const float b02 = + a10*(a21*a33 - a23*a31) - a11*(a20*a33 - a23*a30) + a13*(a20*a31 - a21*a30);
-        const float b03 = -( + a10*(a21*a32 - a22*a31) - a11*(a20*a32 - a22*a30) + a12*(a20*a31 - a21*a30));
-
-        const float b10 = -( + a01*(a22*a33 - a23*a32) - a02*(a21*a33 - a23*a31) + a03*(a21*a32 - a22*a31));
-        const float b11 = + a00*(a22*a33 - a23*a32) - a02*(a20*a33 - a23*a30) + a03*(a20*a32 - a22*a30);
-        const float b12 = -( + a00*(a21*a33 - a23*a31) - a01*(a20*a33 - a23*a30) + a03*(a20*a31 - a21*a30));
-        const float b13 = + a00*(a21*a32 - a22*a31) - a01*(a20*a32 - a22*a30) + a02*(a20*a31 - a21*a30);
-
-        const float b20 = + a01*(a12*a33 - a13*a32) - a02*(a11*a33 - a13*a31) + a03*(a11*a32 - a12*a31);
-        const float b21 = -( + a00*(a12*a33 - a13*a32) - a02*(a10*a33 - a13*a30) + a03*(a10*a32 - a12*a30));
-        const float b22 = + a00*(a11*a33 - a13*a31) - a01*(a10*a33 - a13*a30) + a03*(a10*a31 - a11*a30);
-        const float b23 = -( + a00*(a11*a32 - a12*a31) - a01*(a10*a32 - a12*a30) + a02*(a10*a31 - a11*a30));
-
-        const float b30 = -( + a01*(a12*a23 - a13*a22) - a02*(a11*a23 - a13*a21) + a03*(a11*a22 - a12*a21));
-        const float b31 = + a00*(a12*a23 - a13*a22) - a02*(a10*a23 - a13*a20) + a03*(a10*a22 - a12*a20);
-        const float b32 = -( + a00*(a11*a23 - a13*a21) - a01*(a10*a23 - a13*a20) + a03*(a10*a21 - a11*a20));
-        const float b33 = + a00*(a11*a22 - a12*a21) - a01*(a10*a22 - a12*a20) + a02*(a10*a21 - a11*a20);
-
-        const float det = (a00*b00 + a01*b01 + a02*b02 + a03*b03) / scale;
+        const value_type scale = one/amax;
+        const value_type a00 = m00*scale;
+        const value_type a10 = m10*scale;
+        const value_type a20 = m20*scale;
+        const value_type a30 = m30*scale;
+
+        const value_type a01 = m01*scale;
+        const value_type a11 = m11*scale;
+        const value_type a21 = m21*scale;
+        const value_type a31 = m31*scale;
+
+        const value_type a02 = m02*scale;
+        const value_type a12 = m12*scale;
+        const value_type a22 = m22*scale;
+        const value_type a32 = m32*scale;
+
+        const value_type a03 = m03*scale;
+        const value_type a13 = m13*scale;
+        const value_type a23 = m23*scale;
+        const value_type a33 = m33*scale;
+
+        const value_type b00 = + a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31);
+        const value_type b01 = -( + a10*(a22*a33 - a23*a32) - a12*(a20*a33 - a23*a30) + a13*(a20*a32 - a22*a30));
+        const value_type b02 = + a10*(a21*a33 - a23*a31) - a11*(a20*a33 - a23*a30) + a13*(a20*a31 - a21*a30);
+        const value_type b03 = -( + a10*(a21*a32 - a22*a31) - a11*(a20*a32 - a22*a30) + a12*(a20*a31 - a21*a30));
+
+        const value_type b10 = -( + a01*(a22*a33 - a23*a32) - a02*(a21*a33 - a23*a31) + a03*(a21*a32 - a22*a31));
+        const value_type b11 = + a00*(a22*a33 - a23*a32) - a02*(a20*a33 - a23*a30) + a03*(a20*a32 - a22*a30);
+        const value_type b12 = -( + a00*(a21*a33 - a23*a31) - a01*(a20*a33 - a23*a30) + a03*(a20*a31 - a21*a30));
+        const value_type b13 = + a00*(a21*a32 - a22*a31) - a01*(a20*a32 - a22*a30) + a02*(a20*a31 - a21*a30);
+
+        const value_type b20 = + a01*(a12*a33 - a13*a32) - a02*(a11*a33 - a13*a31) + a03*(a11*a32 - a12*a31);
+        const value_type b21 = -( + a00*(a12*a33 - a13*a32) - a02*(a10*a33 - a13*a30) + a03*(a10*a32 - a12*a30));
+        const value_type b22 = + a00*(a11*a33 - a13*a31) - a01*(a10*a33 - a13*a30) + a03*(a10*a31 - a11*a30);
+        const value_type b23 = -( + a00*(a11*a32 - a12*a31) - a01*(a10*a32 - a12*a30) + a02*(a10*a31 - a11*a30));
+
+        const value_type b30 = -( + a01*(a12*a23 - a13*a22) - a02*(a11*a23 - a13*a21) + a03*(a11*a22 - a12*a21));
+        const value_type b31 = + a00*(a12*a23 - a13*a22) - a02*(a10*a23 - a13*a20) + a03*(a10*a22 - a12*a20);
+        const value_type b32 = -( + a00*(a11*a23 - a13*a21) - a01*(a10*a23 - a13*a20) + a03*(a10*a21 - a11*a20));
+        const value_type b33 = + a00*(a11*a22 - a12*a21) - a01*(a10*a22 - a12*a20) + a02*(a10*a21 - a11*a20);
+
+        const value_type det = (a00*b00 + a01*b01 + a02*b02 + a03*b03) / scale;
         if( 0 == det ) {
             return false;
         }
-        const float invdet = 1.0f / det;
+        const value_type invdet = one / det;
 
         m00 = b00 * invdet;
         m10 = b01 * invdet;
@@ -749,58 +617,58 @@ class Mat4f {
      * @param src the source matrix, which values are to be inverted
      * @return false if {@code src} matrix is singular and inversion not possible, otherwise true
      */
-    bool invert(const Mat4f& src) noexcept {
-        const float amax = src.absMax();
-        if( 0.0f == amax ) {
+    bool invert(const Matrix4& src) noexcept {
+        const value_type amax = src.absMax();
+        if( zero == amax ) {
             return false;
         }
-        const float scale = 1.0f/amax;
-        const float a00 = src.m00*scale;
-        const float a10 = src.m10*scale;
-        const float a20 = src.m20*scale;
-        const float a30 = src.m30*scale;
-
-        const float a01 = src.m01*scale;
-        const float a11 = src.m11*scale;
-        const float a21 = src.m21*scale;
-        const float a31 = src.m31*scale;
-
-        const float a02 = src.m02*scale;
-        const float a12 = src.m12*scale;
-        const float a22 = src.m22*scale;
-        const float a32 = src.m32*scale;
-
-        const float a03 = src.m03*scale;
-        const float a13 = src.m13*scale;
-        const float a23 = src.m23*scale;
-        const float a33 = src.m33*scale;
-
-        const float b00 = + a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31);
-        const float b01 = -( + a10*(a22*a33 - a23*a32) - a12*(a20*a33 - a23*a30) + a13*(a20*a32 - a22*a30));
-        const float b02 = + a10*(a21*a33 - a23*a31) - a11*(a20*a33 - a23*a30) + a13*(a20*a31 - a21*a30);
-        const float b03 = -( + a10*(a21*a32 - a22*a31) - a11*(a20*a32 - a22*a30) + a12*(a20*a31 - a21*a30));
-
-        const float b10 = -( + a01*(a22*a33 - a23*a32) - a02*(a21*a33 - a23*a31) + a03*(a21*a32 - a22*a31));
-        const float b11 = + a00*(a22*a33 - a23*a32) - a02*(a20*a33 - a23*a30) + a03*(a20*a32 - a22*a30);
-        const float b12 = -( + a00*(a21*a33 - a23*a31) - a01*(a20*a33 - a23*a30) + a03*(a20*a31 - a21*a30));
-        const float b13 = + a00*(a21*a32 - a22*a31) - a01*(a20*a32 - a22*a30) + a02*(a20*a31 - a21*a30);
-
-        const float b20 = + a01*(a12*a33 - a13*a32) - a02*(a11*a33 - a13*a31) + a03*(a11*a32 - a12*a31);
-        const float b21 = -( + a00*(a12*a33 - a13*a32) - a02*(a10*a33 - a13*a30) + a03*(a10*a32 - a12*a30));
-        const float b22 = + a00*(a11*a33 - a13*a31) - a01*(a10*a33 - a13*a30) + a03*(a10*a31 - a11*a30);
-        const float b23 = -( + a00*(a11*a32 - a12*a31) - a01*(a10*a32 - a12*a30) + a02*(a10*a31 - a11*a30));
-
-        const float b30 = -( + a01*(a12*a23 - a13*a22) - a02*(a11*a23 - a13*a21) + a03*(a11*a22 - a12*a21));
-        const float b31 = + a00*(a12*a23 - a13*a22) - a02*(a10*a23 - a13*a20) + a03*(a10*a22 - a12*a20);
-        const float b32 = -( + a00*(a11*a23 - a13*a21) - a01*(a10*a23 - a13*a20) + a03*(a10*a21 - a11*a20));
-        const float b33 = + a00*(a11*a22 - a12*a21) - a01*(a10*a22 - a12*a20) + a02*(a10*a21 - a11*a20);
-
-        const float det = (a00*b00 + a01*b01 + a02*b02 + a03*b03) / scale;
+        const value_type scale = one/amax;
+        const value_type a00 = src.m00*scale;
+        const value_type a10 = src.m10*scale;
+        const value_type a20 = src.m20*scale;
+        const value_type a30 = src.m30*scale;
+
+        const value_type a01 = src.m01*scale;
+        const value_type a11 = src.m11*scale;
+        const value_type a21 = src.m21*scale;
+        const value_type a31 = src.m31*scale;
+
+        const value_type a02 = src.m02*scale;
+        const value_type a12 = src.m12*scale;
+        const value_type a22 = src.m22*scale;
+        const value_type a32 = src.m32*scale;
+
+        const value_type a03 = src.m03*scale;
+        const value_type a13 = src.m13*scale;
+        const value_type a23 = src.m23*scale;
+        const value_type a33 = src.m33*scale;
+
+        const value_type b00 = + a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31);
+        const value_type b01 = -( + a10*(a22*a33 - a23*a32) - a12*(a20*a33 - a23*a30) + a13*(a20*a32 - a22*a30));
+        const value_type b02 = + a10*(a21*a33 - a23*a31) - a11*(a20*a33 - a23*a30) + a13*(a20*a31 - a21*a30);
+        const value_type b03 = -( + a10*(a21*a32 - a22*a31) - a11*(a20*a32 - a22*a30) + a12*(a20*a31 - a21*a30));
+
+        const value_type b10 = -( + a01*(a22*a33 - a23*a32) - a02*(a21*a33 - a23*a31) + a03*(a21*a32 - a22*a31));
+        const value_type b11 = + a00*(a22*a33 - a23*a32) - a02*(a20*a33 - a23*a30) + a03*(a20*a32 - a22*a30);
+        const value_type b12 = -( + a00*(a21*a33 - a23*a31) - a01*(a20*a33 - a23*a30) + a03*(a20*a31 - a21*a30));
+        const value_type b13 = + a00*(a21*a32 - a22*a31) - a01*(a20*a32 - a22*a30) + a02*(a20*a31 - a21*a30);
+
+        const value_type b20 = + a01*(a12*a33 - a13*a32) - a02*(a11*a33 - a13*a31) + a03*(a11*a32 - a12*a31);
+        const value_type b21 = -( + a00*(a12*a33 - a13*a32) - a02*(a10*a33 - a13*a30) + a03*(a10*a32 - a12*a30));
+        const value_type b22 = + a00*(a11*a33 - a13*a31) - a01*(a10*a33 - a13*a30) + a03*(a10*a31 - a11*a30);
+        const value_type b23 = -( + a00*(a11*a32 - a12*a31) - a01*(a10*a32 - a12*a30) + a02*(a10*a31 - a11*a30));
+
+        const value_type b30 = -( + a01*(a12*a23 - a13*a22) - a02*(a11*a23 - a13*a21) + a03*(a11*a22 - a12*a21));
+        const value_type b31 = + a00*(a12*a23 - a13*a22) - a02*(a10*a23 - a13*a20) + a03*(a10*a22 - a12*a20);
+        const value_type b32 = -( + a00*(a11*a23 - a13*a21) - a01*(a10*a23 - a13*a20) + a03*(a10*a21 - a11*a20));
+        const value_type b33 = + a00*(a11*a22 - a12*a21) - a01*(a10*a22 - a12*a20) + a02*(a10*a21 - a11*a20);
+
+        const value_type det = (a00*b00 + a01*b01 + a02*b02 + a03*b03) / scale;
 
         if( 0 == det ) {
             return false;
         }
-        const float invdet = 1.0f / det;
+        const value_type invdet = one / det;
 
         m00 = b00 * invdet;
         m10 = b01 * invdet;
@@ -825,8 +693,8 @@ class Mat4f {
     }
 
   private:
-    float absMax() const noexcept {
-        float max = std::abs(m00);
+    value_type absMax() const noexcept {
+        value_type max = std::abs(m00);
         max = std::max(max, std::abs(m01));
         max = std::max(max, std::abs(m02));
         max = std::max(max, std::abs(m03));
@@ -850,17 +718,30 @@ class Mat4f {
 
   public:
     /**
+     * Multiply matrix with scalar: [this] = [this] x [s]
+     * @param s a scalar
+     * @return this matrix for chaining
+     */
+    constexpr Matrix4& operator*=( const value_type s ) noexcept {
+        m00 *= s; m10 *= s; m20 *= s; m30 *= s;
+        m01 *= s; m11 *= s; m21 *= s; m31 *= s;
+        m02 *= s; m12 *= s; m22 *= s; m32 *= s;
+        m03 *= s; m13 *= s; m23 *= s; m33 *= s;
+        return *this;
+    }
+
+    /**
      * Multiply matrix: [this] = [this] x [b]
      * @param b 4x4 matrix
      * @return this matrix for chaining
      * @see #mul(mat4f, mat4f)
      */
-    Mat4f& mul(const Mat4f& b) noexcept {
+    constexpr Matrix4& mul(const Matrix4& b) noexcept {
         // return mul(new mat4f(this), b); // <- roughly half speed
-        float ai0=m00; // row-0, m[0+0*4]
-        float ai1=m01;
-        float ai2=m02;
-        float ai3=m03;
+        value_type ai0=m00; // row-0, m[0+0*4]
+        value_type ai1=m01;
+        value_type ai2=m02;
+        value_type ai3=m03;
         m00 = ai0 * b.m00  +  ai1 * b.m10  +  ai2 * b.m20  +  ai3 * b.m30 ;
         m01 = ai0 * b.m01  +  ai1 * b.m11  +  ai2 * b.m21  +  ai3 * b.m31 ;
         m02 = ai0 * b.m02  +  ai1 * b.m12  +  ai2 * b.m22  +  ai3 * b.m32 ;
@@ -894,6 +775,15 @@ class Mat4f {
         m33 = ai0 * b.m03  +  ai1 * b.m13  +  ai2 * b.m23  +  ai3 * b.m33 ;
         return *this;
     }
+    /**
+     * Multiply matrix: [this] = [this] x [b]
+     * @param b 4x4 matrix
+     * @return this matrix for chaining
+     * @see #mul(mat4f, mat4f)
+     */
+    constexpr Matrix4& operator*=( const Matrix4& rhs ) noexcept {
+        return mul( rhs );
+    }
 
     /**
      * Multiply matrix: [this] = [a] x [b]
@@ -902,7 +792,7 @@ class Mat4f {
      * @return this matrix for chaining
      * @see #mul(mat4f)
      */
-    Mat4f& mul(const Mat4f& a, const Mat4f& b) noexcept {
+    constexpr Matrix4& mul(const Matrix4& a, const Matrix4& b) noexcept {
         // row-0, m[0+0*4]
         m00 = a.m00 * b.m00  +  a.m01 * b.m10  +  a.m02 * b.m20  +  a.m03 * b.m30 ;
         m01 = a.m00 * b.m01  +  a.m01 * b.m11  +  a.m02 * b.m21  +  a.m03 * b.m31 ;
@@ -935,23 +825,24 @@ class Mat4f {
      * @param v_out this * v_in
      * @returns v_out for chaining
      */
-    Vec4f& mulVec4f(const Vec4f& v_in, Vec4f& v_out) const noexcept {
+    constexpr Vec4& mulVec4(const Vec4& v_in, Vec4& v_out) const noexcept {
         // (one matrix row in column-major order) X (column vector)
-        const float x = v_in.x, y = v_in.y, z = v_in.z, w = v_in.w;
+        const value_type x = v_in.x, y = v_in.y, z = v_in.z, w = v_in.w;
         v_out.set( x * m00 + y * m01 + z * m02 + w * m03,
                    x * m10 + y * m11 + z * m12 + w * m13,
                    x * m20 + y * m21 + z * m22 + w * m23,
                    x * m30 + y * m31 + z * m32 + w * m33 );
         return v_out;
     }
+
     /**
-     * Returns new Vec4f holding this * v_in result
+     * Returns new Vec4, with this * v_in
      * @param v_in 4-component column-vector
      */
-    Vec4f operator*(const Vec4f& rhs) const noexcept {
+    constexpr Vec4 operator*(const Vec4& rhs) const noexcept {
         // (one matrix row in column-major order) X (column vector)
-        const float x = rhs.x, y = rhs.y, z = rhs.z, w = rhs.w;
-        return Vec4f( x * m00 + y * m01 + z * m02 + w * m03,
+        const value_type x = rhs.x, y = rhs.y, z = rhs.z, w = rhs.w;
+        return Vec4( x * m00 + y * m01 + z * m02 + w * m03,
                       x * m10 + y * m11 + z * m12 + w * m13,
                       x * m20 + y * m21 + z * m22 + w * m23,
                       x * m30 + y * m31 + z * m32 + w * m33 );
@@ -961,9 +852,9 @@ class Mat4f {
      * @param v_inout 4-component column-vector input and output, i.e. in-place transformation
      * @returns v_inout for chaining
      */
-    Vec4f& mulVec4f(Vec4f& v_inout) const noexcept {
+    constexpr Vec4& mulVec4(Vec4& v_inout) const noexcept {
         // (one matrix row in column-major order) X (column vector)
-        const float x = v_inout.x, y = v_inout.y, z = v_inout.z, w = v_inout.w;
+        const value_type x = v_inout.x, y = v_inout.y, z = v_inout.z, w = v_inout.w;
         v_inout.set( x * m00 + y * m01 + z * m02 + w * m03,
                      x * m10 + y * m11 + z * m12 + w * m13,
                      x * m20 + y * m21 + z * m22 + w * m23,
@@ -982,14 +873,30 @@ class Mat4f {
      * @param v_out m_in * v_in, 3-component column-vector {@link vec3f}
      * @returns v_out for chaining
      */
-    Vec3f& mulVec3f(const Vec3f& v_in, Vec3f& v_out) const noexcept {
+    constexpr Vec3f& mulVec3f(const Vec3f& v_in, Vec3f& v_out) const noexcept {
         // (one matrix row in column-major order) X (column vector)
-        const float x = v_in.x, y = v_in.y, z = v_in.z;
-        v_out.set( x * m00 + y * m01 + z * m02 + 1.0f * m03,
-                   x * m10 + y * m11 + z * m12 + 1.0f * m13,
-                   x * m20 + y * m21 + z * m22 + 1.0f * m23 );
+        const value_type x = v_in.x, y = v_in.y, z = v_in.z;
+        v_out.set( x * m00 + y * m01 + z * m02 + one * m03,
+                   x * m10 + y * m11 + z * m12 + one * m13,
+                   x * m20 + y * m21 + z * m22 + one * m23 );
         return v_out;
     }
+    /**
+     * Returns new Vec3f, with affine 3f-vector transformation by this 4x4 matrix: this * v_in
+     *
+     * 4x4 matrix multiplication with 3-component vector,
+     * using {@code 1} for for {@code v_in.w()} and dropping {@code v_out.w()},
+     * which shall be {@code 1}.
+     *
+     * @param v_in 3-component column-vector {@link vec3f}
+     */
+    constexpr Vec3f operator*(const Vec3f& rhs) const noexcept {
+        // (one matrix row in column-major order) X (column vector)
+        const value_type x = rhs.x, y = rhs.y, z = rhs.z;
+        return Vec3f( x * m00 + y * m01 + z * m02 + one * m03,
+                      x * m10 + y * m11 + z * m12 + one * m13,
+                      x * m20 + y * m21 + z * m22 + one * m23 );
+    }
 
     /**
      * Affine 3f-vector transformation by 4x4 matrix
@@ -1001,12 +908,12 @@ class Mat4f {
      * @param v_inout 3-component column-vector {@link vec3f} input and output, i.e. in-place transformation
      * @returns v_inout for chaining
      */
-    Vec3f& mulVec3f(Vec3f& v_inout) const noexcept {
+    constexpr Vec3f& mulVec3f(Vec3f& v_inout) const noexcept {
         // (one matrix row in column-major order) X (column vector)
-        const float x = v_inout.x, y = v_inout.y, z = v_inout.z;
-        v_inout.set( x * m00 + y * m01 + z * m02 + 1.0f * m03,
-                     x * m10 + y * m11 + z * m12 + 1.0f * m13,
-                     x * m20 + y * m21 + z * m22 + 1.0f * m23 );
+        const value_type x = v_inout.x, y = v_inout.y, z = v_inout.z;
+        v_inout.set( x * m00 + y * m01 + z * m02 + one * m03,
+                     x * m10 + y * m11 + z * m12 + one * m13,
+                     x * m20 + y * m21 + z * m22 + one * m23 );
         return v_inout;
     }
 
@@ -1028,15 +935,15 @@ class Mat4f {
      * @param z z-axis translate
      * @return this matrix for chaining
      */
-    Mat4f& setToTranslation(const float x, const float y, const float z) noexcept {
-        m00 = m11 = m22 = m33 = 1.0f;
+    Matrix4& setToTranslation(const value_type x, const value_type y, const value_type z) noexcept {
+        m00 = m11 = m22 = m33 = one;
         m03 = x;
         m13 = y;
         m23 = z;
         m01 = m02 =
         m10 = m12 =
         m20 = m21 =
-        m30 = m31 = m32 = 0.0f;
+        m30 = m31 = m32 = zero;
         return *this;
     }
 
@@ -1052,7 +959,7 @@ class Mat4f {
      * @param t translate vec3f
      * @return this matrix for chaining
      */
-    Mat4f& setToTranslation(const Vec3f& t) noexcept {
+    Matrix4& setToTranslation(const Vec3f& t) noexcept {
         return setToTranslation(t.x, t.y, t.z);
     }
 
@@ -1070,15 +977,15 @@ class Mat4f {
      * @param z z-axis scale
      * @return this matrix for chaining
      */
-    Mat4f& setToScale(const float x, const float y, const float z) noexcept {
-        m33 = 1.0f;
+    constexpr Matrix4& setToScale(const value_type x, const value_type y, const value_type z) noexcept {
+        m33 = one;
         m00 = x;
         m11 = y;
         m22 = z;
         m01 = m02 = m03 =
         m10 = m12 = m13 =
         m20 = m21 = m23 =
-        m30 = m31 = m32 = 0.0f;
+        m30 = m31 = m32 = zero;
         return *this;
     }
 
@@ -1094,7 +1001,7 @@ class Mat4f {
      * @param s scale vec3f
      * @return this matrix for chaining
      */
-    Mat4f& setToScale(const Vec3f& s) noexcept {
+    constexpr Matrix4& setToScale(const Vec3f& s) noexcept {
         return setToScale(s.x, s.y, s.z);
     }
 
@@ -1114,39 +1021,39 @@ class Mat4f {
      * @param z z of rotation axis
      * @return this matrix for chaining
      */
-    Mat4f& setToRotationAxis(const float ang_rad, float x, float y, float z) noexcept {
-        const float c = std::cos(ang_rad);
-        const float ic= 1.0f - c;
-        const float s = std::sin(ang_rad);
+    Matrix4& setToRotationAxis(const value_type ang_rad, value_type x, value_type y, value_type z) noexcept {
+        const value_type c = std::cos(ang_rad);
+        const value_type ic= one - c;
+        const value_type s = std::sin(ang_rad);
 
         Vec3f tmp(x, y, z); tmp.normalize();
         x = tmp.x; y = tmp.y; z = tmp.z;
 
-        const float xy = x*y;
-        const float xz = x*z;
-        const float xs = x*s;
-        const float ys = y*s;
-        const float yz = y*z;
-        const float zs = z*s;
+        const value_type xy = x*y;
+        const value_type xz = x*z;
+        const value_type xs = x*s;
+        const value_type ys = y*s;
+        const value_type yz = y*z;
+        const value_type zs = z*s;
         m00 = x*x*ic+c;
         m10 = xy*ic+zs;
         m20 = xz*ic-ys;
-        m30 = 0.0f;
+        m30 = zero;
 
         m01 = xy*ic-zs;
         m11 = y*y*ic+c;
         m21 = yz*ic+xs;
-        m31 = 0.0f;
+        m31 = zero;
 
         m02 = xz*ic+ys;
         m12 = yz*ic-xs;
         m22 = z*z*ic+c;
-        m32 = 0.0f;
+        m32 = zero;
 
         m03 = 0.9f;
-        m13 = 0.0f;
-        m23 = 0.0f;
-        m33 = 1.0f;
+        m13 = zero;
+        m23 = zero;
+        m33 = one;
 
         return *this;
     }
@@ -1165,7 +1072,7 @@ class Mat4f {
      * @param axis rotation axis
      * @return this matrix for chaining
      */
-    Mat4f& setToRotationAxis(const float ang_rad, const Vec3f& axis) noexcept {
+    Matrix4& setToRotationAxis(const value_type ang_rad, const Vec3f& axis) noexcept {
         return setToRotationAxis(ang_rad, axis.x, axis.y, axis.z);
     }
 
@@ -1192,34 +1099,34 @@ class Mat4f {
      * @see <a href="http://www.euclideanspace.com/maths/geometry/rotations/conversions/eulerToMatrix/index.htm">euclideanspace.com-eulerToMatrix</a>
      * @see #setToRotation(Quaternion)
      */
-    Mat4f& setToRotationEuler(const float bankX, const float headingY, const float attitudeZ) noexcept {
+    Matrix4& setToRotationEuler(const value_type bankX, const value_type headingY, const value_type attitudeZ) noexcept {
         // Assuming the angles are in radians.
-        const float ch = std::cos(headingY);
-        const float sh = std::sin(headingY);
-        const float ca = std::cos(attitudeZ);
-        const float sa = std::sin(attitudeZ);
-        const float cb = std::cos(bankX);
-        const float sb = std::sin(bankX);
+        const value_type ch = std::cos(headingY);
+        const value_type sh = std::sin(headingY);
+        const value_type ca = std::cos(attitudeZ);
+        const value_type sa = std::sin(attitudeZ);
+        const value_type cb = std::cos(bankX);
+        const value_type sb = std::sin(bankX);
 
         m00 =  ch*ca;
         m10 =  sa;
         m20 = -sh*ca;
-        m30 =  0.0f;
+        m30 =  zero;
 
         m01 =  sh*sb    - ch*sa*cb;
         m11 =  ca*cb;
         m21 =  sh*sa*cb + ch*sb;
-        m31 =  0.0f;
+        m31 =  zero;
 
         m02 =  ch*sa*sb + sh*cb;
         m12 = -ca*sb;
         m22 = -sh*sa*sb + ch*cb;
-        m32 =  0.0f;
+        m32 =  zero;
 
-        m03 =  0.0f;
-        m13 =  0.0f;
-        m23 =  0.0f;
-        m33 =  1.0f;
+        m03 =  zero;
+        m13 =  zero;
+        m23 =  zero;
+        m33 =  one;
 
         return *this;
     }
@@ -1245,38 +1152,11 @@ class Mat4f {
      * @see <a href="http://www.euclideanspace.com/maths/geometry/rotations/conversions/eulerToMatrix/index.htm">euclideanspace.com-eulerToMatrix</a>
      * @see #setToRotation(Quaternion)
      */
-    Mat4f& setToRotationEuler(const Vec3f& angradXYZ) noexcept {
+    Matrix4& setToRotationEuler(const Vec3f& angradXYZ) noexcept {
         return setToRotationEuler(angradXYZ.x, angradXYZ.y, angradXYZ.z);
     }
 
     /**
-     * Set this matrix to rotation using the given Quaternion.
-     * <p>
-     * Implementation Details:
-     * <ul>
-     *   <li> makes identity matrix if {@link #magnitudeSquared()} is {@link FloatUtil#isZero(float, float) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
-     *   <li> The fields [m00 .. m22] define the rotation</li>
-     * </ul>
-     * </p>
-     *
-     * @param q the Quaternion representing the rotation
-     * @return this matrix for chaining
-     * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q54">Matrix-FAQ Q54</a>
-     * @see Quaternion#toMatrix(float[])
-     * @see #getRotation()
-     */
-    Mat4f& setToRotation(const Quaternion& q) noexcept;
-
-    /**
-     * Returns the rotation [m00 .. m22] fields converted to a Quaternion.
-     * @param res resulting Quaternion
-     * @return the resulting Quaternion for chaining.
-     * @see Quaternion#setFromMatrix(float, float, float, float, float, float, float, float, float)
-     * @see #setToRotation(Quaternion)
-     */
-    Quaternion& getRotation(Quaternion& res) const noexcept;
-
-    /**
      * Set this matrix to orthogonal projection.
      * <pre>
       Ortho matrix (Column Order):
@@ -1293,31 +1173,31 @@ class Mat4f {
      * @param zFar
      * @return this matrix for chaining
      */
-    Mat4f& setToOrtho(const float left, const float right,
-                      const float bottom, const float top,
-                      const float zNear, const float zFar) noexcept {
+    Matrix4& setToOrtho(const value_type left, const value_type right,
+                      const value_type bottom, const value_type top,
+                      const value_type zNear, const value_type zFar) noexcept {
         {
-            // m00 = m11 = m22 = m33 = 1.0f;
-            m10 = m20 = m30 = 0.0f;
-            m01 = m21 = m31 = 0.0f;
-            m02 = m12 = m32 = 0.0f;
-            // m03 = m13 = m23 = 0.0f;
+            // m00 = m11 = m22 = m33 = one;
+            m10 = m20 = m30 = zero;
+            m01 = m21 = m31 = zero;
+            m02 = m12 = m32 = zero;
+            // m03 = m13 = m23 = zero;
         }
-        const float dx=right-left;
-        const float dy=top-bottom;
-        const float dz=zFar-zNear;
-        const float tx=-1.0f*(right+left)/dx;
-        const float ty=-1.0f*(top+bottom)/dy;
-        const float tz=-1.0f*(zFar+zNear)/dz;
+        const value_type dx=right-left;
+        const value_type dy=top-bottom;
+        const value_type dz=zFar-zNear;
+        const value_type tx=-one*(right+left)/dx;
+        const value_type ty=-one*(top+bottom)/dy;
+        const value_type tz=-one*(zFar+zNear)/dz;
 
-        m00 =  2.0f/dx;
-        m11 =  2.0f/dy;
-        m22 = -2.0f/dz;
+        m00 =  two/dx;
+        m11 =  two/dy;
+        m22 = -two/dz;
 
         m03 = tx;
         m13 = ty;
         m23 = tz;
-        m33 = 1.0f;
+        m33 = one;
 
         return *this;
     }
@@ -1341,10 +1221,10 @@ class Mat4f {
      * @throws IllegalArgumentException if {@code zNear <= 0} or {@code zFar <= zNear}
      *                                  or {@code left == right}, or {@code bottom == top}.
      */
-    Mat4f& setToFrustum(const float left, const float right,
-                        const float bottom, const float top,
-                        const float zNear, const float zFar) {
-        if( zNear <= 0.0f || zFar <= zNear ) {
+    Matrix4& setToFrustum(const value_type left, const value_type right,
+                        const value_type bottom, const value_type top,
+                        const value_type zNear, const value_type zFar) {
+        if( zNear <= zero || zFar <= zNear ) {
             throw jau::IllegalArgumentException("Requirements zNear > 0 and zFar > zNear, but zNear "+std::to_string(zNear)+", zFar "+std::to_string(zFar), E_FILE_LINE);
         }
         if( left == right || top == bottom) {
@@ -1352,18 +1232,18 @@ class Mat4f {
         }
         {
             // m00 = m11 = m22 = m33 = 1f;
-            m10 = m20 = m30 = 0.0f;
-            m01 = m21 = m31 = 0.0f;
-            m03 = m13 = 0.0f;
+            m10 = m20 = m30 = zero;
+            m01 = m21 = m31 = zero;
+            m03 = m13 = zero;
         }
-        const float zNear2 = 2.0f*zNear;
-        const float dx=right-left;
-        const float dy=top-bottom;
-        const float dz=zFar-zNear;
-        const float A=(right+left)/dx;
-        const float B=(top+bottom)/dy;
-        const float C=-1.0f*(zFar+zNear)/dz;
-        const float D=-2.0f*(zFar*zNear)/dz;
+        const value_type zNear2 = two*zNear;
+        const value_type dx=right-left;
+        const value_type dy=top-bottom;
+        const value_type dz=zFar-zNear;
+        const value_type A=(right+left)/dx;
+        const value_type B=(top+bottom)/dy;
+        const value_type C=-one*(zFar+zNear)/dz;
+        const value_type D=-two*(zFar*zNear)/dz;
 
         m00 = zNear2/dx;
         m11 = zNear2/dy;
@@ -1371,16 +1251,16 @@ class Mat4f {
         m02 = A;
         m12 = B;
         m22 = C;
-        m32 = -1.0f;
+        m32 = -one;
 
         m23 = D;
-        m33 = 0.0f;
+        m33 = zero;
 
         return *this;
     }
 
     /**
-     * Set this matrix to perspective {@link #setToFrustum(float, float, float, float, float, float) frustum} projection.
+     * Set this matrix to perspective {@link #setToFrustum(value_type, value_type, value_type, value_type, value_type, value_type) frustum} projection.
      *
      * @param fovy_rad angle in radians
      * @param aspect aspect ratio width / height
@@ -1388,50 +1268,37 @@ class Mat4f {
      * @param zFar
      * @return this matrix for chaining
      * @throws IllegalArgumentException if {@code zNear <= 0} or {@code zFar <= zNear}
-     * @see #setToFrustum(float, float, float, float, float, float)
+     * @see #setToFrustum(value_type, value_type, value_type, value_type, value_type, value_type)
      */
-    Mat4f& setToPerspective(const float fovy_rad, const float aspect, const float zNear, const float zFar) {
-        const float top    =  std::tan(fovy_rad/2.0f) * zNear; // use tangent of half-fov !
-        const float bottom =  -1.0f * top;    //          -1f * fovhvTan.top * zNear
-        const float left   = aspect * bottom; // aspect * -1f * fovhvTan.top * zNear
-        const float right  = aspect * top;    // aspect * fovhvTan.top * zNear
+    Matrix4& setToPerspective(const value_type fovy_rad, const value_type aspect, const value_type zNear, const value_type zFar) {
+        const value_type top    =  std::tan(fovy_rad/two) * zNear; // use tangent of half-fov !
+        const value_type bottom =  -one * top;    //          -1f * fovhvTan.top * zNear
+        const value_type left   = aspect * bottom; // aspect * -1f * fovhvTan.top * zNear
+        const value_type right  = aspect * top;    // aspect * fovhvTan.top * zNear
         return setToFrustum(left, right, bottom, top, zNear, zFar);
     }
 
     /**
-     * Set this matrix to perspective {@link #setToFrustum(float, float, float, float, float, float) frustum} projection.
+     * Set this matrix to perspective {@link #setToFrustum(value_type, value_type, value_type, value_type, value_type, value_type) frustum} projection.
      *
      * @param fovhv {@link FovHVHalves} field of view in both directions, may not be centered, either in radians or tangent
      * @param zNear
      * @param zFar
      * @return this matrix for chaining
      * @throws IllegalArgumentException if {@code zNear <= 0} or {@code zFar <= zNear}
-     * @see #setToFrustum(float, float, float, float, float, float)
+     * @see #setToFrustum(value_type, value_type, value_type, value_type, value_type, value_type)
      * @see Frustum#updateByFovDesc(mat4f, com.jogamp.math.geom.Frustum.FovDesc)
      */
-    Mat4f& setToPerspective(const FovHVHalves& fovhv, const float zNear, const float zFar) {
+    Matrix4& setToPerspective(const FovHVHalves& fovhv, const value_type zNear, const value_type zFar) {
         const FovHVHalves fovhvTan = fovhv.toTangents();  // use tangent of half-fov !
-        const float top    =         fovhvTan.top    * zNear;
-        const float bottom = -1.0f * fovhvTan.bottom * zNear;
-        const float left   = -1.0f * fovhvTan.left   * zNear;
-        const float right  =         fovhvTan.right  * zNear;
+        const value_type top    =         fovhvTan.top    * zNear;
+        const value_type bottom = -one * fovhvTan.bottom * zNear;
+        const value_type left   = -one * fovhvTan.left   * zNear;
+        const value_type right  =         fovhvTan.right  * zNear;
         return setToFrustum(left, right, bottom, top, zNear, zFar);
     }
 
     /**
-     * Calculate the frustum planes in world coordinates
-     * using this column major order matrix, usually a projection (P) or premultiplied P*MV matrix.
-     *
-     * Frustum plane's normals will point to the inside of the viewing frustum,
-     * as required by the {@link Frustum} class.
-     *
-     * May use geom::Frustum::setFromMat4() directly.
-     *
-     * @see geom::Frustum::setFromMat4()
-     */
-    geom::Frustum& getFrustum(geom::Frustum& frustum) noexcept;
-
-    /**
      * Set this matrix to the <i>look-at</i> matrix based on given parameters.
      * <p>
      * Consist out of two matrix multiplications:
@@ -1451,7 +1318,7 @@ class Mat4f {
      * @param tmp temporary mat4f used for multiplication
      * @return this matrix for chaining
      */
-    Mat4f& setToLookAt(const Vec3f& eye, const Vec3f& center, const Vec3f& up, Mat4f& tmp) noexcept {
+    Matrix4& setToLookAt(const Vec3f& eye, const Vec3f& center, const Vec3f& up, Matrix4& tmp) noexcept {
         // normalized forward!
         const Vec3f fwd = ( center - eye ).normalize();
 
@@ -1503,8 +1370,8 @@ class Mat4f {
      * </p>
      * <p>
      * To effectively use the generated pick matrix for picking,
-     * call {@link #setToPick(float, float, float, float, Recti, mat4f) setToPick(..)}
-     * and multiply a {@link #setToPerspective(float, float, float, float) custom perspective matrix}
+     * call {@link #setToPick(value_type, value_type, value_type, value_type, Recti, mat4f) setToPick(..)}
+     * and multiply a {@link #setToPerspective(value_type, value_type, value_type, value_type) custom perspective matrix}
      * by this pick matrix. Then you may load the result onto the perspective matrix stack.
      * </p>
      * @param x the center x-component of a picking region in window coordinates
@@ -1515,16 +1382,16 @@ class Mat4f {
      * @param mat4Tmp temp storage
      * @return true if successful or false if either delta value is <= zero.
      */
-    bool setToPick(const float x, const float y, const float deltaX, const float deltaY,
-                   const Recti& viewport, Mat4f& mat4Tmp) noexcept {
+    bool setToPick(const value_type x, const value_type y, const value_type deltaX, const value_type deltaY,
+                   const Recti& viewport, Matrix4& mat4Tmp) noexcept {
         if (deltaX <= 0 || deltaY <= 0) {
             return false;
         }
         /* Translate and scale the picked region to the entire window */
-        setToTranslation( ( viewport.width()  - 2.0f * ( x - viewport.x() ) ) / deltaX,
-                          ( viewport.height() - 2.0f * ( y - viewport.y() ) ) / deltaY,
+        setToTranslation( ( viewport.width()  - two * ( x - viewport.x() ) ) / deltaX,
+                          ( viewport.height() - two * ( y - viewport.y() ) ) / deltaY,
                           0);
-        mat4Tmp.setToScale( viewport.width() / deltaX, viewport.height() / deltaY, 1.0f );
+        mat4Tmp.setToScale( viewport.width() / deltaX, viewport.height() / deltaY, one );
         mul(mat4Tmp);
         return true;
     }
@@ -1534,7 +1401,7 @@ class Mat4f {
     //
 
     /**
-     * Rotate this matrix about give axis and angle in radians, i.e. multiply by {@link #setToRotationAxis(float, float, float, float) axis-rotation matrix}.
+     * Rotate this matrix about give axis and angle in radians, i.e. multiply by {@link #setToRotationAxis(value_type, value_type, value_type, value_type) axis-rotation matrix}.
      * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q38">Matrix-FAQ Q38</a>
      * @param angrad angle in radians
      * @param x x of rotation axis
@@ -1543,40 +1410,31 @@ class Mat4f {
      * @param tmp temporary mat4f used for multiplication
      * @return this matrix for chaining
      */
-    Mat4f& rotate(const float ang_rad, const float x, const float y, const float z, Mat4f& tmp) noexcept {
+    Matrix4& rotate(const value_type ang_rad, const value_type x, const value_type y, const value_type z, Matrix4& tmp) noexcept {
         return mul( tmp.setToRotationAxis(ang_rad, x, y, z) );
     }
 
     /**
-     * Rotate this matrix about give axis and angle in radians, i.e. multiply by {@link #setToRotationAxis(float, vec3f) axis-rotation matrix}.
+     * Rotate this matrix about give axis and angle in radians, i.e. multiply by {@link #setToRotationAxis(value_type, vec3f) axis-rotation matrix}.
      * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q38">Matrix-FAQ Q38</a>
      * @param angrad angle in radians
      * @param axis rotation axis
      * @param tmp temporary mat4f used for multiplication
      * @return this matrix for chaining
      */
-    Mat4f& rotate(const float ang_rad, const Vec3f& axis, Mat4f& tmp) noexcept {
+    Matrix4& rotate(const value_type ang_rad, const Vec3f& axis, Matrix4& tmp) noexcept {
         return mul( tmp.setToRotationAxis(ang_rad, axis) );
     }
 
     /**
-     * Rotate this matrix with the given {@link Quaternion}, i.e. multiply by {@link #setToRotation(Quaternion) Quaternion's rotation matrix}.
-     * @param tmp temporary mat4f used for multiplication
-     * @return this matrix for chaining
-     */
-    Mat4f& rotate(const Quaternion& quat, Mat4f& tmp) noexcept {
-        return mul( tmp.setToRotation(quat) );
-    }
-
-    /**
-     * Translate this matrix, i.e. multiply by {@link #setToTranslation(float, float, float) translation matrix}.
+     * Translate this matrix, i.e. multiply by {@link #setToTranslation(value_type, value_type, value_type) translation matrix}.
      * @param x x translation
      * @param y y translation
      * @param z z translation
      * @param tmp temporary mat4f used for multiplication
      * @return this matrix for chaining
      */
-    Mat4f& translate(const float x, const float y, const float z, Mat4f& tmp) noexcept {
+    Matrix4& translate(const value_type x, const value_type y, const value_type z, Matrix4& tmp) noexcept {
         return mul( tmp.setToTranslation(x, y, z) );
     }
 
@@ -1587,29 +1445,29 @@ class Mat4f {
      * @param tmp temporary mat4f used for multiplication
      * @return this matrix for chaining
      */
-    public final Mat4f translate(final Vec3f t, final Mat4f tmp) {
+    public final Matrix4 translate(final Vec3f t, final Matrix4 tmp) {
         return mul( tmp.setToTranslation(t) );
     }
 
     /**
-     * Scale this matrix, i.e. multiply by {@link #setToScale(float, float, float) scale matrix}.
+     * Scale this matrix, i.e. multiply by {@link #setToScale(value_type, value_type, value_type) scale matrix}.
      * @param x x scale
      * @param y y scale
      * @param z z scale
      * @param tmp temporary mat4f used for multiplication
      * @return this matrix for chaining
      */
-    public final Mat4f scale(const float x, const float y, const float z, final Mat4f tmp) {
+    public final Matrix4 scale(const value_type x, const value_type y, const value_type z, final Matrix4 tmp) {
         return mul( tmp.setToScale(x, y, z) );
     }
 
     /**
-     * Scale this matrix, i.e. multiply by {@link #setToScale(float, float, float) scale matrix}.
+     * Scale this matrix, i.e. multiply by {@link #setToScale(value_type, value_type, value_type) scale matrix}.
      * @param s scale for x-, y- and z-axis
      * @param tmp temporary mat4f used for multiplication
      * @return this matrix for chaining
      */
-    public final Mat4f scale(const float s, final Mat4f tmp) {
+    public final Matrix4 scale(const value_type s, final Matrix4 tmp) {
         return mul( tmp.setToScale(s, s, s) );
     }
 
@@ -1638,7 +1496,7 @@ class Mat4f {
     //
 
     /**
-     * Equals check using a given {@link FloatUtil#EPSILON} value and {@link FloatUtil#isEqual(float, float, float)}.
+     * Equals check using a given {@link FloatUtil#EPSILON} value and {@link FloatUtil#isEqual(value_type, value_type, value_type)}.
      * <p>
      * Implementation considers following corner cases:
      * <ul>
@@ -1650,7 +1508,7 @@ class Mat4f {
      * @param epsilon consider using {@link FloatUtil#EPSILON}
      * @return true if all components differ less than {@code epsilon}, otherwise false.
      */
-    public boolean isEqual(final Mat4f o, const float epsilon) {
+    public boolean isEqual(final Matrix4 o, const value_type epsilon) {
         if( this == o ) {
             return true;
         } else {
@@ -1674,7 +1532,7 @@ class Mat4f {
     }
 
     /**
-     * Equals check using {@link FloatUtil#EPSILON} value and {@link FloatUtil#isEqual(float, float, float)}.
+     * Equals check using {@link FloatUtil#EPSILON} value and {@link FloatUtil#isEqual(value_type, value_type, value_type)}.
      * <p>
      * Implementation considers following corner cases:
      * <ul>
@@ -1685,14 +1543,14 @@ class Mat4f {
      * @param o comparison value
      * @return true if all components differ less than {@link FloatUtil#EPSILON}, otherwise false.
      */
-    public boolean isEqual(final Mat4f o) {
+    public boolean isEqual(final Matrix4 o) {
         return isEqual(o, FloatUtil.EPSILON);
     }
 
     @Override
     public boolean equals(final Object o) {
-        if( o instanceof Mat4f ) {
-            return isEqual((Mat4f)o, FloatUtil.EPSILON);
+        if( o instanceof Matrix4 ) {
+            return isEqual((Matrix4)o, FloatUtil.EPSILON);
         } else {
             return false;
         }
@@ -1715,26 +1573,26 @@ class Mat4f {
      * @param winPos 3 component window coordinate, the result
      * @return true if successful, otherwise false (z is 1)
      */
-    public static boolean mapObjToWin(final Vec3f obj, final Mat4f mMv, final Mat4f mP,
+    public static boolean mapObjToWin(final Vec3f obj, final Matrix4 mMv, final Matrix4 mP,
                                       final Recti viewport, final Vec3f winPos)
     {
-        final Vec4f vec4Tmp1 = new Vec4f(obj, 1f);
+        final Vec4 vec4Tmp1 = new Vec4(obj, 1f);
 
         // vec4Tmp2 = Mv * o
         // rawWinPos = P  * vec4Tmp2
         // rawWinPos = P * ( Mv * o )
         // rawWinPos = P * Mv * o
-        final Vec4f vec4Tmp2 = mMv.mulVec4f(vec4Tmp1, new Vec4f());
-        final Vec4f rawWinPos = mP.mulVec4f(vec4Tmp2, vec4Tmp1);
+        final Vec4 vec4Tmp2 = mMv.mulVec4(vec4Tmp1, new Vec4());
+        final Vec4 rawWinPos = mP.mulVec4(vec4Tmp2, vec4Tmp1);
 
-        if (rawWinPos.w() == 0.0f) {
+        if (rawWinPos.w() == zero) {
             return false;
         }
 
-        const float s = ( 1.0f / rawWinPos.w() ) * 0.5f;
+        const value_type s = ( one / rawWinPos.w() ) * half;
 
         // Map x, y and z to range 0-1 (w is ignored)
-        rawWinPos.scale(s).add(0.5f, 0.5f, 0.5f, 0f);
+        rawWinPos.scale(s).add(half, half, half, 0f);
 
         // Map x,y to viewport
         winPos.set( rawWinPos.x() * viewport.width() +  viewport.x(),
@@ -1756,22 +1614,22 @@ class Mat4f {
      * @param winPos 3 component window coordinate, the result
      * @return true if successful, otherwise false (z is 1)
      */
-    public static boolean mapObjToWin(final Vec3f obj, final Mat4f mPMv,
+    public static boolean mapObjToWin(final Vec3f obj, final Matrix4 mPMv,
                                       final Recti viewport, final Vec3f winPos)
     {
-        final Vec4f vec4Tmp2 = new Vec4f(obj, 1f);
+        final Vec4 vec4Tmp2 = new Vec4(obj, 1f);
 
         // rawWinPos = P * Mv * o
-        final Vec4f rawWinPos = mPMv.mulVec4f(vec4Tmp2, new Vec4f());
+        final Vec4 rawWinPos = mPMv.mulVec4(vec4Tmp2, new Vec4());
 
-        if (rawWinPos.w() == 0.0f) {
+        if (rawWinPos.w() == zero) {
             return false;
         }
 
-        const float s = ( 1.0f / rawWinPos.w() ) * 0.5f;
+        const value_type s = ( one / rawWinPos.w() ) * half;
 
         // Map x, y and z to range 0-1 (w is ignored)
-        rawWinPos.scale(s).add(0.5f, 0.5f, 0.5f, 0f);
+        rawWinPos.scale(s).add(half, half, half, 0f);
 
         // Map x,y to viewport
         winPos.set( rawWinPos.x() * viewport.width() +  viewport.x(),
@@ -1797,19 +1655,19 @@ class Mat4f {
      * @param mat4Tmp 16 component matrix for temp storage
      * @return true if successful, otherwise false (failed to invert matrix, or becomes infinity due to zero z)
      */
-    public static boolean mapWinToObj(const float winx, const float winy, const float winz,
-                                      final Mat4f mMv, final Mat4f mP,
+    public static boolean mapWinToObj(const value_type winx, const value_type winy, const value_type winz,
+                                      final Matrix4 mMv, final Matrix4 mP,
                                       final Recti viewport,
                                       final Vec3f objPos,
-                                      final Mat4f mat4Tmp)
+                                      final Matrix4 mat4Tmp)
     {
         // invPMv = Inv(P x Mv)
-        final Mat4f invPMv = mat4Tmp.mul(mP, mMv);
+        final Matrix4 invPMv = mat4Tmp.mul(mP, mMv);
         if( !invPMv.invert() ) {
             return false;
         }
 
-        final Vec4f winPos = new Vec4f(winx, winy, winz, 1f);
+        final Vec4 winPos = new Vec4(winx, winy, winz, 1f);
 
         // Map x and y from window coordinates
         winPos.add(-viewport.x(), -viewport.y(), 0f, 0f).mul(1f/viewport.width(), 1f/viewport.height(), 1f, 1f);
@@ -1818,9 +1676,9 @@ class Mat4f {
         winPos.mul(2f, 2f, 2f, 1f).add(-1f, -1f, -1f, 0f);
 
         // rawObjPos = Inv(P x Mv) *  winPos
-        final Vec4f rawObjPos = invPMv.mulVec4f(winPos, new Vec4f());
+        final Vec4 rawObjPos = invPMv.mulVec4(winPos, new Vec4());
 
-        if ( rawObjPos.w() == 0.0f ) {
+        if ( rawObjPos.w() == zero ) {
             return false;
         }
         objPos.set( rawObjPos.scale( 1f / rawObjPos.w() ) );
@@ -1842,15 +1700,15 @@ class Mat4f {
      * @param objPos 3 component object coordinate, the result
      * @return true if successful, otherwise false (null invert matrix, or becomes infinity due to zero z)
      */
-    public static boolean mapWinToObj(const float winx, const float winy, const float winz,
-                                      final Mat4f invPMv,
+    public static boolean mapWinToObj(const value_type winx, const value_type winy, const value_type winz,
+                                      final Matrix4 invPMv,
                                       final Recti viewport,
                                       final Vec3f objPos)
     {
         if( null == invPMv ) {
             return false;
         }
-        final Vec4f winPos = new Vec4f(winx, winy, winz, 1f);
+        final Vec4 winPos = new Vec4(winx, winy, winz, 1f);
 
         // Map x and y from window coordinates
         winPos.add(-viewport.x(), -viewport.y(), 0f, 0f).mul(1f/viewport.width(), 1f/viewport.height(), 1f, 1f);
@@ -1859,9 +1717,9 @@ class Mat4f {
         winPos.mul(2f, 2f, 2f, 1f).add(-1f, -1f, -1f, 0f);
 
         // rawObjPos = Inv(P x Mv) *  winPos
-        final Vec4f rawObjPos = invPMv.mulVec4f(winPos, new Vec4f());
+        final Vec4 rawObjPos = invPMv.mulVec4(winPos, new Vec4());
 
-        if ( rawObjPos.w() == 0.0f ) {
+        if ( rawObjPos.w() == zero ) {
             return false;
         }
         objPos.set( rawObjPos.scale( 1f / rawObjPos.w() ) );
@@ -1885,15 +1743,15 @@ class Mat4f {
      * @param objPos1 3 component object coordinate, the result
      * @return true if successful, otherwise false (null invert matrix, or becomes infinity due to zero z)
      */
-    public static boolean mapWinToObj(const float winx, const float winy, const float winz1, const float winz2,
-                                      final Mat4f invPMv,
+    public static boolean mapWinToObj(const value_type winx, const value_type winy, const value_type winz1, const value_type winz2,
+                                      final Matrix4 invPMv,
                                       final Recti viewport,
                                       final Vec3f objPos1, final Vec3f objPos2)
     {
         if( null == invPMv ) {
             return false;
         }
-        final Vec4f winPos = new Vec4f(winx, winy, winz1, 1f);
+        final Vec4 winPos = new Vec4(winx, winy, winz1, 1f);
 
         // Map x and y from window coordinates
         winPos.add(-viewport.x(), -viewport.y(), 0f, 0f).mul(1f/viewport.width(), 1f/viewport.height(), 1f, 1f);
@@ -1902,9 +1760,9 @@ class Mat4f {
         winPos.mul(2f, 2f, 2f, 1f).add(-1f, -1f, -1f, 0f);
 
         // rawObjPos = Inv(P x Mv) *  winPos1
-        final Vec4f rawObjPos = invPMv.mulVec4f(winPos, new Vec4f());
+        final Vec4 rawObjPos = invPMv.mulVec4(winPos, new Vec4());
 
-        if ( rawObjPos.w() == 0.0f ) {
+        if ( rawObjPos.w() == zero ) {
             return false;
         }
         objPos1.set( rawObjPos.scale( 1f / rawObjPos.w() ) );
@@ -1916,9 +1774,9 @@ class Mat4f {
         winPos.setZ( winz2 * 2f - 1f );
 
         // rawObjPos = Inv(P x Mv) *  winPos2
-        invPMv.mulVec4f(winPos, rawObjPos);
+        invPMv.mulVec4(winPos, rawObjPos);
 
-        if ( rawObjPos.w() == 0.0f ) {
+        if ( rawObjPos.w() == zero ) {
             return false;
         }
         objPos2.set( rawObjPos.scale( 1f / rawObjPos.w() ) );
@@ -1945,20 +1803,20 @@ class Mat4f {
      * @param mat4Tmp 16 component matrix for temp storage
      * @return true if successful, otherwise false (failed to invert matrix, or becomes infinity due to zero z)
      */
-    public static boolean mapWinToObj4(const float winx, const float winy, const float winz, const float clipw,
-                                       final Mat4f mMv, final Mat4f mP,
+    public static boolean mapWinToObj4(const value_type winx, const value_type winy, const value_type winz, const value_type clipw,
+                                       final Matrix4 mMv, final Matrix4 mP,
                                        final Recti viewport,
-                                       const float near, const float far,
-                                       final Vec4f objPos,
-                                       final Mat4f mat4Tmp)
+                                       const value_type near, const value_type far,
+                                       final Vec4 objPos,
+                                       final Matrix4 mat4Tmp)
     {
         // invPMv = Inv(P x Mv)
-        final Mat4f invPMv = mat4Tmp.mul(mP, mMv);
+        final Matrix4 invPMv = mat4Tmp.mul(mP, mMv);
         if( !invPMv.invert() ) {
             return false;
         }
 
-        final Vec4f winPos = new Vec4f(winx, winy, winz, clipw);
+        final Vec4 winPos = new Vec4(winx, winy, winz, clipw);
 
         // Map x and y from window coordinates
         winPos.add(-viewport.x(), -viewport.y(), -near, 0f).mul(1f/viewport.width(), 1f/viewport.height(), 1f/(far-near), 1f);
@@ -1967,9 +1825,9 @@ class Mat4f {
         winPos.mul(2f, 2f, 2f, 1f).add(-1f, -1f, -1f, 0f);
 
         // objPos = Inv(P x Mv) *  winPos
-        invPMv.mulVec4f(winPos, objPos);
+        invPMv.mulVec4(winPos, objPos);
 
-        if ( objPos.w() == 0.0f ) {
+        if ( objPos.w() == zero ) {
             return false;
         }
         return true;
@@ -1992,16 +1850,16 @@ class Mat4f {
      * @param obj_pos 4 component object coordinate, the result
      * @return true if successful, otherwise false (null invert matrix, or becomes infinity due to zero z)
      */
-    public static boolean mapWinToObj4(const float winx, const float winy, const float winz, const float clipw,
-                                       final Mat4f invPMv,
+    public static boolean mapWinToObj4(const value_type winx, const value_type winy, const value_type winz, const value_type clipw,
+                                       final Matrix4 invPMv,
                                        final Recti viewport,
-                                       const float near, const float far,
-                                       final Vec4f objPos)
+                                       const value_type near, const value_type far,
+                                       final Vec4 objPos)
     {
         if( null == invPMv ) {
             return false;
         }
-        final Vec4f winPos = new Vec4f(winx, winy, winz, clipw);
+        final Vec4 winPos = new Vec4(winx, winy, winz, clipw);
 
         // Map x and y from window coordinates
         winPos.add(-viewport.x(), -viewport.y(), -near, 0f).mul(1f/viewport.width(), 1f/viewport.height(), 1f/(far-near), 1f);
@@ -2010,9 +1868,9 @@ class Mat4f {
         winPos.mul(2f, 2f, 2f, 1f).add(-1f, -1f, -1f, 0f);
 
         // objPos = Inv(P x Mv) *  winPos
-        invPMv.mulVec4f(winPos, objPos);
+        invPMv.mulVec4(winPos, objPos);
 
-        if ( objPos.w() == 0.0f ) {
+        if ( objPos.w() == zero ) {
             return false;
         }
         return true;
@@ -2021,13 +1879,13 @@ class Mat4f {
     /**
      * Map two window coordinates w/ shared X/Y and distinctive Z
      * to a {@link Ray}. The resulting {@link Ray} maybe used for <i>picking</i>
-     * using a {@link AABBox#getRayIntersection(vec3f, Ray, float, boolean)}.
+     * using a {@link AABBox#getRayIntersection(vec3f, Ray, value_type, boolean)}.
      * <p>
      * Notes for picking <i>winz0</i> and <i>winz1</i>:
      * <ul>
-     *   <li>see {@link FloatUtil#getZBufferEpsilon(int, float, float)}</li>
-     *   <li>see {@link FloatUtil#getZBufferValue(int, float, float, float)}</li>
-     *   <li>see {@link FloatUtil#getOrthoWinZ(float, float, float)}</li>
+     *   <li>see {@link FloatUtil#getZBufferEpsilon(int, value_type, value_type)}</li>
+     *   <li>see {@link FloatUtil#getZBufferValue(int, value_type, value_type, value_type)}</li>
+     *   <li>see {@link FloatUtil#getOrthoWinZ(value_type, value_type, value_type)}</li>
      * </ul>
      * </p>
      * @param winx
@@ -2042,13 +1900,13 @@ class Mat4f {
      * @param mat4Tmp2 16 component matrix for temp storage
      * @return true if successful, otherwise false (failed to invert matrix, or becomes z is infinity)
      */
-    public static boolean mapWinToRay(const float winx, const float winy, const float winz0, const float winz1,
-                                      final Mat4f mMv, final Mat4f mP,
+    public static boolean mapWinToRay(const value_type winx, const value_type winy, const value_type winz0, const value_type winz1,
+                                      final Matrix4 mMv, final Matrix4 mP,
                                       final Recti viewport,
                                       final Ray ray,
-                                      final Mat4f mat4Tmp1, final Mat4f mat4Tmp2) {
+                                      final Matrix4 mat4Tmp1, final Matrix4 mat4Tmp2) {
         // invPMv = Inv(P x Mv)
-        final Mat4f invPMv = mat4Tmp1.mul(mP, mMv);
+        final Matrix4 invPMv = mat4Tmp1.mul(mP, mMv);
         if( !invPMv.invert() ) {
             return false;
         }
@@ -2064,13 +1922,13 @@ class Mat4f {
     /**
      * Map two window coordinates w/ shared X/Y and distinctive Z
      * to a {@link Ray}. The resulting {@link Ray} maybe used for <i>picking</i>
-     * using a {@link AABBox#getRayIntersection(vec3f, Ray, float, boolean)}.
+     * using a {@link AABBox#getRayIntersection(vec3f, Ray, value_type, boolean)}.
      * <p>
      * Notes for picking <i>winz0</i> and <i>winz1</i>:
      * <ul>
-     *   <li>see {@link FloatUtil#getZBufferEpsilon(int, float, float)}</li>
-     *   <li>see {@link FloatUtil#getZBufferValue(int, float, float, float)}</li>
-     *   <li>see {@link FloatUtil#getOrthoWinZ(float, float, float)}</li>
+     *   <li>see {@link FloatUtil#getZBufferEpsilon(int, value_type, value_type)}</li>
+     *   <li>see {@link FloatUtil#getZBufferValue(int, value_type, value_type, value_type)}</li>
+     *   <li>see {@link FloatUtil#getOrthoWinZ(value_type, value_type, value_type)}</li>
      * </ul>
      * </p>
      * @param winx
@@ -2082,8 +1940,8 @@ class Mat4f {
      * @param ray storage for the resulting {@link Ray}
      * @return true if successful, otherwise false (null invert matrix, or becomes z is infinity)
      */
-    public static boolean mapWinToRay(const float winx, const float winy, const float winz0, const float winz1,
-                                      final Mat4f invPMv,
+    public static boolean mapWinToRay(const value_type winx, const value_type winy, const value_type winz0, const value_type winz1,
+                                      final Matrix4 invPMv,
                                       final Recti viewport,
                                       final Ray ray) {
         if( mapWinToObj(winx, winy, winz0, winz1, invPMv, viewport, ray.orig, ray.dir) ) {
@@ -2101,11 +1959,11 @@ class Mat4f {
     /**
      * @param sb optional passed StringBuilder instance to be used
      * @param rowPrefix optional prefix for each row
-     * @param f the format string of one floating point, i.e. "%10.5f", see {@link java.util.Formatter}
+     * @param f the format string of one value_typeing point, i.e. "%10.5f", see {@link java.util.Formatter}
      * @return matrix string representation
      */
     public StringBuilder toString(final StringBuilder sb, final String rowPrefix, final String f) {
-        const float[] tmp = new float[16];
+        const value_type[] tmp = new value_type[16];
         this.get(tmp);
         return FloatUtil.matrixToString(sb, rowPrefix, f,tmp, 0, 4, 4, false /* rowMajorOrder */);
     }
@@ -2119,18 +1977,47 @@ class Mat4f {
     /**
      * Returns a formatted string representation of this matrix
      * @param rowPrefix prefix for each row
-     * @param f format string for each float element, e.g. "%10.5f"
+     * @param f format string for each value_type element, e.g. "%10.5f"
      * @return matrix
      */
-    std::string toString(const std::string& rowPrefix, const std::string& f) const noexcept;
+    std::string toString(const std::string& rowPrefix, const std::string& f) const noexcept {
+        std::string sb;
+        value_type tmp[16];
+        get(tmp);
+        return jau::mat_to_string(sb, rowPrefix, f, tmp, 4, 4, false /* rowMajorOrder */); // creates a copy-out!
+    }
 
     std::string toString() const noexcept { return toString("", "%10.5f"); }
 };
 
-std::ostream& operator<<(std::ostream& out, const Mat4f& v) noexcept {
+template<typename T,
+         std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+constexpr Matrix4<T> operator*(const Matrix4<T>& lhs, const Matrix4<T>& rhs ) noexcept {
+    Matrix4<T> r(lhs); r *= rhs; return r;
+}
+
+template<typename T,
+         std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+constexpr Matrix4<T> operator*(const Matrix4<T>& lhs, const T s ) noexcept {
+    Matrix4<T> r(lhs); r *= s; return r;
+}
+
+template<typename T,
+         std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+constexpr Matrix4<T> operator*(const T s, const Matrix4<T>& rhs) noexcept {
+    Matrix4<T> r(rhs); r *= s; return r;
+}
+
+template<typename T,
+         std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+std::ostream& operator<<(std::ostream& out, const Matrix4<T>& v) noexcept {
     return out << v.toString();
 }
 
+typedef Matrix4<float> Mat4f;
+
+static_assert(alignof(float) == alignof(Mat4f));
+
 /**@}*/
 
 } // namespace jau::math
diff --git a/include/jau/math/quaternion.hpp b/include/jau/math/quaternion.hpp
index 6429ac5..a423c88 100644
--- a/include/jau/math/quaternion.hpp
+++ b/include/jau/math/quaternion.hpp
@@ -57,11 +57,24 @@ namespace jau::math {
  * See <a href="http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/index.htm">euclideanspace.com-Quaternion</a>
  * </p>
  */
-class Quaternion {
-  private:
-    float m_x, m_y, m_z, m_w;
-
+template<typename Value_type,
+         std::enable_if_t<std::is_floating_point_v<Value_type>, bool> = true>
+class alignas(Value_type) Quaternion {
   public:
+    typedef Value_type               value_type;
+    typedef value_type*              pointer;
+    typedef const value_type*        const_pointer;
+    typedef value_type&              reference;
+    typedef const value_type&        const_reference;
+    typedef value_type*              iterator;
+    typedef const value_type*        const_iterator;
+
+    typedef Vector3F<value_type, std::is_floating_point_v<Value_type>> Vec3;
+
+    constexpr static const value_type zero = value_type(0);
+    constexpr static const value_type one  = value_type(1);
+    constexpr static const value_type two  = value_type(2);
+    constexpr static const value_type half = one/two;
 
     /**
      * Quaternion Epsilon, used with equals method to determine if two Quaternions are close enough to be considered equal.
@@ -69,12 +82,17 @@ class Quaternion {
      * Using {@value}, which is ~10 times {@link FloatUtil#EPSILON}.
      * </p>
      */
-    constexpr static const float ALLOWED_DEVIANCE = 1.0E-6f; // FloatUtil.EPSILON == 1.1920929E-7f; double ALLOWED_DEVIANCE: 1.0E-8f
+    constexpr static const value_type allowed_deviation = 1.0E-6f; // FIXME: float EPSILON == 1.1920929E-7f; double ALLOWED_DEVIANCE: 1.0E-8f
+
+  private:
+    value_type m_x, m_y, m_z, m_w;
+
+  public:
 
     constexpr Quaternion() noexcept
     : m_x(0), m_y(0), m_z(0), m_w(1) {}
 
-    constexpr Quaternion(const float x, const float y, const float z, const float w) noexcept
+    constexpr Quaternion(const value_type x, const value_type y, const value_type z, const value_type w) noexcept
     : m_x(x), m_y(y), m_z(z), m_w(w) {}
 
     constexpr Quaternion(const Quaternion& o) noexcept = default;
@@ -87,7 +105,7 @@ class Quaternion {
      * which is not applied here.
      * @return the squared magnitude of this quaternion.
      */
-    float magnitudeSquared() const noexcept {
+    constexpr value_type magnitudeSquared() const noexcept {
         return m_w*m_w + m_x*m_x + m_y*m_y + m_z*m_z;
     }
 
@@ -100,13 +118,13 @@ class Quaternion {
      * <p>
      * Implementation Details:
      * <ul>
-     *   <li> returns 0f if {@link #magnitudeSquared()} is {@link jau::is_zero(float, float) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
-     *   <li> returns 1f if {@link #magnitudeSquared()} is {@link FloatUtil#isEqual(float, float, float) equals 1f} using {@link FloatUtil#EPSILON epsilon}</li>
+     *   <li> returns 0f if {@link #magnitudeSquared()} is {@link jau::is_zero(value_type, value_type) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
+     *   <li> returns 1f if {@link #magnitudeSquared()} is {@link FloatUtil#isEqual(value_type, value_type, value_type) equals 1f} using {@link FloatUtil#EPSILON epsilon}</li>
      * </ul>
      * </p>
      */
-    float magnitude() const noexcept {
-        const float magnitudeSQ = magnitudeSquared();
+    constexpr_cxx26 value_type magnitude() const noexcept {
+        const value_type magnitudeSQ = magnitudeSquared();
         if ( jau::is_zero(magnitudeSQ) ) {
             return 0.0f;
         }
@@ -116,33 +134,33 @@ class Quaternion {
         return std::sqrt(magnitudeSQ);
     }
 
-    float w() const noexcept { return m_w; }
+    constexpr value_type w() const noexcept { return m_w; }
 
-    void set_w(float w) noexcept { m_w = w; }
+    constexpr void set_w(value_type w) noexcept { m_w = w; }
 
-    float x() const noexcept { return m_x; }
+    constexpr value_type x() const noexcept { return m_x; }
 
-    void set_x(float x) noexcept { m_x = x; }
+    constexpr void set_x(value_type x) noexcept { m_x = x; }
 
-    float y() const noexcept { return m_y; }
+    constexpr value_type y() const noexcept { return m_y; }
 
-    void set_y(float y) noexcept { m_y = y; }
+    constexpr void set_y(value_type y) noexcept { m_y = y; }
 
-    float z() const noexcept { return m_z; }
+    constexpr value_type z() const noexcept { return m_z; }
 
-    void set_z(float z) noexcept { m_z = z; }
+    constexpr void set_z(value_type z) noexcept { m_z = z; }
 
     /**
      * Returns the dot product of this quaternion with the given x,y,z and m_w components.
      */
-    float dot(float x, float y, float z, float w) const noexcept {
+    constexpr value_type dot(value_type x, value_type y, value_type z, value_type w) const noexcept {
         return m_x * x + m_y * y + m_z * z + m_w * w;
     }
 
     /**
      * Returns the dot product of this quaternion with the given quaternion
      */
-    float dot(const Quaternion& quat) const noexcept {
+    constexpr value_type dot(const Quaternion& quat) const noexcept {
         return dot(quat.x(), quat.y(), quat.z(), quat.w());
     }
 
@@ -155,7 +173,7 @@ class Quaternion {
      * {@link jau::is_zero3f() against zero}.
      * </p>
      */
-    bool is_identity() const noexcept {
+    constexpr bool is_identity() const noexcept {
         return jau::equals(1.0f, m_w) && jau::is_zero3f(m_x, m_y, m_z);
         // return m_w == 1f && m_x == 0f && m_y == 0f && m_z == 0f;
     }
@@ -164,7 +182,7 @@ class Quaternion {
      * Set this quaternion to identity (x=0,y=0,z=0,w=1)
      * @return this quaternion for chaining.
      */
-    Quaternion& set_identity() noexcept {
+    constexpr Quaternion& set_identity() noexcept {
         m_x = m_y = m_z = 0.0f; m_w = 1.0f;
         return *this;
     }
@@ -179,12 +197,12 @@ class Quaternion {
      * </p>
      * @return this quaternion for chaining.
      */
-    Quaternion& normalize() noexcept {
-        const float norm = magnitude();
+    constexpr Quaternion& normalize() noexcept {
+        const value_type norm = magnitude();
         if ( jau::is_zero(norm) ) {
             set_identity();
         } else {
-            const float invNorm = 1.0f/norm;
+            const value_type invNorm = 1.0f/norm;
             m_w *= invNorm;
             m_x *= invNorm;
             m_y *= invNorm;
@@ -198,7 +216,7 @@ class Quaternion {
      * @return this quaternion for chaining.
      * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q49">Matrix-FAQ Q49</a>
      */
-    Quaternion& conjugate() noexcept {
+    constexpr Quaternion& conjugate() noexcept {
         m_x = -m_x;
         m_y = -m_y;
         m_z = -m_z;
@@ -210,18 +228,18 @@ class Quaternion {
      * <p>
      * Implementation Details:
      * <ul>
-     *   <li> {@link #conjugate() conjugates} if {@link #magnitudeSquared()} is is {@link FloatUtil#isEqual(float, float, float) equals 1f} using {@link FloatUtil#EPSILON epsilon}</li>
+     *   <li> {@link #conjugate() conjugates} if {@link #magnitudeSquared()} is is {@link FloatUtil#isEqual(value_type, value_type, value_type) equals 1f} using {@link FloatUtil#EPSILON epsilon}</li>
      * </ul>
      * </p>
      * @return this quaternion for chaining.
      * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q50">Matrix-FAQ Q50</a>
      */
-    Quaternion& invert() noexcept {
-        const float magnitudeSQ = magnitudeSquared();
+    constexpr Quaternion& invert() noexcept {
+        const value_type magnitudeSQ = magnitudeSquared();
         if ( jau::equals(1.0f, magnitudeSQ) ) {
             conjugate();
         } else {
-            const float invmsq = 1.0f/magnitudeSQ;
+            const value_type invmsq = 1.0f/magnitudeSQ;
             m_w *= invmsq;
             m_x = -m_x * invmsq;
             m_y = -m_y * invmsq;
@@ -234,7 +252,7 @@ class Quaternion {
      * Set all values of this quaternion using the given components.
      * @return this quaternion for chaining.
      */
-    constexpr Quaternion& set(const float x, const float y, const float z, const float w) noexcept {
+    constexpr Quaternion& set(const value_type x, const value_type y, const value_type z, const value_type w) noexcept {
         m_x = x;
         m_y = y;
         m_z = z;
@@ -290,11 +308,11 @@ class Quaternion {
     /**
      * Scale this quaternion by a scalar: this = this * rhs, returns this
      *
-     * @param n a float constant
+     * @param n a value_type constant
      * @return this quaternion for chaining.
      * @see <a href="http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/code/index.htm#scale">euclideanspace.com-QuaternionScale</a>
      */
-    constexpr Quaternion& operator*=(const float rhs) noexcept {
+    constexpr Quaternion& operator*=(const value_type rhs) noexcept {
         m_x *= rhs;
         m_y *= rhs;
         m_z *= rhs;
@@ -317,17 +335,17 @@ class Quaternion {
      * @param axisZ m_z-coord of rotation axis
      * @return this quaternion for chaining.
      */
-    Quaternion& rotateByAngleNormalAxis(float angle, float axisX, float axisY, float axisZ) noexcept {
+    constexpr_cxx26 Quaternion& rotateByAngleNormalAxis(const value_type angle, const value_type axisX, const value_type axisY, const value_type axisZ) noexcept {
         if( jau::is_zero3f(axisX, axisY, axisZ) ) {
             // no change
             return *this;
         }
-        const float halfAngle = 0.5f * angle;
-        const float sin = std::sin(halfAngle);
-        const float qw = std::cos(halfAngle);
-        const float qx = sin * axisX;
-        const float qy = sin * axisY;
-        const float qz = sin * axisZ;
+        const value_type halfAngle = 0.5f * angle;
+        const value_type sin = std::sin(halfAngle);
+        const value_type qw = std::cos(halfAngle);
+        const value_type qx = sin * axisX;
+        const value_type qy = sin * axisY;
+        const value_type qz = sin * axisZ;
         return set( m_x * qw + m_y * qz - m_z * qy + m_w * qx,
                    -m_x * qz + m_y * qw + m_z * qx + m_w * qy,
                     m_x * qy - m_y * qx + m_z * qw + m_w * qz,
@@ -344,10 +362,10 @@ class Quaternion {
      * </p>
      *
      * @param angle in radians
-     * @param axis  Vec3f coord of rotation axis
+     * @param axis  Vec3 coord of rotation axis
      * @return this quaternion for chaining.
      */
-    Quaternion& rotateByAngleNormalAxis(float angle, const Vec3f& axis) noexcept {
+    constexpr_cxx26 Quaternion& rotateByAngleNormalAxis(const value_type angle, const Vec3& axis) noexcept {
         return rotateByAngleNormalAxis(angle, axis.x, axis.y, axis.z);
     }
 
@@ -357,10 +375,12 @@ class Quaternion {
      * @param angle in radians
      * @return this quaternion for chaining.
      */
-    Quaternion& rotateByAngleX(float angle) noexcept {
-        const float halfAngle = 0.5f * angle;
-        const float sin = std::sin(halfAngle);
-        const float cos = std::cos(halfAngle);
+    constexpr_cxx26 Quaternion& rotateByAngleX(const value_type angle) noexcept {
+        const value_type halfAngle = 0.5f * angle;
+        return rotateByAngleX(std::sin(halfAngle), std::cos(halfAngle));
+    }
+    /** Rotate this quaternion around X axis with the given angle's sin + cos values */
+    constexpr Quaternion& rotateByAngleX(const value_type sin, const value_type cos) noexcept {
         return set( m_x * cos + m_w * sin,
                     m_y * cos + m_z * sin,
                    -m_y * sin + m_z * cos,
@@ -373,10 +393,12 @@ class Quaternion {
      * @param angle in radians
      * @return this quaternion for chaining.
      */
-    Quaternion& rotateByAngleY(float angle) noexcept {
-        const float halfAngle = 0.5f * angle;
-        const float sin = std::sin(halfAngle);
-        const float cos = std::cos(halfAngle);
+    constexpr_cxx26 Quaternion& rotateByAngleY(value_type angle) noexcept {
+        const value_type halfAngle = 0.5f * angle;
+        return rotateByAngleY(std::sin(halfAngle), std::cos(halfAngle));
+    }
+    /** Rotate this quaternion around Y axis with the given angle's sin + cos values */
+    constexpr Quaternion& rotateByAngleY(const value_type sin, const value_type cos) noexcept {
         return set( m_x * cos - m_z * sin,
                     m_y * cos + m_w * sin,
                     m_x * sin + m_z * cos,
@@ -389,10 +411,12 @@ class Quaternion {
      * @param angle in radians
      * @return this quaternion for chaining.
      */
-    Quaternion& rotateByAngleZ(float angle) noexcept {
-        const float halfAngle = 0.5f * angle;
-        const float sin = std::sin(halfAngle);
-        const float cos = std::cos(halfAngle);
+    constexpr_cxx26 Quaternion& rotateByAngleZ(value_type angle) noexcept {
+        const value_type halfAngle = 0.5f * angle;
+        return rotateByAngleZ(std::sin(halfAngle), std::cos(halfAngle));
+    }
+    /** Rotate this quaternion around Y axis with the given angle's sin + cos values */
+    constexpr Quaternion& rotateByAngleZ(const value_type sin, const value_type cos) noexcept {
         return set( m_x * cos + m_y * sin,
                    -m_x * sin + m_y * cos,
                     m_z * cos + m_w * sin,
@@ -409,12 +433,12 @@ class Quaternion {
      *  <li>z - attitude</li>
      * </ul>
      * </p>
-     * For details see {@link #rotateByEuler(float, float, float)}.
+     * For details see {@link #rotateByEuler(value_type, value_type, value_type)}.
      * @param angradXYZ euler angle array in radians
      * @return this quaternion for chaining.
-     * @see #rotateByEuler(float, float, float)
+     * @see #rotateByEuler(value_type, value_type, value_type)
      */
-    Quaternion& rotateByEuler(const Vec3f& angradXYZ) noexcept {
+    Quaternion& rotateByEuler(const Vec3& angradXYZ) noexcept {
         return rotateByEuler(angradXYZ.x, angradXYZ.y, angradXYZ.z);
     }
 
@@ -423,15 +447,15 @@ class Quaternion {
      * <p>
      * The rotations are applied in the given order and using chained rotation per axis:
      * <ul>
-     *  <li>y - heading  - {@link #rotateByAngleY(float)}</li>
-     *  <li>z - attitude - {@link #rotateByAngleZ(float)}</li>
-     *  <li>x - bank     - {@link #rotateByAngleX(float)}</li>
+     *  <li>y - heading  - {@link #rotateByAngleY(value_type)}</li>
+     *  <li>z - attitude - {@link #rotateByAngleZ(value_type)}</li>
+     *  <li>x - bank     - {@link #rotateByAngleX(value_type)}</li>
      * </ul>
      * </p>
      * <p>
      * Implementation Details:
      * <ul>
-     *   <li> NOP if all angles are {@link jau::is_zero(float, float) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
+     *   <li> NOP if all angles are {@link jau::is_zero(value_type, value_type) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
      *   <li> result is {@link #normalize()}ed</li>
      * </ul>
      * </p>
@@ -439,12 +463,12 @@ class Quaternion {
      * @param headingY the Euler yaw angle in radians. (rotation about the Y axis)
      * @param attitudeZ the Euler roll angle in radians. (rotation about the Z axis)
      * @return this quaternion for chaining.
-     * @see #rotateByAngleY(float)
-     * @see #rotateByAngleZ(float)
-     * @see #rotateByAngleX(float)
-     * @see #setFromEuler(float, float, float)
+     * @see #rotateByAngleY(value_type)
+     * @see #rotateByAngleZ(value_type)
+     * @see #rotateByAngleX(value_type)
+     * @see #setFromEuler(value_type, value_type, value_type)
      */
-    Quaternion& rotateByEuler(float bankX, float headingY, float attitudeZ) noexcept {
+    constexpr_cxx26 Quaternion& rotateByEuler(const value_type bankX, const value_type headingY, const value_type attitudeZ) noexcept {
         if ( jau::is_zero3f(bankX, headingY, attitudeZ) ) {
             return *this;
         } else {
@@ -461,18 +485,18 @@ class Quaternion {
      * @return the given out store for chaining
      * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q63">Matrix-FAQ Q63</a>
      */
-    Vec3f rotateVector(const Vec3f& in) noexcept {
-        Vec3f out;
+    Vec3 rotateVector(const Vec3& in) noexcept {
+        Vec3 out;
         if( in.is_zero() ) {
             out.set(0, 0, 0);
         } else {
-            const float vecX = in.x;
-            const float vecY = in.y;
-            const float vecZ = in.z;
-            const float x_x = m_x*m_x;
-            const float y_y = m_y*m_y;
-            const float z_z = m_z*m_z;
-            const float w_w = m_w*m_w;
+            const value_type vecX = in.x;
+            const value_type vecY = in.y;
+            const value_type vecZ = in.z;
+            const value_type x_x = m_x*m_x;
+            const value_type y_y = m_y*m_y;
+            const value_type z_z = m_z*m_z;
+            const value_type w_w = m_w*m_w;
 
             out.x =     w_w * vecX
                          + x_x * vecX
@@ -505,26 +529,26 @@ class Quaternion {
      *
      * @param a start quaternion
      * @param b end  quaternion
-     * @param changeAmnt float between 0 and 1 representing interpolation.
+     * @param changeAmnt value_type between 0 and 1 representing interpolation.
      * @return this quaternion for chaining.
      * @see <a href="http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/slerp/">euclideanspace.com-QuaternionSlerp</a>
      */
-    Quaternion& set_slerp(const Quaternion& a, const Quaternion& b, const float changeAmnt) noexcept {
+    constexpr_cxx26 Quaternion& set_slerp(const Quaternion& a, const Quaternion& b, const value_type changeAmnt) noexcept {
         // std::cerr << "Slerp.0: A " << a << ", B " << b << ", t " << changeAmnt << std::endl;
         if (changeAmnt == 0.0f) {
             *this = a;
         } else if (changeAmnt == 1.0f) {
             *this = b;
         } else {
-            float bx = b.m_x;
-            float by = b.m_y;
-            float bz = b.m_z;
-            float bw = b.m_w;
+            value_type bx = b.m_x;
+            value_type by = b.m_y;
+            value_type bz = b.m_z;
+            value_type bw = b.m_w;
 
             // Calculate angle between them (quat dot product)
-            float cosHalfTheta = a.m_x * bx + a.m_y * by + a.m_z * bz + a.m_w * bw;
+            value_type cosHalfTheta = a.m_x * bx + a.m_y * by + a.m_z * bz + a.m_w * bw;
 
-            float scale0, scale1;
+            value_type scale0, scale1;
 
             if( cosHalfTheta >= 0.95f ) {
                 // quaternions are close, just use linear interpolation
@@ -536,7 +560,7 @@ class Quaternion {
                 scale0 = 0.5f;
                 scale1 = 0.5f;
             } else {
-                if( cosHalfTheta <= -std::numeric_limits<float>::epsilon() ) { // FIXME: .. or shall we use the upper bound 'cosHalfTheta < EPSILON' ?
+                if( cosHalfTheta <= -std::numeric_limits<value_type>::epsilon() ) { // FIXME: .. or shall we use the upper bound 'cosHalfTheta < EPSILON' ?
                     // Negate the second quaternion and the result of the dot product (Inversion)
                     bx *= -1.0f;
                     by *= -1.0f;
@@ -544,8 +568,8 @@ class Quaternion {
                     bw *= -1.0f;
                     cosHalfTheta *= -1.0f;
                 }
-                const float halfTheta = std::acos(cosHalfTheta);
-                const float sinHalfTheta = std::sqrt(1.0f - cosHalfTheta*cosHalfTheta);
+                const value_type halfTheta = std::acos(cosHalfTheta);
+                const value_type sinHalfTheta = std::sqrt(1.0f - cosHalfTheta*cosHalfTheta);
                 // if theta = 180 degrees then result is not fully defined
                 // we could rotate around any axis normal to qa or qb
                 if ( std::abs(sinHalfTheta) < 0.001f ){ // fabs is floating point absolute
@@ -588,8 +612,8 @@ class Quaternion {
      * @return this quaternion for chaining.
      * @see <a href="http://www.euclideanspace.com/maths/algebra/vectors/lookat/index.htm">euclideanspace.com-LookUp</a>
      */
-    Quaternion& setLookAt(const Vec3f& directionIn, const Vec3f& upIn,
-                          Vec3f& xAxisOut, Vec3f& yAxisOut, Vec3f& zAxisOut) noexcept {
+    Quaternion& setLookAt(const Vec3& directionIn, const Vec3& upIn,
+                          Vec3& xAxisOut, Vec3& yAxisOut, Vec3& zAxisOut) noexcept {
         // Z = norm(dir)
         (zAxisOut = directionIn).normalize();
 
@@ -603,15 +627,15 @@ class Quaternion {
         yAxisOut.cross(zAxisOut, xAxisOut).normalize();
 
         /**
-            const float m00 = xAxisOut[0];
-            const float m01 = yAxisOut[0];
-            const float m02 = zAxisOut[0];
-            const float m10 = xAxisOut[1];
-            const float m11 = yAxisOut[1];
-            const float m12 = zAxisOut[1];
-            const float m20 = xAxisOut[2];
-            const float m21 = yAxisOut[2];
-            const float m22 = zAxisOut[2];
+            const value_type m00 = xAxisOut[0];
+            const value_type m01 = yAxisOut[0];
+            const value_type m02 = zAxisOut[0];
+            const value_type m10 = xAxisOut[1];
+            const value_type m11 = yAxisOut[1];
+            const value_type m12 = zAxisOut[1];
+            const value_type m20 = xAxisOut[2];
+            const value_type m21 = yAxisOut[2];
+            const value_type m22 = zAxisOut[2];
          */
         return setFromAxes(xAxisOut, yAxisOut, zAxisOut).normalize();
     }
@@ -630,22 +654,22 @@ class Quaternion {
      * <p>
      * Implementation Details:
      * <ul>
-     *   <li> set_identity() if square vector-length is jau::is_zero2f(float, float) using epsilon</li>
+     *   <li> set_identity() if square vector-length is jau::is_zero2f(value_type, value_type) using epsilon</li>
      * </ul>
      * </p>
      * @param v1 not normalized
      * @param v2 not normalized
      * @return this quaternion for chaining.
      */
-    Quaternion& setFromVectors(const Vec3f& v1, const Vec3f& v2) noexcept {
-        const float factor = v1.length() * v2.length();
+    constexpr_cxx26 Quaternion& setFromVectors(const Vec3& v1, const Vec3& v2) noexcept {
+        const value_type factor = v1.length() * v2.length();
         if ( jau::is_zero(factor) ) {
             return set_identity();
         } else {
-            const float dot = v1.dot(v2) / factor; // normalize
-            const float theta = std::acos(std::max(-1.0f, std::min(dot, 1.0f))); // clipping [-1..1]
+            const value_type dot = v1.dot(v2) / factor; // normalize
+            const value_type theta = std::acos(std::max(-1.0f, std::min(dot, 1.0f))); // clipping [-1..1]
 
-            Vec3f tmpPivotVec = v1.cross(v2);
+            Vec3 tmpPivotVec = v1.cross(v2);
 
             if ( dot < 0.0f && jau::is_zero( tmpPivotVec.length() ) ) {
                 // Vectors parallel and opposite direction, therefore a rotation of 180 degrees about any vector
@@ -684,22 +708,22 @@ class Quaternion {
      * <p>
      * Implementation Details:
      * <ul>
-     *   <li> {@link #setIdentity()} if square vector-length is {@link jau::is_zero(float, float) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
+     *   <li> {@link #setIdentity()} if square vector-length is {@link jau::is_zero(value_type, value_type) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
      * </ul>
      * </p>
      * @param v1 normalized
      * @param v2 normalized
      * @return this quaternion for chaining.
      */
-    Quaternion& setFromNormalVectors(const Vec3f& v1, const Vec3f& v2) noexcept {
-        const float factor = v1.length() * v2.length();
+    constexpr_cxx26 Quaternion& setFromNormalVectors(const Vec3& v1, const Vec3& v2) noexcept {
+        const value_type factor = v1.length() * v2.length();
         if ( jau::is_zero(factor) ) {
             return set_identity();
         } else {
-            const float dot = v1.dot(v2) / factor; // normalize
-            const float theta = std::acos(std::max(-1.0f, std::min(dot, 1.0f))); // clipping [-1..1]
+            const value_type dot = v1.dot(v2) / factor; // normalize
+            const value_type theta = std::acos(std::max(-1.0f, std::min(dot, 1.0f))); // clipping [-1..1]
 
-            Vec3f tmpPivotVec = v1.cross(v2);
+            Vec3 tmpPivotVec = v1.cross(v2);
 
             if ( dot < 0.0f && jau::is_zero( tmpPivotVec.length() ) ) {
                 // Vectors parallel and opposite direction, therefore a rotation of 180 degrees about any vector
@@ -733,7 +757,7 @@ class Quaternion {
      * <p>
      * Implementation Details:
      * <ul>
-     *   <li> {@link #setIdentity()} if axis is {@link jau::is_zero(float, float) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
+     *   <li> {@link #setIdentity()} if axis is {@link jau::is_zero(value_type, value_type) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
      * </ul>
      * </p>
      * @param angle rotation angle (rads)
@@ -741,10 +765,10 @@ class Quaternion {
      * @return this quaternion for chaining.
      *
      * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q56">Matrix-FAQ Q56</a>
-     * @see #toAngleAxis(Vec3f)
+     * @see #toAngleAxis(Vec3)
      */
-    Quaternion& setFromAngleAxis(const float angle, const Vec3f& vector) noexcept {
-        return setFromAngleNormalAxis(angle, Vec3f(vector).normalize());
+    Quaternion& setFromAngleAxis(const value_type angle, const Vec3& vector) noexcept {
+        return setFromAngleNormalAxis(angle, Vec3(vector).normalize());
     }
 
     /***
@@ -752,7 +776,7 @@ class Quaternion {
      * <p>
      * Implementation Details:
      * <ul>
-     *   <li> {@link #setIdentity()} if axis is {@link jau::is_zero(float, float) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
+     *   <li> {@link #setIdentity()} if axis is {@link jau::is_zero(value_type, value_type) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
      * </ul>
      * </p>
      * @param angle rotation angle (rads)
@@ -760,14 +784,14 @@ class Quaternion {
      * @return this quaternion for chaining.
      *
      * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q56">Matrix-FAQ Q56</a>
-     * @see #toAngleAxis(Vec3f)
+     * @see #toAngleAxis(Vec3)
      */
-    Quaternion& setFromAngleNormalAxis(const float angle, const Vec3f& vector) noexcept {
+    constexpr_cxx26 Quaternion& setFromAngleNormalAxis(const value_type angle, const Vec3& vector) noexcept {
         if( vector.is_zero() ) {
             set_identity();
         } else {
-            const float halfangle = angle * 0.5f;
-            const float sin = std::sin(halfangle);
+            const value_type halfangle = angle * 0.5f;
+            const value_type sin = std::sin(halfangle);
             m_x = vector.x * sin;
             m_y = vector.y * sin;
             m_z = vector.z * sin;
@@ -781,17 +805,17 @@ class Quaternion {
      *
      * @param axis storage for computed axis
      * @return the rotation angle in radians
-     * @see #setFromAngleAxis(float, Vec3f, Vec3f)
+     * @see #setFromAngleAxis(value_type, Vec3, Vec3)
      */
-    float toAngleAxis(Vec3f& axis) const noexcept {
-        const float sqrLength = m_x*m_x + m_y*m_y + m_z*m_z;
-        float angle;
+    constexpr_cxx26 value_type toAngleAxis(Vec3& axis) const noexcept {
+        const value_type sqrLength = m_x*m_x + m_y*m_y + m_z*m_z;
+        value_type angle;
         if ( jau::is_zero(sqrLength) ) { // length is ~0
             angle = 0.0f;
             axis.set( 1.0f, 0.0f, 0.0f );
         } else {
             angle = std::acos(m_w) * 2.0f;
-            const float invLength = 1.0f / std::sqrt(sqrLength);
+            const value_type invLength = 1.0f / std::sqrt(sqrLength);
             axis.set( m_x * invLength,
                       m_y * invLength,
                       m_z * invLength );
@@ -809,12 +833,12 @@ class Quaternion {
      *  <li>z - attitude</li>
      * </ul>
      * </p>
-     * For details see {@link #setFromEuler(float, float, float)}.
+     * For details see {@link #setFromEuler(value_type, value_type, value_type)}.
      * @param angradXYZ euler angle vector in radians holding x-bank, m_y-heading and m_z-attitude
      * @return this quaternion for chaining.
-     * @see #setFromEuler(float, float, float)
+     * @see #setFromEuler(value_type, value_type, value_type)
      */
-    Quaternion& setFromEuler(const Vec3f& angradXYZ) noexcept {
+    constexpr_cxx26 Quaternion& setFromEuler(const Vec3& angradXYZ) noexcept {
         return setFromEuler(angradXYZ.x, angradXYZ.y, angradXYZ.z);
     }
 
@@ -831,7 +855,7 @@ class Quaternion {
      * <p>
      * Implementation Details:
      * <ul>
-     *   <li> {@link #setIdentity()} if all angles are {@link jau::is_zero(float, float) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
+     *   <li> {@link #setIdentity()} if all angles are {@link jau::is_zero(value_type, value_type) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
      *   <li> result is {@link #normalize()}ed</li>
      * </ul>
      * </p>
@@ -843,27 +867,27 @@ class Quaternion {
      * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q60">Matrix-FAQ Q60</a>
      * @see <a href="http://vered.rose.utoronto.ca/people/david_dir/GEMS/GEMS.html">Gems</a>
      * @see <a href="http://www.euclideanspace.com/maths/geometry/rotations/conversions/eulerToQuaternion/index.htm">euclideanspace.com-eulerToQuaternion</a>
-     * @see #toEuler(Vec3f)
+     * @see #toEuler(Vec3)
      */
-    Quaternion& setFromEuler(const float bankX, const float headingY, const float attitudeZ) noexcept {
+    constexpr_cxx26 Quaternion& setFromEuler(const value_type bankX, const value_type headingY, const value_type attitudeZ) noexcept {
         if ( jau::is_zero3f(bankX, headingY, attitudeZ) ) {
             return set_identity();
         } else {
-            float angle = headingY * 0.5f;
-            const float sinHeadingY = std::sin(angle);
-            const float cosHeadingY = std::cos(angle);
+            value_type angle = headingY * 0.5f;
+            const value_type sinHeadingY = std::sin(angle);
+            const value_type cosHeadingY = std::cos(angle);
             angle = attitudeZ * 0.5f;
-            const float sinAttitudeZ = std::sin(angle);
-            const float cosAttitudeZ = std::cos(angle);
+            const value_type sinAttitudeZ = std::sin(angle);
+            const value_type cosAttitudeZ = std::cos(angle);
             angle = bankX * 0.5f;
-            const float sinBankX = std::sin(angle);
-            const float cosBankX = std::cos(angle);
+            const value_type sinBankX = std::sin(angle);
+            const value_type cosBankX = std::cos(angle);
 
             // variables used to reduce multiplication calls.
-            const float cosHeadingXcosAttitude = cosHeadingY * cosAttitudeZ;
-            const float sinHeadingXsinAttitude = sinHeadingY * sinAttitudeZ;
-            const float cosHeadingXsinAttitude = cosHeadingY * sinAttitudeZ;
-            const float sinHeadingXcosAttitude = sinHeadingY * cosAttitudeZ;
+            const value_type cosHeadingXcosAttitude = cosHeadingY * cosAttitudeZ;
+            const value_type sinHeadingXsinAttitude = sinHeadingY * sinAttitudeZ;
+            const value_type cosHeadingXsinAttitude = cosHeadingY * sinAttitudeZ;
+            const value_type sinHeadingXcosAttitude = sinHeadingY * cosAttitudeZ;
 
             m_w = cosHeadingXcosAttitude * cosBankX - sinHeadingXsinAttitude * sinBankX;
             m_x = cosHeadingXcosAttitude * sinBankX + sinHeadingXsinAttitude * cosBankX;
@@ -884,28 +908,28 @@ class Quaternion {
      * </ul>
      * </p>
      *
-     * @return new Vec3f euler angle result vector filled with x-bank, y-heading and m_z-attitude
+     * @return new Vec3 euler angle result vector filled with x-bank, y-heading and m_z-attitude
      * @see <a href="http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToEuler/index.htm">euclideanspace.com-quaternionToEuler</a>
-     * @see #setFromEuler(float, float, float)
+     * @see #setFromEuler(value_type, value_type, value_type)
      */
-    Vec3f toEuler() noexcept {
-        const float sqw = m_w*m_w;
-        const float sqx = m_x*m_x;
-        const float sqy = m_y*m_y;
-        const float sqz = m_z*m_z;
-        const float unit = sqx + sqy + sqz + sqw; // if normalized is one, otherwise, is correction factor
-        const float test = m_x*m_y + m_z*m_w;
+    constexpr_cxx26 Vec3 toEuler() noexcept {
+        const value_type sqw = m_w*m_w;
+        const value_type sqx = m_x*m_x;
+        const value_type sqy = m_y*m_y;
+        const value_type sqz = m_z*m_z;
+        const value_type unit = sqx + sqy + sqz + sqw; // if normalized is one, otherwise, is correction factor
+        const value_type test = m_x*m_y + m_z*m_w;
 
         if (test > 0.499f * unit) { // singularity at north pole
-            return Vec3f( 0.0f,                        // m_x-bank
+            return Vec3( 0.0f,                        // m_x-bank
                           2.0f * std::atan2(m_x, m_w), // y-heading
                           M_PI_2 );                    // z-attitude
         } else if (test < -0.499f * unit) { // singularity at south pole
-            return Vec3f( 0.0f,                        // m_x-bank
+            return Vec3( 0.0f,                        // m_x-bank
                          -2.0 * std::atan2(m_x, m_w),  // m_y-heading
                          -M_PI_2 );                    // m_z-attitude
         } else {
-            return Vec3f( std::atan2(2.0f * m_x * m_w - 2.0f * m_y * m_z, -sqx + sqy - sqz + sqw), // m_x-bank
+            return Vec3( std::atan2(2.0f * m_x * m_w - 2.0f * m_y * m_z, -sqx + sqy - sqz + sqw), // m_x-bank
                           std::atan2(2.0f * m_y * m_w - 2.0f * m_x * m_z,  sqx - sqy - sqz + sqw), // m_y-heading
                           std::asin( 2.0f * test / unit) );                                        // z-attitude
         }
@@ -924,35 +948,35 @@ class Quaternion {
      * @return this quaternion for chaining.
      * @see #setFromMatrix(Matrix4f)
      */
-    Quaternion& setFromMat3(const float m00, const float m01, const float m02,
-                            const float m10, const float m11, const float m12,
-                            const float m20, const float m21, const float m22) noexcept {
+    constexpr Quaternion& setFromMat(const value_type m00, const value_type m01, const value_type m02,
+                                     const value_type m10, const value_type m11, const value_type m12,
+                                     const value_type m20, const value_type m21, const value_type m22) noexcept {
         // Note: Other implementations uses 'T' w/o '+1f' and compares 'T >= 0' while adding missing 1f in sqrt expr.
         //       However .. this causes setLookAt(..) to fail and actually violates the 'trace definition'.
 
         // The trace T is the sum of the diagonal elements; see
         // http://mathworld.wolfram.com/MatrixTrace.html
-        const float T = m00 + m11 + m22 + 1.0f;
+        const value_type T = m00 + m11 + m22 + 1.0f;
         if ( T > 0.0f ) {
-            const float S = 0.5f / std::sqrt(T);  // S = 1 / ( 2 t )
+            const value_type S = 0.5f / std::sqrt(T);  // S = 1 / ( 2 t )
             m_w = 0.25f / S;                      // m_w = 1 / ( 4 S ) = t / 2
             m_x = ( m21 - m12 ) * S;
             m_y = ( m02 - m20 ) * S;
             m_z = ( m10 - m01 ) * S;
         } else if ( m00 > m11 && m00 > m22) {
-            const float S = 0.5f / std::sqrt(1.0f + m00 - m11 - m22); // S=4*qx
+            const value_type S = 0.5f / std::sqrt(1.0f + m00 - m11 - m22); // S=4*qx
             m_w = ( m21 - m12 ) * S;
             m_x = 0.25f / S;
             m_y = ( m10 + m01 ) * S;
             m_z = ( m02 + m20 ) * S;
         } else if ( m11 > m22 ) {
-            const float S = 0.5f / std::sqrt(1.0f + m11 - m00 - m22); // S=4*qy
+            const value_type S = 0.5f / std::sqrt(1.0f + m11 - m00 - m22); // S=4*qy
             m_w = ( m02 - m20 ) * S;
             m_x = ( m20 + m01 ) * S;
             m_y = 0.25f / S;
             m_z = ( m21 + m12 ) * S;
         } else {
-            const float S = 0.5f / std::sqrt(1.0f + m22 - m00 - m11); // S=4*qz
+            const value_type S = 0.5f / std::sqrt(1.0f + m22 - m00 - m11); // S=4*qz
             m_w = ( m10 - m01 ) * S;
             m_x = ( m02 + m20 ) * S;
             m_y = ( m21 + m12 ) * S;
@@ -973,10 +997,27 @@ class Quaternion {
      *
      * @return this quaternion for chaining.
      * @see Matrix4f#getRotation(Quaternion)
-     * @see #setFromMatrix(float, float, float, float, float, float, float, float, float)
+     * @see #setFromMatrix(value_type, value_type, value_type, value_type, value_type, value_type, value_type, value_type, value_type)
      */
-    Quaternion& setFromMat3(const Mat4f& m) noexcept {
-        return setFromMat3(m.m00, m.m01, m.m02, m.m10, m.m11, m.m12, m.m20, m.m21, m.m22);
+    constexpr Quaternion& setFromMat(const Mat4f& m) noexcept {
+        return setFromMat(m.m00, m.m01, m.m02, m.m10, m.m11, m.m12, m.m20, m.m21, m.m22);
+    }
+
+    /**
+     * Initializes this quaternion to represent a rotation formed by the given three <i>orthogonal</i> axes.
+     * <p>
+     * No validation whether the axes are <i>orthogonal</i> is performed.
+     * </p>
+     *
+     * @param xAxis vector representing the <i>orthogonal</i> x-axis of the coordinate system.
+     * @param yAxis vector representing the <i>orthogonal</i> y-axis of the coordinate system.
+     * @param zAxis vector representing the <i>orthogonal</i> m_z-axis of the coordinate system.
+     * @return this quaternion for chaining.
+     */
+    constexpr Quaternion& setFromAxes(const Vec3& xAxis, const Vec3& yAxis, const Vec3& zAxis) noexcept {
+        return setFromMat(xAxis.x, yAxis.x, zAxis.x,
+                          xAxis.y, yAxis.y, zAxis.y,
+                          xAxis.z, yAxis.z, zAxis.z);
     }
 
     /**
@@ -984,20 +1025,19 @@ class Quaternion {
      * <p>
      * Implementation Details:
      * <ul>
-     *   <li> makes identity matrix if {@link #magnitudeSquared()} is {@link jau::is_zero(float, float) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
+     *   <li> makes identity matrix if {@link #magnitudeSquared()} is {@link jau::is_zero(value_type, value_type) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
      * </ul>
      * </p>
      *
      * @return resulting normalized column matrix 4x4
      * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q54">Matrix-FAQ Q54</a>
-     * @see #setFromMatrix(float, float, float, float, float, float, float, float, float)
+     * @see #setFromMatrix(value_type, value_type, value_type, value_type, value_type, value_type, value_type, value_type, value_type)
      * @see Matrix4f#setToRotation(Quaternion)
      */
     Mat4f toMatrix() const noexcept {
-        Mat4f m;
-        m.setToRotation(*this);
-        return m;
+        Mat4f m; toMatrix(m); return m;
     }
+
     /**
      * Transform this quaternion to a normalized 4x4 column matrix representing the rotation,
      * see toMatrix().
@@ -1005,25 +1045,72 @@ class Quaternion {
      * @param out store for the resulting normalized column matrix 4x4
      * @return the given matrix store
      */
-    Mat4f& toMatrix(Mat4f& out) const noexcept {
-        return out.setToRotation(*this);
+    Mat4f& toMatrix(Mat4f& m) const noexcept {
+        // pre-multiply scaled-reciprocal-magnitude to reduce multiplications
+        const value_type norm = magnitudeSquared();
+        if ( jau::is_zero(norm) ) {
+            // identity matrix -> srecip = 0f
+            m.loadIdentity();
+            return m;
+        }
+        value_type srecip;
+        if ( jau::equals(1.0f, norm) ) {
+            srecip = 2.0f;
+        } else {
+            srecip = 2.0f / norm;
+        }
+        const value_type x = m_x;
+        const value_type y = m_y;
+        const value_type z = m_z;
+        const value_type w = m_w;
+
+        const value_type xs = srecip * x;
+        const value_type ys = srecip * y;
+        const value_type zs = srecip * z;
+
+        const value_type xx = x  * xs;
+        const value_type xy = x  * ys;
+        const value_type xz = x  * zs;
+        const value_type xw = xs * w;
+        const value_type yy = y  * ys;
+        const value_type yz = y  * zs;
+        const value_type yw = ys * w;
+        const value_type zz = z  * zs;
+        const value_type zw = zs * w;
+
+        m.m00 = 1.0f - ( yy + zz );
+        m.m01 =        ( xy - zw );
+        m.m02 =        ( xz + yw );
+        m.m03 = 0.0f;
+
+        m.m10 =        ( xy + zw );
+        m.m11 = 1.0f - ( xx + zz );
+        m.m12 =        ( yz - xw );
+        m.m13 = 0.0f;
+
+        m.m20 =        ( xz - yw );
+        m.m21 =        ( yz + xw );
+        m.m22 = 1.0f - ( xx + yy );
+        m.m23 = 0.0f;
+
+        m.m30 = m.m31 = m.m32 = 0.0f;
+        m.m33 = 1.0f;
+        return m;
     }
 
     /**
-     * Initializes this quaternion to represent a rotation formed by the given three <i>orthogonal</i> axes.
-     * <p>
-     * No validation whether the axes are <i>orthogonal</i> is performed.
-     * </p>
+     * Extracts this quaternion's <i>orthogonal</i> rotation axes.
      *
      * @param xAxis vector representing the <i>orthogonal</i> x-axis of the coordinate system.
      * @param yAxis vector representing the <i>orthogonal</i> y-axis of the coordinate system.
      * @param zAxis vector representing the <i>orthogonal</i> m_z-axis of the coordinate system.
-     * @return this quaternion for chaining.
+     * @param tmp temporary Matrix4 used for toMatrix()
      */
-    Quaternion& setFromAxes(const Vec3f& xAxis, const Vec3f& yAxis, const Vec3f& zAxis) noexcept {
-        return setFromMat3(xAxis.x, yAxis.x, zAxis.x,
-                           xAxis.y, yAxis.y, zAxis.y,
-                           xAxis.z, yAxis.z, zAxis.z);
+    void toAxes(Vec3& xAxis, Vec3& yAxis, Vec3& zAxis, Matrix4<value_type>& tmp) const noexcept {
+        toMatrix(tmp);
+        tmp.getColumn(2, zAxis);
+        tmp.getColumn(1, yAxis);
+        tmp.getColumn(0, xAxis);
     }
 
     /**
@@ -1033,12 +1120,9 @@ class Quaternion {
      * @param yAxis vector representing the <i>orthogonal</i> y-axis of the coordinate system.
      * @param zAxis vector representing the <i>orthogonal</i> m_z-axis of the coordinate system.
      */
-    void toAxes(Vec3f& xAxis, Vec3f& yAxis, Vec3f& zAxis) const noexcept {
-        Mat4f tmpMat4;
-        tmpMat4.setToRotation(*this);
-        tmpMat4.getColumn(2, zAxis);
-        tmpMat4.getColumn(1, yAxis);
-        tmpMat4.getColumn(0, xAxis);
+    void toAxes(Vec3& xAxis, Vec3& yAxis, Vec3& zAxis) const noexcept {
+        Matrix4<value_type> tmp;
+        toAxes(xAxis, yAxis, zAxis, tmp);
     }
 
     //
@@ -1049,14 +1133,14 @@ class Quaternion {
      * @param o the object to compare for equality
      * @return true if this quaternion and the provided quaternion have roughly the same x, m_y, m_z and m_w values.
      */
-    constexpr bool equals(const Quaternion& o) const noexcept {
+    constexpr bool operator==(const Quaternion& o) const noexcept {
         if (this == &o) {
             return true;
         }
-        return std::abs(m_x - o.m_x) <= ALLOWED_DEVIANCE &&
-               std::abs(m_y - o.m_y) <= ALLOWED_DEVIANCE &&
-               std::abs(m_z - o.m_z) <= ALLOWED_DEVIANCE &&
-               std::abs(m_w - o.m_w) <= ALLOWED_DEVIANCE;
+        return std::abs(m_x - o.m_x) <= allowed_deviation &&
+               std::abs(m_y - o.m_y) <= allowed_deviation &&
+               std::abs(m_z - o.m_z) <= allowed_deviation &&
+               std::abs(m_w - o.m_w) <= allowed_deviation;
     }
 
     std::string toString() const noexcept {
@@ -1064,35 +1148,45 @@ class Quaternion {
     }
 };
 
-constexpr bool operator==(const Quaternion& lhs, const Quaternion& rhs ) noexcept {
-    return lhs.equals(rhs);
-}
-
-constexpr Quaternion operator+(const Quaternion& lhs, const Quaternion& rhs ) noexcept {
-    return Quaternion(lhs) += rhs;
+template<typename T,
+         std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+constexpr Quaternion<T> operator+(const Quaternion<T>& lhs, const Quaternion<T>& rhs ) noexcept {
+    Quaternion<T> r(lhs); r += rhs; return r;
 }
 
-constexpr Quaternion operator-(const Quaternion& lhs, const Quaternion& rhs ) noexcept {
-    return Quaternion(lhs) -= rhs;
+template<typename T,
+         std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+constexpr Quaternion<T> operator-(const Quaternion<T>& lhs, const Quaternion<T>& rhs ) noexcept {
+    Quaternion<T> r(lhs); r -= rhs; return r;
 }
 
-constexpr Quaternion operator*(const Quaternion& lhs, const Quaternion& rhs ) noexcept {
-    return Quaternion(lhs) *= rhs;
+template<typename T,
+         std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+constexpr Quaternion<T> operator*(const Quaternion<T>& lhs, const Quaternion<T>& rhs ) noexcept {
+    Quaternion<T> r(lhs); r *= rhs; return r;
 }
 
-constexpr Quaternion operator*(const Quaternion& lhs, const float s ) noexcept {
-    return Quaternion(lhs) *= s;
+template<typename T,
+         std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+constexpr Quaternion<T> operator*(const Quaternion<T>& lhs, const T s ) noexcept {
+    Quaternion<T> r(lhs); r *= s; return r;
 }
 
-constexpr Quaternion operator*(const float s, const Quaternion& rhs) noexcept {
-    return Quaternion(rhs) *= s;
+template<typename T,
+         std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+constexpr Quaternion<T> operator*(const T s, const Quaternion<T>& rhs) noexcept {
+    Quaternion<T> r(rhs); r *= s; return r;
 }
 
-
-std::ostream& operator<<(std::ostream& out, const Quaternion& v) noexcept {
+template<typename T,
+         std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+std::ostream& operator<<(std::ostream& out, const Quaternion<T>& v) noexcept {
     return out << v.toString();
 }
 
+typedef Quaternion<float> Quat4f;
+static_assert(alignof(float) == alignof(Quat4f));
+
 /**@}*/
 
 } // namespace jau::math
diff --git a/include/jau/math/recti.hpp b/include/jau/math/recti.hpp
index 43a9aba..846be22 100644
--- a/include/jau/math/recti.hpp
+++ b/include/jau/math/recti.hpp
@@ -38,31 +38,42 @@ namespace jau::math {
     /**
      * Rectangle with x, y, width and height integer components.
      */
-    class Recti {
+    template<typename Value_type,
+             std::enable_if_t<std::is_integral_v<Value_type>, bool> = true>
+    class RectI {
+      public:
+        typedef Value_type                  value_type;
+        typedef value_type*                 pointer;
+        typedef const value_type*           const_pointer;
+        typedef value_type&                 reference;
+        typedef const value_type&           const_reference;
+        typedef value_type*                 iterator;
+        typedef const value_type*           const_iterator;
+
       private:
-        int m_x;
-        int m_y;
-        int m_width;
-        int m_height;
+        value_type m_x;
+        value_type m_y;
+        value_type m_width;
+        value_type m_height;
 
       public:
-        constexpr Recti() noexcept
+        constexpr RectI() noexcept
         : m_x(0), m_y(0), m_width(0), m_height(0) {}
 
-        constexpr Recti(const int xywh[/*4*/]) noexcept{
+        constexpr RectI(const value_type xywh[/*4*/]) noexcept{
             set(xywh);
         }
 
-        constexpr Recti(const int x, const int y, const int width, const int height) noexcept {
+        constexpr RectI(const value_type x, const value_type y, const value_type width, const value_type height) noexcept {
             set(x, y, width, height);
         }
 
-        constexpr Recti(const Recti& o) noexcept = default;
-        constexpr Recti(Recti&& o) noexcept = default;
-        constexpr Recti& operator=(const Recti&) noexcept = default;
-        constexpr Recti& operator=(Recti&&) noexcept = default;
+        constexpr RectI(const RectI& o) noexcept = default;
+        constexpr RectI(RectI&& o) noexcept = default;
+        constexpr RectI& operator=(const RectI&) noexcept = default;
+        constexpr RectI& operator=(RectI&&) noexcept = default;
 
-        constexpr bool operator==(const Recti& rhs ) const noexcept {
+        constexpr bool operator==(const RectI& rhs ) const noexcept {
             if( this == &rhs ) {
                 return true;
             }
@@ -76,41 +87,26 @@ namespace jau::math {
         } */
 
         /** this = { x, y, width, height }, returns this. */
-        constexpr Recti& set(const int x, const int y, const int width, const int height) noexcept {
-            m_x = x;
-            m_y = y;
-            m_width = width;
-            m_height= height;
-            return *this;
-        }
+        constexpr RectI& set(const value_type x, const value_type y, const value_type width, const value_type height) noexcept
+        { m_x = x; m_y = y; m_width = width; m_height= height; return *this; }
 
         /** this = xywh, returns this. */
-        constexpr Recti& set(const int xywh[/*4*/]) noexcept {
-            m_x = xywh[0];
-            m_y = xywh[1];
-            m_width = xywh[2];
-            m_height= xywh[3];
-            return *this;
-        }
+        constexpr RectI& set(const_iterator xywh) noexcept
+        { m_x = xywh[0]; m_y = xywh[1]; m_width = xywh[2]; m_height= xywh[3]; return *this; }
 
         /** xywh = this, returns xywh. */
-        int* get(int xywh[/*4*/]) const noexcept {
-            xywh[0] = m_x;
-            xywh[1] = m_y;
-            xywh[2] = m_width;
-            xywh[3] = m_height;
-            return xywh;
-        }
+        iterator get(iterator xywh) const noexcept
+        { xywh[0] = m_x; xywh[1] = m_y; xywh[2] = m_width; xywh[3] = m_height; return xywh; }
 
-        int x() const noexcept { return m_x; }
-        int y() const noexcept { return m_y; }
-        int width() const noexcept { return m_width; }
-        int height() const noexcept { return m_height; }
+        value_type x() const noexcept { return m_x; }
+        value_type y() const noexcept { return m_y; }
+        value_type width() const noexcept { return m_width; }
+        value_type height() const noexcept { return m_height; }
 
-        void setX(const int x) noexcept { m_x = x; }
-        void setY(const int y) noexcept { m_y = y; }
-        void setWidth(const int width) noexcept { m_width = width; }
-        void setHeight(const int height) noexcept { m_height = height; }
+        void setX(const value_type x) noexcept { m_x = x; }
+        void setY(const value_type y) noexcept { m_y = y; }
+        void setWidth(const value_type width) noexcept { m_width = width; }
+        void setHeight(const value_type height) noexcept { m_height = height; }
 
         /** Return true if area is zero. */
         bool is_zero() const noexcept {
@@ -121,10 +117,15 @@ namespace jau::math {
         { return std::to_string(m_x)+"/"+std::to_string(m_y)+" "+std::to_string(m_width)+"x"+std::to_string(m_height); }
     };
 
-    std::ostream& operator<<(std::ostream& out, const Recti& v) noexcept {
+    template<typename T,
+             std::enable_if_t<std::numeric_limits<T>::is_integer, bool> = true>
+    std::ostream& operator<<(std::ostream& out, const RectI<T>& v) noexcept {
         return out << v.toString();
     }
 
+    typedef RectI<int> Recti;
+    static_assert(alignof(int) == alignof(Recti));
+
 /**@}*/
 
 } // namespace jau::math
diff --git a/include/jau/math/vec2f.hpp b/include/jau/math/vec2f.hpp
index 20943bd..e4e954a 100644
--- a/include/jau/math/vec2f.hpp
+++ b/include/jau/math/vec2f.hpp
@@ -27,6 +27,7 @@
 #include <cmath>
 #include <cstdarg>
 #include <cstdint>
+#include <cassert>
 #include <limits>
 #include <string>
 #include <iostream>
@@ -41,46 +42,64 @@ namespace jau::math {
      */
 
     /**
-     * 2D vector using two float components.
+     * 2D vector using two value_type components.
      */
-    class Vec2f {
-    public:
-        float x;
-        float y;
-
-        static constexpr Vec2f from_length_angle(const float magnitude, const float radians) noexcept {
-            return Vec2f(magnitude * std::cos(radians), magnitude * std::sin(radians));
+    template<typename Value_type,
+             std::enable_if_t<std::is_floating_point_v<Value_type>, bool> = true>
+    class alignas(Value_type) Vector2F {
+      public:
+        typedef Value_type                  value_type;
+        typedef value_type*                 pointer;
+        typedef const value_type*           const_pointer;
+        typedef value_type&                 reference;
+        typedef const value_type&           const_reference;
+        typedef value_type*                 iterator;
+        typedef const value_type*           const_iterator;
+
+        constexpr static const value_type zero = value_type(0);
+        constexpr static const value_type one  = value_type(1);
+
+        value_type x;
+        value_type y;
+
+        static constexpr_cxx26 Vector2F from_length_angle(const value_type magnitude, const value_type radians) noexcept {
+            return Vector2F(magnitude * std::cos(radians), magnitude * std::sin(radians));
         }
 
-        constexpr Vec2f() noexcept
-        : x(0), y(0) {}
+        constexpr Vector2F() noexcept
+        : x(zero), y(zero) {}
 
-        constexpr Vec2f(const float x_, const float y_) noexcept
+        constexpr Vector2F(const value_type v) noexcept
+        : x(v), y(v) {}
+
+        constexpr Vector2F(const value_type x_, const value_type y_) noexcept
         : x(x_), y(y_) {}
 
-        constexpr Vec2f(const Vec2f& o) noexcept = default;
-        constexpr Vec2f(Vec2f&& o) noexcept = default;
-        constexpr Vec2f& operator=(const Vec2f&) noexcept = default;
-        constexpr Vec2f& operator=(Vec2f&&) noexcept = default;
+        constexpr Vector2F(const Vector2F& o) noexcept = default;
+        constexpr Vector2F(Vector2F&& o) noexcept = default;
+        constexpr Vector2F& operator=(const Vector2F&) noexcept = default;
+        constexpr Vector2F& operator=(Vector2F&&) noexcept = default;
 
-        /** Returns read-only component w/o boundary check */
-        float operator[](size_t i) const noexcept {
-            return reinterpret_cast<const float*>(this)[i];
+        /** Returns read-only component */
+        constexpr value_type operator[](size_t i) const noexcept {
+            assert(i < 2);
+            return (&x)[i];
         }
 
-        /** Returns writeable reference to component w/o boundary check */
-        float& operator[](size_t i) noexcept {
-            return reinterpret_cast<float*>(this)[i];
+        /** Returns writeable reference to component */
+        constexpr reference operator[](size_t i) noexcept {
+            assert(i < 2);
+            return (&x)[i];
         }
 
         /** xy = this, returns xy. */
-        float* get(float xy[/*2*/]) const noexcept {
+        constexpr iterator get(iterator xy) const noexcept {
             xy[0] = x;
             xy[1] = y;
             return xy;
         }
 
-        constexpr bool operator==(const Vec2f& rhs ) const noexcept {
+        constexpr bool operator==(const Vector2F& rhs ) const noexcept {
             if( this == &rhs ) {
                 return true;
             }
@@ -91,18 +110,22 @@ namespace jau::math {
             return ...
         } */
 
-        constexpr Vec2f& set(const float vx, const float vy) noexcept
+        constexpr Vector2F& set(const value_type vx, const value_type vy) noexcept
         { x=vx; y=vy; return *this; }
 
-        constexpr Vec2f& add(const float dx, const float dy) noexcept
+        /** this = xy, returns this. */
+        constexpr Vector2F& set(const_iterator xy) noexcept
+        { x=xy[0]; y=xy[1]; return *this; }
+
+        constexpr Vector2F& add(const value_type dx, const value_type dy) noexcept
         { x+=dx; y+=dy; return *this; }
 
-        constexpr Vec2f& operator+=(const Vec2f& rhs ) noexcept {
+        constexpr Vector2F& operator+=(const Vector2F& rhs ) noexcept {
             x+=rhs.x; y+=rhs.y;
             return *this;
         }
 
-        constexpr Vec2f& operator-=(const Vec2f& rhs ) noexcept {
+        constexpr Vector2F& operator-=(const Vector2F& rhs ) noexcept {
             x-=rhs.x; y-=rhs.y;
             return *this;
         }
@@ -112,7 +135,7 @@ namespace jau::math {
          * @param s scale factor
          * @return this instance
          */
-        constexpr Vec2f& operator*=(const float s ) noexcept {
+        constexpr Vector2F& operator*=(const value_type s ) noexcept {
             x*=s; y*=s;
             return *this;
         }
@@ -122,33 +145,33 @@ namespace jau::math {
          * @param s scale factor
          * @return this instance
          */
-        constexpr Vec2f& operator/=(const float s ) noexcept {
+        constexpr Vector2F& operator/=(const value_type s ) noexcept {
             x/=s; y/=s;
             return *this;
         }
 
         /** Rotates this vector in place, returns *this */
-        Vec2f& rotate(const float radians, const Vec2f& ctr) noexcept {
+        constexpr_cxx26 Vector2F& rotate(const value_type radians, const Vector2F& ctr) noexcept {
             return rotate(std::sin(radians), std::cos(radians), ctr);
         }
 
         /** Rotates this vector in place, returns *this */
-        constexpr Vec2f& rotate(const float sin, const float cos, const Vec2f& ctr) noexcept {
-            const float x0 = x - ctr.x;
-            const float y0 = y - ctr.y;
+        constexpr Vector2F& rotate(const value_type sin, const value_type cos, const Vector2F& ctr) noexcept {
+            const value_type x0 = x - ctr.x;
+            const value_type y0 = y - ctr.y;
             x = x0 * cos - y0 * sin + ctr.x;
             y = x0 * sin + y0 * cos + ctr.y;
             return *this;
         }
 
         /** Rotates this vector in place, returns *this */
-        Vec2f& rotate(const float radians) noexcept {
+        constexpr_cxx26 Vector2F& rotate(const value_type radians) noexcept {
             return rotate(std::sin(radians), std::cos(radians));
         }
 
         /** Rotates this vector in place, returns *this */
-        constexpr Vec2f& rotate(const float sin, const float cos) noexcept {
-            const float x0 = x;
+        constexpr Vector2F& rotate(const value_type sin, const value_type cos) noexcept {
+            const value_type x0 = x;
             x = x0 * cos - y * sin;
             y = x0 * sin + y * cos;
             return *this;
@@ -163,33 +186,25 @@ namespace jau::math {
         /**
          * Return the squared length of this vector, a.k.a the squared <i>norm</i> or squared <i>magnitude</i>
          */
-        constexpr float length_sq() const noexcept {
+        constexpr value_type length_sq() const noexcept {
             return x*x + y*y;
         }
 
         /**
          * Return the length of this vector, a.k.a the <i>norm</i> or <i>magnitude</i>
          */
-        constexpr float length() const noexcept {
+        constexpr value_type length() const noexcept {
             return std::sqrt(length_sq());
         }
 
-        /**
-         * Return the direction angle of this vector in radians
-         */
-        float angle() const noexcept {
-            // Utilize atan2 taking y=sin(a) and x=cos(a), resulting in proper direction angle for all quadrants.
-            return std::atan2( y, x );
-        }
-
         /** Normalize this vector in place, returns *this */
-        constexpr Vec2f& normalize() noexcept {
-            const float lengthSq = length_sq();
+        constexpr Vector2F& normalize() noexcept {
+            const value_type lengthSq = length_sq();
             if ( jau::is_zero( lengthSq ) ) {
-                x = 0.0f;
-                y = 0.0f;
+                x = zero;
+                y = zero;
             } else {
-                const float invSqr = 1.0f / std::sqrt(lengthSq);
+                const value_type invSqr = one / std::sqrt(lengthSq);
                 x *= invSqr;
                 y *= invSqr;
             }
@@ -197,30 +212,38 @@ namespace jau::math {
         }
 
         /**
+         * Return the direction angle of this vector in radians
+         */
+        constexpr_cxx26 value_type angle() const noexcept {
+            // Utilize atan2 taking y=sin(a) and x=cos(a), resulting in proper direction angle for all quadrants.
+            return std::atan2( y, x );
+        }
+
+        /**
          * Return the squared distance between this vector and the given one.
          * <p>
          * When comparing the relative distance between two points it is usually sufficient to compare the squared
          * distances, thus avoiding an expensive square root operation.
          * </p>
          */
-        constexpr float dist_sq(const Vec2f& o) const noexcept {
-            const float dx = x - o.x;
-            const float dy = y - o.y;
+        constexpr value_type dist_sq(const Vector2F& o) const noexcept {
+            const value_type dx = x - o.x;
+            const value_type dy = y - o.y;
             return dx*dx + dy*dy;
         }
 
         /**
          * Return the distance between this vector and the given one.
          */
-        constexpr float dist(const Vec2f& o) const noexcept {
+        constexpr value_type dist(const Vector2F& o) const noexcept {
             return std::sqrt(dist_sq(o));
         }
 
         /**
          * Return the dot product of this vector and the given one
-         * @return the dot product as float
+         * @return the dot product as value_type
          */
-        constexpr float dot(const Vec2f& o) const noexcept {
+        constexpr value_type dot(const Vector2F& o) const noexcept {
             return x*o.x + y*o.y;
         }
 
@@ -231,65 +254,95 @@ namespace jau::math {
          *
          * @return the resulting scalar
          */
-        constexpr float cross(const Vec2f& o) const noexcept {
+        constexpr value_type cross(const Vector2F& o) const noexcept {
             return x * o.y - y * o.x;
         }
 
         /**
          * Return the cosines of the angle between two vectors
          */
-        constexpr float cos_angle(const Vec2f& o) const noexcept {
+        constexpr value_type cos_angle(const Vector2F& o) const noexcept {
             return dot(o) / ( length() * o.length() ) ;
         }
 
         /**
          * Return the angle between two vectors in radians
          */
-        float angle(const Vec2f& o) const noexcept {
+        constexpr_cxx26 value_type angle(const Vector2F& o) const noexcept {
             return std::acos( cos_angle(o) );
         }
 
         /**
          * Return the counter-clock-wise (CCW) normal of this vector, i.e. perp(endicular) vector
          */
-        Vec2f normal_ccw() const noexcept {
-            return Vec2f(-y, x);
+        constexpr Vector2F normal_ccw() const noexcept {
+            return Vector2F(-y, x);
         }
 
-        bool intersects(const Vec2f& o) const noexcept {
-            const float eps = std::numeric_limits<float>::epsilon();
+        constexpr_cxx23 bool intersects(const Vector2F& o) const noexcept {
+            const value_type eps = std::numeric_limits<value_type>::epsilon();
             if( std::abs(x-o.x) >= eps || std::abs(y-o.y) >= eps ) {
                 return false;
             }
             return true;
         }
     };
-    typedef Vec2f Point2f;
 
-    constexpr Vec2f operator+(const Vec2f& lhs, const Vec2f& rhs ) noexcept {
-        return Vec2f(lhs) += rhs;
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    constexpr Vector2F<T> operator+(const Vector2F<T>& lhs, const Vector2F<T>& rhs ) noexcept {
+        // Returning a Vector2F<T> object from the returned reference of operator+=()
+        // may hinder copy-elision or "named return value optimization" (NRVO).
+        // return Vector2F<T>(lhs) += rhs;
+
+        // Returning named object allows copy-elision (NRVO),
+        // only one object is created 'on target'.
+        Vector2F<T> r(lhs); r += rhs; return r;
     }
 
-    constexpr Vec2f operator-(const Vec2f& lhs, const Vec2f& rhs ) noexcept {
-        return Vec2f(lhs) -= rhs;
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    constexpr Vector2F<T> operator-(const Vector2F<T>& lhs, const Vector2F<T>& rhs ) noexcept {
+        Vector2F<T> r(lhs); r -= rhs; return r;
     }
 
-    constexpr Vec2f operator*(const Vec2f& lhs, const float s ) noexcept {
-        return Vec2f(lhs) *= s;
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    constexpr Vector2F<T> operator*(const Vector2F<T>& lhs, const T s ) noexcept {
+        Vector2F<T> r(lhs); r *= s; return r;
     }
 
-    constexpr Vec2f operator*(const float s, const Vec2f& rhs) noexcept {
-        return Vec2f(rhs) *= s;
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    constexpr Vector2F<T> operator*(const T s, const Vector2F<T>& rhs) noexcept {
+        Vector2F<T> r(rhs); r *= s; return r;
     }
 
-    constexpr Vec2f operator/(const Vec2f& lhs, const float s ) noexcept {
-        return Vec2f(lhs) /= s;
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    constexpr Vector2F<T> operator/(const Vector2F<T>& lhs, const T s ) noexcept {
+        Vector2F<T> r(lhs); r /= s; return r;
     }
 
-    std::ostream& operator<<(std::ostream& out, const Vec2f& v) noexcept {
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    std::ostream& operator<<(std::ostream& out, const Vector2F<T>& v) noexcept {
         return out << v.toString();
     }
 
+    typedef Vector2F<float> Vec2f;
+    static_assert(alignof(float) == alignof(Vec2f));
+
+    /**
+     * Point2F alias of Vector2F
+     */
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    using Point2F = Vector2F<T>;
+
+    typedef Point2F<float> Point2f;
+    static_assert(alignof(float) == alignof(Point2f));
+
     /**
      * Simple compound denoting a ray.
      * <p>
@@ -301,21 +354,28 @@ namespace jau::math {
      * </pre>
      * </p>
      */
-    class Ray2f {
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    class alignas(T) Ray2F {
     public:
         /** Origin of Ray. */
-        Point2f orig;
+        Point2F<T> orig;
 
         /** Normalized direction vector of ray. */
-        Vec2f dir;
+        Vector2F<T> dir;
 
         std::string toString() const noexcept { return "Ray[orig "+orig.toString()+", dir "+dir.toString() +"]"; }
     };
 
-    std::ostream& operator<<(std::ostream& out, const Ray2f& v) noexcept {
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    std::ostream& operator<<(std::ostream& out, const Ray2F<T>& v) noexcept {
         return out << v.toString();
     }
 
+    typedef Ray2F<float> Ray2f;
+    static_assert(alignof(float) == alignof(Ray2f));
+
     /**@}*/
 
 } // namespace jau::math
diff --git a/include/jau/math/vec2i.hpp b/include/jau/math/vec2i.hpp
index 8b6934c..40f3ecd 100644
--- a/include/jau/math/vec2i.hpp
+++ b/include/jau/math/vec2i.hpp
@@ -27,6 +27,7 @@
 #include <cmath>
 #include <cstdarg>
 #include <cstdint>
+#include <cassert>
 #include <limits>
 #include <string>
 #include <iostream>
@@ -43,32 +44,59 @@ namespace jau::math {
     /**
      * 2D vector using two integer components.
      */
-    struct Vec2i {
-        int x;
-        int y;
+    template<typename Value_type,
+             std::enable_if_t<std::is_integral_v<Value_type>, bool> = true>
+    struct alignas(Value_type) Vector2I {
+        typedef Value_type                  value_type;
+        typedef value_type*                 pointer;
+        typedef const value_type*           const_pointer;
+        typedef value_type&                 reference;
+        typedef const value_type&           const_reference;
+        typedef value_type*                 iterator;
+        typedef const value_type*           const_iterator;
 
-        constexpr Vec2i() noexcept
-        : x(0), y(0) {}
+        typedef typename jau::float_bytes<sizeof(value_type)>::type float_type;
 
-        constexpr Vec2i(const int x_, const int y_) noexcept
+        constexpr static const value_type zero = value_type(0);
+        constexpr static const value_type one  = value_type(1);
+
+        value_type x;
+        value_type y;
+
+        constexpr Vector2I() noexcept
+        : x(zero), y(zero) {}
+
+        constexpr Vector2I(const value_type v) noexcept
+        : x(v), y(v) {}
+
+        constexpr Vector2I(const value_type x_, const value_type y_) noexcept
         : x(x_), y(y_) {}
 
-        constexpr Vec2i(const Vec2i& o) noexcept = default;
-        constexpr Vec2i(Vec2i&& o) noexcept = default;
-        constexpr Vec2i& operator=(const Vec2i&) noexcept = default;
-        constexpr Vec2i& operator=(Vec2i&&) noexcept = default;
+        constexpr Vector2I(const Vector2I& o) noexcept = default;
+        constexpr Vector2I(Vector2I&& o) noexcept = default;
+        constexpr Vector2I& operator=(const Vector2I&) noexcept = default;
+        constexpr Vector2I& operator=(Vector2I&&) noexcept = default;
 
-        /** Returns read-only component w/o boundary check */
-        int operator[](size_t i) const noexcept {
-            return reinterpret_cast<const int*>(this)[i];
+        /** Returns read-only component */
+        constexpr value_type operator[](size_t i) const noexcept {
+            assert(i < 2);
+            return (&x)[i];
         }
 
-        /** Returns writeable reference to component w/o boundary check */
-        int& operator[](size_t i) noexcept {
-            return reinterpret_cast<int*>(this)[i];
+        /** Returns writeable reference to component */
+        constexpr reference operator[](size_t i) noexcept {
+            assert(i < 2);
+            return (&x)[i];
         }
 
-        constexpr bool operator==(const Vec2i& rhs ) const noexcept {
+        /** xy = this, returns xy. */
+        constexpr iterator get(iterator xy) const noexcept {
+            xy[0] = x;
+            xy[1] = y;
+            return xy;
+        }
+
+        constexpr bool operator==(const Vector2I& rhs ) const noexcept {
             if( this == &rhs ) {
                 return true;
             }
@@ -79,18 +107,23 @@ namespace jau::math {
             return ...
         } */
 
-        constexpr Vec2i& set(const int vx, const int vy) noexcept
+        constexpr Vector2I& set(const value_type vx, const value_type vy) noexcept
         { x=vx; y=vy; return *this; }
 
-        constexpr Vec2i& add(const int dx, const int dy) noexcept
+        /** this = xy, returns this. */
+        constexpr Vector2I& set(const_iterator xy) noexcept
+        { x=xy[0]; y=xy[1]; return *this; }
+
+
+        constexpr Vector2I& add(const value_type dx, const value_type dy) noexcept
         { x+=dx; y+=dy; return *this; }
 
-        constexpr Vec2i& operator+=(const Vec2i& rhs ) noexcept {
+        constexpr Vector2I& operator+=(const Vector2I& rhs ) noexcept {
             x+=rhs.x; y+=rhs.y;
             return *this;
         }
 
-        constexpr Vec2i& operator-=(const Vec2i& rhs ) noexcept {
+        constexpr Vector2I& operator-=(const Vector2I& rhs ) noexcept {
             x-=rhs.x; y-=rhs.y;
             return *this;
         }
@@ -100,7 +133,7 @@ namespace jau::math {
          * @param s scale factor
          * @return this instance
          */
-        constexpr Vec2i& operator*=(const int s ) noexcept {
+        constexpr Vector2I& operator*=(const value_type s ) noexcept {
             x*=s; y*=s;
             return *this;
         }
@@ -110,56 +143,110 @@ namespace jau::math {
          * @param s scale factor
          * @return this instance
          */
-        constexpr Vec2i& operator/=(const int s ) noexcept {
+        constexpr Vector2I& operator/=(const value_type s ) noexcept {
             x/=s; y/=s;
             return *this;
         }
 
-        void rotate(const float radians, const Vec2i& ctr) {
-            const float cos = std::cos(radians);
-            const float sin = std::sin(radians);
+        constexpr_cxx26 void rotate(const float_type radians, const Vector2I& ctr) {
+            const float_type cos = std::cos(radians);
+            const float_type sin = std::sin(radians);
             rotate(sin, cos, ctr);
         }
-        void rotate(const float sin, const float cos, const Vec2i& ctr) {
-            const float x0 = (float)(x - ctr.x);
-            const float y0 = (float)(y - ctr.y);
-            const int tmp = jau::round_to_int( x0 * cos - y0 * sin ) + ctr.x;
-                        y = jau::round_to_int( x0 * sin + y0 * cos ) + ctr.y;
+        void rotate(const float_type sin, const float_type cos, const Vector2I& ctr) {
+            const float_type x0 = static_cast<float_type>(x - ctr.x);
+            const float_type y0 = static_cast<float_type>(y - ctr.y);
+            const value_type tmp = jau::round_to_int<float_type>( x0 * cos - y0 * sin ) + ctr.x;
+                         y = jau::round_to_int<float_type>( x0 * sin + y0 * cos ) + ctr.y;
             x = tmp;
         }
 
         std::string toString() const noexcept { return std::to_string(x)+"/"+std::to_string(y); }
 
-        bool intersects(const Vec2i& o) {
+        constexpr bool is_zero() const noexcept {
+            return jau::is_zero(x) && jau::is_zero(y);
+        }
+
+        /**
+         * Return the squared length of this vector, a.k.a the squared <i>norm</i> or squared <i>magnitude</i>
+         */
+        constexpr value_type length_sq() const noexcept {
+            return x*x + y*y;
+        }
+
+        /**
+         * Return the length of this vector, a.k.a the <i>norm</i> or <i>magnitude</i>
+         */
+        constexpr value_type length() const noexcept {
+            // return jau::round_to_int<float_type>( std::sqrt<float_type>( length_sq() ) );
+            return jau::round_to_int<float_type>( std::sqrt( static_cast<float_type> ( length_sq() )) );
+        }
+
+        /** Normalize this vector in place, returns *this */
+        constexpr Vector2I& normalize() noexcept {
+            const value_type lengthSq = length_sq();
+            if ( jau::is_zero( lengthSq ) ) {
+                x = zero;
+                y = zero;
+            } else {
+                const float_type invSqr = one / static_cast<float_type>( std::sqrt(lengthSq) );
+                x = jau::round_to_int<float_type>( x * invSqr );
+                y = jau::round_to_int<float_type>( x * invSqr );
+            }
+            return *this;
+        }
+
+
+        bool intersects(const Vector2I& o) {
             return x == o.x && y == o.y;
         }
     };
-    typedef Vec2i Point2i;
 
-    constexpr Vec2i operator+(const Vec2i& lhs, const Vec2i& rhs ) noexcept {
-        return Vec2i(lhs) += rhs;
+    template<typename T,
+             std::enable_if_t<std::numeric_limits<T>::is_integer, bool> = true>
+    constexpr Vector2I<T> operator+(const Vector2I<T>& lhs, const Vector2I<T>& rhs ) noexcept {
+        // Returning a Vector2I<T> object from the returned reference of operator+=()
+        // may hinder copy-elision or "named return value optimization" (NRVO).
+        // return Vector2I<T>(lhs) += rhs;
+
+        // Returning named object allows copy-elision (NRVO),
+        // only one object is created 'on target'.
+        Vector2I<T> r(lhs); r += rhs; return r;
     }
 
-    constexpr Vec2i operator-(const Vec2i& lhs, const Vec2i& rhs ) noexcept {
-        return Vec2i(lhs) -= rhs;
+    template<typename T,
+             std::enable_if_t<std::numeric_limits<T>::is_integer, bool> = true>
+    constexpr Vector2I<T> operator-(const Vector2I<T>& lhs, const Vector2I<T>& rhs ) noexcept {
+        Vector2I<T> r(lhs); r -= rhs; return r;
     }
 
-    constexpr Vec2i operator*(const Vec2i& lhs, const float s ) noexcept {
-        return Vec2i(lhs) *= s;
+    template<typename T,
+             std::enable_if_t<std::numeric_limits<T>::is_integer, bool> = true>
+    constexpr Vector2I<T> operator*(const Vector2I<T>& lhs, const float s ) noexcept {
+        Vector2I<T> r(lhs); r *= s; return r;
     }
 
-    constexpr Vec2i operator*(const float s, const Vec2i& rhs) noexcept {
-        return Vec2i(rhs) *= s;
+    template<typename T,
+             std::enable_if_t<std::numeric_limits<T>::is_integer, bool> = true>
+    constexpr Vector2I<T> operator*(const float s, const Vector2I<T>& rhs) noexcept {
+        Vector2I<T> r(rhs); r *= s; return r;
     }
 
-    constexpr Vec2i operator/(const Vec2i& lhs, const float s ) noexcept {
-        return Vec2i(lhs) /= s;
+    template<typename T,
+             std::enable_if_t<std::numeric_limits<T>::is_integer, bool> = true>
+    constexpr Vector2I<T> operator/(const Vector2I<T>& lhs, const float s ) noexcept {
+        Vector2I<T> r(lhs); r /= s; return r;
     }
 
-    std::ostream& operator<<(std::ostream& out, const Vec2i& v) noexcept {
+    template<typename T,
+             std::enable_if_t<std::numeric_limits<T>::is_integer, bool> = true>
+    std::ostream& operator<<(std::ostream& out, const Vector2I<T>& v) noexcept {
         return out << v.toString();
     }
 
+    typedef Vector2I<int> Vec2i;
+    static_assert(alignof(int) == alignof(Vec2i));
+
     /**@}*/
 
 } // namespace jau::math
diff --git a/include/jau/math/vec3f.hpp b/include/jau/math/vec3f.hpp
index e94ad40..74a2034 100644
--- a/include/jau/math/vec3f.hpp
+++ b/include/jau/math/vec3f.hpp
@@ -27,8 +27,10 @@
 #include <cmath>
 #include <cstdarg>
 #include <cstdint>
+#include <cassert>
 #include <limits>
 #include <string>
+#include <initializer_list>
 #include <iostream>
 
 #include <jau/float_math.hpp>
@@ -41,47 +43,69 @@ namespace jau::math {
      *  @{
      */
 
-    class Vec4f; // forward
-
     /**
-     * 3D vector using three float components.
+     * 3D vector using three value_type components.
      */
-    struct Vec3f {
+    template<typename Value_type,
+             std::enable_if_t<std::is_floating_point_v<Value_type>, bool> = true>
+    class alignas(Value_type) Vector3F {
         public:
-            float x;
-            float y;
-            float z;
+            typedef Value_type                  value_type;
+            typedef value_type*                 pointer;
+            typedef const value_type*           const_pointer;
+            typedef value_type&                 reference;
+            typedef const value_type&           const_reference;
+            typedef value_type*                 iterator;
+            typedef const value_type*           const_iterator;
+
+            constexpr static const value_type zero = value_type(0);
+            constexpr static const value_type one  = value_type(1);
+
+            value_type x;
+            value_type y;
+            value_type z;
 
-            constexpr Vec3f() noexcept
-            : x(0), y(0), z(0) {}
+            constexpr Vector3F() noexcept
+            : x(zero), y(zero), z(zero) {}
 
-            constexpr Vec3f(const float x_, const float y_, const float z_) noexcept
+            constexpr Vector3F(const value_type v) noexcept
+            : x(v), y(v), z(v) {}
+
+            constexpr Vector3F(const value_type x_, const value_type y_, const value_type z_) noexcept
             : x(x_), y(y_), z(z_) {}
 
-            constexpr Vec3f(const Vec3f& o) noexcept = default;
-            constexpr Vec3f(Vec3f&& o) noexcept = default;
-            constexpr Vec3f& operator=(const Vec3f&) noexcept = default;
-            constexpr Vec3f& operator=(Vec3f&&) noexcept = default;
+            constexpr Vector3F(const_iterator v) noexcept
+            : x(v[0]), y(v[1]), z(v[2]) {}
+
+            constexpr Vector3F(std::initializer_list<value_type> v) noexcept
+            : x(v[0]), y(v[1]), z(v[2]) {}
 
-            /** Returns read-only component w/o boundary check */
-            float operator[](size_t i) const noexcept {
-                return reinterpret_cast<const float*>(this)[i];
+            constexpr Vector3F(const Vector3F& o) noexcept = default;
+            constexpr Vector3F(Vector3F&& o) noexcept = default;
+            constexpr Vector3F& operator=(const Vector3F&) noexcept = default;
+            constexpr Vector3F& operator=(Vector3F&&) noexcept = default;
+
+            /** Returns read-only component */
+            constexpr value_type operator[](size_t i) const noexcept {
+                assert(i < 3);
+                return (&x)[i];
             }
 
-            /** Returns writeable reference to component w/o boundary check */
-            float& operator[](size_t i) noexcept {
-                return reinterpret_cast<float*>(this)[i];
+            /** Returns writeable reference to component  */
+            constexpr reference operator[](size_t i) noexcept {
+                assert(i < 3);
+                return (&x)[i];
             }
 
             /** xyz = this, returns xyz. */
-            float* get(float xyz[/*3*/]) const noexcept {
+            constexpr iterator get(iterator xyz) const noexcept {
                 xyz[0] = x;
                 xyz[1] = y;
                 xyz[2] = z;
                 return xyz;
             }
 
-            constexpr bool equals(const Vec3f& o, const float epsilon=std::numeric_limits<float>::epsilon()) const noexcept {
+            constexpr bool equals(const Vector3F& o, const value_type epsilon=std::numeric_limits<value_type>::epsilon()) const noexcept {
                 if( this == &o ) {
                     return true;
                 } else {
@@ -90,7 +114,7 @@ namespace jau::math {
                            jau::equals(z, o.z, epsilon);
                 }
             }
-            constexpr bool operator==(const Vec3f& rhs) const noexcept {
+            constexpr bool operator==(const Vector3F& rhs) const noexcept {
                 return equals(rhs);
             }
             /** TODO
@@ -99,27 +123,28 @@ namespace jau::math {
             } */
 
             /** this = { o, z }, returns this. */
-            constexpr Vec3f& set(const Vec2f& o, const float z_) noexcept {
+            constexpr Vector3F& set(const Vec2f& o, const value_type z_) noexcept {
                 x = o.x;
                 y = o.y;
                 z = z_;
                 return *this;
             }
 
-            /** this = o while dropping w, returns this. */
-            Vec3f& set(const Vec4f& o) noexcept;
-
-            constexpr Vec3f& set(const float vx, const float vy, const float vz) noexcept
+            constexpr Vector3F& set(const value_type vx, const value_type vy, const value_type vz) noexcept
             { x=vx; y=vy; z=vz; return *this; }
 
-            constexpr void add(const float dx, const float dy, const float dz) noexcept { x+=dx; y+=dy; z+=dz; }
+            /** this = xyz, returns this. */
+            constexpr Vector3F& set(const_iterator xyz) noexcept
+            { x=xyz[0]; y=xyz[1]; z=xyz[2]; return *this; }
 
-            constexpr Vec3f& operator+=(const Vec3f& rhs ) noexcept {
+            constexpr void add(const value_type dx, const value_type dy, const value_type dz) noexcept { x+=dx; y+=dy; z+=dz; }
+
+            constexpr Vector3F& operator+=(const Vector3F& rhs ) noexcept {
                 x+=rhs.x; y+=rhs.y; z+=rhs.z;
                 return *this;
             }
 
-            constexpr Vec3f& operator-=(const Vec3f& rhs ) noexcept {
+            constexpr Vector3F& operator-=(const Vector3F& rhs ) noexcept {
                 x-=rhs.x; y-=rhs.y; z-=rhs.z;
                 return *this;
             }
@@ -129,7 +154,7 @@ namespace jau::math {
              * @param s scale factor
              * @return this instance
              */
-            constexpr Vec3f& operator*=(const float s ) noexcept {
+            constexpr Vector3F& operator*=(const value_type s ) noexcept {
                 x*=s; y*=s; z*=s;
                 return *this;
             }
@@ -139,28 +164,22 @@ namespace jau::math {
              * @param s scale factor
              * @return this instance
              */
-            constexpr Vec3f& operator/=(const float s ) noexcept {
+            constexpr Vector3F& operator/=(const value_type s ) noexcept {
                 x/=s; y/=s; z/=s;
                 return *this;
             }
 
-            void rotate(const float radians, const Vec3f& ctr) noexcept {
-                // FIXME z, consider using a matrix or quaternions
-                const float cos = std::cos(radians);
-                const float sin = std::sin(radians);
-                const float x0 = x - ctr.x;
-                const float y0 = y - ctr.y;
-                const float tmp = x0 * cos - y0 * sin + ctr.x;
-                              y = x0 * sin + y0 * cos + ctr.y;
-                x = tmp;
+            /** Rotates this vector around the Z-axis in place, returns *this */
+            constexpr_cxx26 Vector3F& rotateZ(const value_type radians) noexcept {
+                return rotateZ(std::sin(radians), std::cos(radians));
             }
-            constexpr void rotate(const float sin, const float cos, const Vec3f& ctr) noexcept {
-                // FIXME z, consider using a matrix or quaternions
-                const float x0 = x - ctr.x;
-                const float y0 = y - ctr.y;
-                const float tmp = x0 * cos - y0 * sin + ctr.x;
-                              y = x0 * sin + y0 * cos + ctr.y;
-                x = tmp;
+
+            /** Rotates this vector in place, returns *this */
+            constexpr Vector3F& rotateZ(const value_type sin, const value_type cos) noexcept {
+                const value_type x0 = x;
+                x = x0 * cos - y * sin;
+                y = x0 * sin + y * cos;
+                return *this;
             }
 
             std::string toString() const noexcept { return std::to_string(x)+"/"+std::to_string(y)+"/"+std::to_string(z); }
@@ -172,28 +191,28 @@ namespace jau::math {
             /**
              * Return the squared length of a vector, a.k.a the squared <i>norm</i> or squared <i>magnitude</i>
              */
-            constexpr float length_sq() const noexcept {
+            constexpr value_type length_sq() const noexcept {
                 return x*x + y*y + z*z;
             }
 
             /**
              * Return the length of a vector, a.k.a the <i>norm</i> or <i>magnitude</i>
              */
-            constexpr float length() const noexcept {
+            constexpr value_type length() const noexcept {
                 return std::sqrt(length_sq());
             }
 
             /**
              * Normalize this vector in place
              */
-            constexpr Vec3f& normalize() noexcept {
-                const float lengthSq = length_sq();
+            constexpr Vector3F& normalize() noexcept {
+                const value_type lengthSq = length_sq();
                 if ( jau::is_zero( lengthSq ) ) {
-                    x = 0.0f;
-                    y = 0.0f;
-                    z = 0.0f;
+                    x = zero;
+                    y = zero;
+                    z = zero;
                 } else {
-                    const float invSqr = 1.0f / std::sqrt(lengthSq);
+                    const value_type invSqr = one / std::sqrt(lengthSq);
                     x *= invSqr;
                     y *= invSqr;
                     z *= invSqr;
@@ -208,25 +227,25 @@ namespace jau::math {
              * distances, thus avoiding an expensive square root operation.
              * </p>
              */
-            constexpr float dist_sq(const Vec3f& o) const noexcept {
-                const float dx = x - o.x;
-                const float dy = y - o.y;
-                const float dz = z - o.z;
+            constexpr value_type dist_sq(const Vector3F& o) const noexcept {
+                const value_type dx = x - o.x;
+                const value_type dy = y - o.y;
+                const value_type dz = z - o.z;
                 return dx*dx + dy*dy + dz*dz;
             }
 
             /**
              * Return the distance between this vector and the given one.
              */
-            constexpr float dist(const Vec3f& o) const noexcept {
+            constexpr value_type dist(const Vector3F& o) const noexcept {
                 return std::sqrt(dist_sq(o));
             }
 
             /**
              * Return the dot product of this vector and the given one
-             * @return the dot product as float
+             * @return the dot product as value_type
              */
-            constexpr float dot(const Vec3f& o) const noexcept {
+            constexpr value_type dot(const Vector3F& o) const noexcept {
                 return x*o.x + y*o.y + z*o.z;
             }
 
@@ -234,8 +253,8 @@ namespace jau::math {
              * cross product this x b
              * @return new resulting vector
              */
-            constexpr Vec3f cross(const Vec3f& b) const noexcept {
-                return Vec3f( y * b.z - z * b.y,
+            constexpr Vector3F cross(const Vector3F& b) const noexcept {
+                return Vector3F( y * b.z - z * b.y,
                               z * b.x - x * b.z,
                               x * b.y - y * b.y);
             }
@@ -244,7 +263,7 @@ namespace jau::math {
              * cross product this = a x b, with a, b different from this
              * @return this for chaining
              */
-            constexpr Vec3f& cross(const Vec3f& a, const Vec3f& b) noexcept {
+            constexpr Vector3F& cross(const Vector3F& a, const Vector3F& b) noexcept {
                 x = a.y * b.z - a.z * b.y;
                 y = a.z * b.x - a.x * b.z;
                 z = a.x * b.y - a.y * b.y;
@@ -256,7 +275,7 @@ namespace jau::math {
              * @param vec1 vector 1
              * @param vec2 vector 2
              */
-            constexpr float cos_angle(const Vec3f& o) const noexcept {
+            constexpr value_type cos_angle(const Vector3F& o) const noexcept {
                 return dot(o) / ( length() * o.length() ) ;
             }
 
@@ -265,44 +284,74 @@ namespace jau::math {
              * @param vec1 vector 1
              * @param vec2 vector 2
              */
-            float angle(const Vec3f& o) const noexcept {
+            constexpr_cxx26 value_type angle(const Vector3F& o) const noexcept {
                 return std::acos( cos_angle(o) );
             }
 
-            bool intersects(const Vec3f& o) const noexcept {
-                const float eps = std::numeric_limits<float>::epsilon();
+            constexpr_cxx23 bool intersects(const Vector3F& o) const noexcept {
+                const value_type eps = std::numeric_limits<value_type>::epsilon();
                 if( std::abs(x-o.x) >= eps || std::abs(y-o.y) >= eps || std::abs(z-o.z) >= eps ) {
                     return false;
                 }
                 return true;
             }
     };
-    typedef Vec3f Point3f;
 
-    constexpr Vec3f operator+(const Vec3f& lhs, const Vec3f& rhs ) noexcept {
-        return Vec3f(lhs) += rhs;
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    constexpr Vector3F<T> operator+(const Vector3F<T>& lhs, const Vector3F<T>& rhs ) noexcept {
+        // Returning a Vector3 object from the returned reference of operator+=()
+        // may hinder copy-elision or "named return value optimization" (NRVO).
+        // return Vector3<T>(lhs) += rhs;
+
+        // Returning named object allows copy-elision (NRVO),
+        // only one object is created 'on target'.
+        Vector3F<T> r(lhs); r += rhs; return r;
     }
 
-    constexpr Vec3f operator-(const Vec3f& lhs, const Vec3f& rhs ) noexcept {
-        return Vec3f(lhs) -= rhs;
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    constexpr Vector3F<T> operator-(const Vector3F<T>& lhs, const Vector3F<T>& rhs ) noexcept {
+        Vector3F<T> r(lhs); r -= rhs; return r;
     }
 
-    constexpr Vec3f operator*(const Vec3f& lhs, const float s ) noexcept {
-        return Vec3f(lhs) *= s;
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    constexpr Vector3F<T> operator*(const Vector3F<T>& lhs, const T s ) noexcept {
+        Vector3F<T> r(lhs); r *= s; return r;
     }
 
-    constexpr Vec3f operator*(const float s, const Vec3f& rhs) noexcept {
-        return Vec3f(rhs) *= s;
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    constexpr Vector3F<T> operator*(const T s, const Vector3F<T>& rhs) noexcept {
+        Vector3F<T> r(rhs); r *= s; return r;
     }
 
-    constexpr Vec3f operator/(const Vec3f& lhs, const float s ) noexcept {
-        return Vec3f(lhs) /= s;
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    constexpr Vector3F<T> operator/(const Vector3F<T>& lhs, const T s ) noexcept {
+        Vector3F<T> r(lhs); r /= s; return r;
     }
 
-    std::ostream& operator<<(std::ostream& out, const Vec3f& v) noexcept {
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    std::ostream& operator<<(std::ostream& out, const Vector3F<T>& v) noexcept {
         return out << v.toString();
     }
 
+    typedef Vector3F<float> Vec3f;
+    static_assert(alignof(float) == alignof(Vec3f));
+
+    /**
+     * Point3F alias of Vector3F
+     */
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    using Point3F = Vector3F<T>;
+
+    typedef Point3F<float> Point3f;
+    static_assert(alignof(float) == alignof(Point3f));
+
     /**
      * Simple compound denoting a ray.
      * <p>
@@ -314,20 +363,27 @@ namespace jau::math {
      * </pre>
      * </p>
      */
-    struct Ray3f {
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    struct alignas(T) Ray3F {
         /** Origin of Ray. */
-        Point3f orig;
+        Point3F<T> orig;
 
         /** Normalized direction vector of ray. */
-        Vec3f dir;
+        Vector3F<T> dir;
 
         std::string toString() const noexcept { return "Ray[orig "+orig.toString()+", dir "+dir.toString() +"]"; }
     };
 
-    std::ostream& operator<<(std::ostream& out, const Ray3f& v) noexcept {
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    std::ostream& operator<<(std::ostream& out, const Ray3F<T>& v) noexcept {
         return out << v.toString();
     }
 
+    typedef Ray3F<float> Ray3f;
+    static_assert(alignof(float) == alignof(Ray3f));
+
     /**@}*/
 
 } // namespace jau::math
diff --git a/include/jau/math/vec4f.hpp b/include/jau/math/vec4f.hpp
index ba23056..dc936b3 100644
--- a/include/jau/math/vec4f.hpp
+++ b/include/jau/math/vec4f.hpp
@@ -27,8 +27,10 @@
 #include <cmath>
 #include <cstdarg>
 #include <cstdint>
+#include <cassert>
 #include <limits>
 #include <string>
+#include <initializer_list>
 #include <iostream>
 
 #include <jau/float_math.hpp>
@@ -42,38 +44,64 @@ namespace jau::math {
      */
 
     /**
-     * 4D vector using four float components.
+     * 4D vector using four value_type components.
      */
-    class Vec4f {
+    template<typename Value_type,
+             std::enable_if_t<std::is_floating_point_v<Value_type>, bool> = true>
+    class alignas(Value_type) Vector4F {
         public:
-            float x;
-            float y;
-            float z;
-            float w;
+            typedef Value_type               value_type;
+            typedef value_type*              pointer;
+            typedef const value_type*        const_pointer;
+            typedef value_type&              reference;
+            typedef const value_type&        const_reference;
+            typedef value_type*              iterator;
+            typedef const value_type*        const_iterator;
 
-            constexpr Vec4f() noexcept
-            : x(0), y(0), z(0), w(0) {}
+            typedef Vector3F<value_type, std::is_floating_point_v<Value_type>> Vec3;
 
-            constexpr Vec4f(const float x_, const float y_, const float z_, const float w_) noexcept
+            constexpr static const value_type zero = value_type(0);
+            constexpr static const value_type one  = value_type(1);
+
+            value_type x;
+            value_type y;
+            value_type z;
+            value_type w;
+
+            constexpr Vector4F() noexcept
+            : x(zero), y(zero), z(zero), w(zero) {}
+
+            constexpr Vector4F(const value_type v) noexcept
+            : x(v), y(v), z(v), w(v) {}
+
+            constexpr Vector4F(const value_type x_, const value_type y_, const value_type z_, const value_type w_) noexcept
             : x(x_), y(y_), z(z_), w(w_) {}
 
-            constexpr Vec4f(const Vec4f& o) noexcept = default;
-            constexpr Vec4f(Vec4f&& o) noexcept = default;
-            constexpr Vec4f& operator=(const Vec4f&) noexcept = default;
-            constexpr Vec4f& operator=(Vec4f&&) noexcept = default;
+            constexpr Vector4F(const_iterator v) noexcept
+            : x(v[0]), y(v[1]), z(v[2]), w(v[3]) {}
+
+            constexpr Vector4F(std::initializer_list<value_type> v) noexcept
+            : x(v[0]), y(v[1]), z(v[2]), w(v[3]) {}
 
-            /** Returns read-only component w/o boundary check */
-            float operator[](size_t i) const noexcept {
-                return reinterpret_cast<const float*>(this)[i];
+            constexpr Vector4F(const Vector4F& o) noexcept = default;
+            constexpr Vector4F(Vector4F&& o) noexcept = default;
+            constexpr Vector4F& operator=(const Vector4F&) noexcept = default;
+            constexpr Vector4F& operator=(Vector4F&&) noexcept = default;
+
+            /** Returns read-only component */
+            constexpr value_type operator[](size_t i) const noexcept {
+                assert(i < 4);
+                return (&x)[i];
             }
 
-            /** Returns writeable reference to component w/o boundary check */
-            float& operator[](size_t i) noexcept {
-                return reinterpret_cast<float*>(this)[i];
+            /** Returns writeable reference to component */
+            constexpr reference operator[](size_t i) noexcept {
+                assert(i < 4);
+                return (&x)[i];
             }
 
             /** xyzw = this, returns xyzw. */
-            float* get(float xyzw[/*4*/]) const noexcept {
+            constexpr iterator get(iterator xyzw) const noexcept {
                 xyzw[0] = x;
                 xyzw[1] = y;
                 xyzw[2] = z;
@@ -81,7 +109,15 @@ namespace jau::math {
                 return xyzw;
             }
 
-            constexpr bool operator==(const Vec4f& rhs ) const noexcept {
+            /** out = { this.x, this.y, this.z } dropping w, returns out. */
+            constexpr Vec3& getVec3(Vec3& out) const noexcept {
+                out.x = x;
+                out.y = y;
+                out.z = z;
+                return out;
+            }
+
+            constexpr bool operator==(const Vector4F& rhs ) const noexcept {
                 if( this == &rhs ) {
                     return true;
                 }
@@ -94,27 +130,26 @@ namespace jau::math {
             } */
 
             /** this = { o, w }, returns this. */
-            constexpr Vec4f& set(const Vec3f& o, const float w_) noexcept {
-                x = o.x;
-                y = o.y;
-                z = o.z;
-                w = w_;
-                return *this;
-            }
+            constexpr Vector4F& set(const Vec3f& o, const value_type w_) noexcept
+            { x = o.x; y = o.y; z = o.z; w = w_; return *this; }
 
-            constexpr Vec4f& set(const float vx, const float vy, const float vz, const float vw) noexcept
+            constexpr Vector4F& set(const value_type vx, const value_type vy, const value_type vz, const value_type vw) noexcept
             { x=vx; y=vy; z=vz; w=vw; return *this; }
 
-            constexpr void add(const float dx, const float dy, const float dz, const float dw) noexcept {
+            /** this = xyzw, returns this. */
+            constexpr Vector4F& set(const_iterator xyzw) noexcept
+            { x=xyzw[0]; y=xyzw[1]; z=xyzw[2]; z=xyzw[3]; return *this; }
+
+            constexpr void add(const value_type dx, const value_type dy, const value_type dz, const value_type dw) noexcept {
                 x+=dx; y+=dy; z+=dz; w+=dw;
             }
 
-            constexpr Vec4f& operator+=(const Vec4f& rhs ) noexcept {
+            constexpr Vector4F& operator+=(const Vector4F& rhs ) noexcept {
                 x+=rhs.x; y+=rhs.y; z+=rhs.z; w+=rhs.w;
                 return *this;
             }
 
-            constexpr Vec4f& operator-=(const Vec4f& rhs ) noexcept {
+            constexpr Vector4F& operator-=(const Vector4F& rhs ) noexcept {
                 x-=rhs.x; y-=rhs.y; z-=rhs.z; w-=rhs.w;
                 return *this;
             }
@@ -124,7 +159,7 @@ namespace jau::math {
              * @param s scale factor
              * @return this instance
              */
-            constexpr Vec4f& operator*=(const float s ) noexcept {
+            constexpr Vector4F& operator*=(const value_type s ) noexcept {
                 x*=s; y*=s; z*=s; w*=s;
                 return *this;
             }
@@ -134,7 +169,7 @@ namespace jau::math {
              * @param s scale factor
              * @return this instance
              */
-            constexpr Vec4f& operator/=(const float s ) noexcept {
+            constexpr Vector4F& operator/=(const value_type s ) noexcept {
                 x/=s; y/=s; z/=s; w/=s;
                 return *this;
             }
@@ -148,29 +183,29 @@ namespace jau::math {
             /**
              * Return the squared length of a vector, a.k.a the squared <i>norm</i> or squared <i>magnitude</i>
              */
-            constexpr float length_sq() const noexcept {
+            constexpr value_type length_sq() const noexcept {
                 return x*x + y*y + z*z + w*w;
             }
 
             /**
              * Return the length of a vector, a.k.a the <i>norm</i> or <i>magnitude</i>
              */
-            constexpr float length() const noexcept {
+            constexpr value_type length() const noexcept {
                 return std::sqrt(length_sq());
             }
 
             /**
              * Normalize this vector in place
              */
-            constexpr Vec4f& normalize() noexcept {
-                const float lengthSq = length_sq();
+            constexpr Vector4F& normalize() noexcept {
+                const value_type lengthSq = length_sq();
                 if ( jau::is_zero( lengthSq ) ) {
-                    x = 0.0f;
-                    y = 0.0f;
-                    z = 0.0f;
-                    w = 0.0f;
+                    x = zero;
+                    y = zero;
+                    z = zero;
+                    w = zero;
                 } else {
-                    const float invSqr = 1.0f / std::sqrt(lengthSq);
+                    const value_type invSqr = one / std::sqrt(lengthSq);
                     x *= invSqr;
                     y *= invSqr;
                     z *= invSqr;
@@ -186,23 +221,23 @@ namespace jau::math {
              * distances, thus avoiding an expensive square root operation.
              * </p>
              */
-            constexpr float dist_sq(const Vec4f& o) const noexcept {
-                const float dx = x - o.x;
-                const float dy = y - o.y;
-                const float dz = z - o.z;
-                const float dw = w - o.w;
+            constexpr value_type dist_sq(const Vector4F& o) const noexcept {
+                const value_type dx = x - o.x;
+                const value_type dy = y - o.y;
+                const value_type dz = z - o.z;
+                const value_type dw = w - o.w;
                 return dx*dx + dy*dy + dz*dz + dw*dw;
             }
 
             /**
              * Return the distance between this vector and the given one.
              */
-            constexpr float dist(const Vec4f& o) const noexcept {
+            constexpr value_type dist(const Vector4F& o) const noexcept {
                 return std::sqrt(dist_sq(o));
             }
 
-            bool intersects(const Vec4f& o) const noexcept {
-                const float eps = std::numeric_limits<float>::epsilon();
+            constexpr_cxx23 bool intersects(const Vector4F& o) const noexcept {
+                const value_type eps = std::numeric_limits<value_type>::epsilon();
                 if( std::abs(x-o.x) >= eps || std::abs(y-o.y) >= eps ||
                     std::abs(z-o.z) >= eps || std::abs(w-o.w) >= eps ) {
                     return false;
@@ -210,32 +245,69 @@ namespace jau::math {
                 return true;
             }
     };
-    typedef Vec4f Point4f;
 
-    constexpr Vec4f operator+(const Vec4f& lhs, const Vec4f& rhs ) noexcept {
-        return Vec4f(lhs) += rhs;
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    constexpr Vector4F<T> operator+(const Vector4F<T>& lhs, const Vector4F<T>& rhs ) noexcept {
+        // Returning a Vector4 object from the returned reference of operator+=()
+        // may hinder copy-elision or "named return value optimization" (NRVO).
+        // return Vector4<T>(lhs) += rhs;
+
+        // Returning named object allows copy-elision (NRVO),
+        // only one object is created 'on target'.
+        Vector4F<T> r(lhs); r += rhs; return r;
     }
 
-    constexpr Vec4f operator-(const Vec4f& lhs, const Vec4f& rhs ) noexcept {
-        return Vec4f(lhs) -= rhs;
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    constexpr Vector4F<T> operator-(const Vector4F<T>& lhs, const Vector4F<T>& rhs ) noexcept {
+        Vector4F<T> r(lhs); r -= rhs; return r;
     }
 
-    constexpr Vec4f operator*(const Vec4f& lhs, const float s ) noexcept {
-        return Vec4f(lhs) *= s;
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    constexpr Vector4F<T> operator*(const Vector4F<T>& lhs, const T s ) noexcept {
+        Vector4F<T> r(lhs); r *= s; return r;
     }
 
-    constexpr Vec4f operator*(const float s, const Vec4f& rhs) noexcept {
-        return Vec4f(rhs) *= s;
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    constexpr Vector4F<T> operator*(const T s, const Vector4F<T>& rhs) noexcept {
+        Vector4F<T> r(rhs); r *= s; return r;
     }
 
-    constexpr Vec4f operator/(const Vec4f& lhs, const float s ) noexcept {
-        return Vec4f(lhs) /= s;
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    constexpr Vector4F<T> operator/(const Vector4F<T>& lhs, const T s ) noexcept {
+        Vector4F<T> r(lhs); r /= s; return r;
     }
 
-    std::ostream& operator<<(std::ostream& out, const Vec4f& v) noexcept {
+    /** out = { this.x, this.y, this.z } dropping w, returns out. */
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    constexpr Vector3F<T> to_vec3(const Vector4F<T>& v) noexcept {
+        Vector3F<T> r; v.getVec3(r); return r;
+    }
+
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    std::ostream& operator<<(std::ostream& out, const Vector4F<T>& v) noexcept {
         return out << v.toString();
     }
 
+    typedef Vector4F<float> Vec4f;
+    static_assert(alignof(float) == alignof(Vec4f));
+
+    /**
+     * Point4F alias of Vector4F
+     */
+    template<typename T,
+             std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+    using Point4F = Vector4F<T>;
+
+    typedef Point4F<float> Point4f;
+    static_assert(alignof(float) == alignof(Point4f));
+
     /**@}*/
 
 } // namespace jau::math
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 29f20a1..6eb11f3 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -14,7 +14,6 @@ set (jaulib_LIB_SRCS
   cpuid.cpp
   debug.cpp
   basic_types.cpp
-  math.cpp
   base_codec.cpp
   byte_stream.cpp
   eui48.cpp
diff --git a/src/math.cpp b/src/math.cpp
deleted file mode 100644
index 8f26d10..0000000
--- a/src/math.cpp
+++ /dev/null
@@ -1,181 +0,0 @@
-/*
- * Author: Sven Gothel <[email protected]>
- * Copyright (c) 2014-2024 Gothel Software e.K.
- *
- * Permission is hereby granted, free of charge, to any person obtaining
- * a copy of this software and associated documentation files (the
- * "Software"), to deal in the Software without restriction, including
- * without limitation the rights to use, copy, modify, merge, publish,
- * distribute, sublicense, and/or sell copies of the Software, and to
- * permit persons to whom the Software is furnished to do so, subject to
- * the following conditions:
- *
- * The above copyright notice and this permission notice shall be
- * included in all copies or substantial portions of the Software.
- *
- * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
- * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
- * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
- * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
- * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
- * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
- */
-
-#include <cstdint>
-#include <cinttypes>
-#include <cstring>
-
-#include <ctime>
-
-#include <algorithm>
-
-#include <jau/debug.hpp>
-#include <jau/string_util.hpp>
-#include <jau/math/vec3f.hpp>
-#include <jau/math/mat4f.hpp>
-#include <jau/math/aabbox3f.hpp>
-#include <jau/math/quaternion.hpp>
-#include <jau/math/geom/frustum.hpp>
-#include <jau/math/geom/plane/affine_transform.hpp>
-
-jau::math::Vec3f& jau::math::Vec3f::set(const Vec4f& o) noexcept {
-    x = o.x;
-    y = o.y;
-    z = o.z;
-    return *this;
-}
-
-jau::math::Mat4f& jau::math::Mat4f::setToRotation(const jau::math::Quaternion& q) noexcept {
-    // pre-multiply scaled-reciprocal-magnitude to reduce multiplications
-    const float norm = q.magnitudeSquared();
-    if ( jau::is_zero(norm) ) {
-        // identity matrix -> srecip = 0f
-        loadIdentity();
-        return *this;
-    }
-    float srecip;
-    if ( jau::equals(1.0f, norm) ) {
-        srecip = 2.0f;
-    } else {
-        srecip = 2.0f / norm;
-    }
-
-    const float x = q.x();
-    const float y = q.y();
-    const float z = q.z();
-    const float w = q.w();
-
-    const float xs = srecip * x;
-    const float ys = srecip * y;
-    const float zs = srecip * z;
-
-    const float xx = x  * xs;
-    const float xy = x  * ys;
-    const float xz = x  * zs;
-    const float xw = xs * w;
-    const float yy = y  * ys;
-    const float yz = y  * zs;
-    const float yw = ys * w;
-    const float zz = z  * zs;
-    const float zw = zs * w;
-
-    m00 = 1.0f - ( yy + zz );
-    m01 =        ( xy - zw );
-    m02 =        ( xz + yw );
-    m03 = 0.0f;
-
-    m10 =        ( xy + zw );
-    m11 = 1.0f - ( xx + zz );
-    m12 =        ( yz - xw );
-    m13 = 0.0f;
-
-    m20 =        ( xz - yw );
-    m21 =        ( yz + xw );
-    m22 = 1.0f - ( xx + yy );
-    m23 = 0.0f;
-
-    m30 = m31 = m32 = 0.0f;
-    m33 = 1.0f;
-    return *this;
-}
-
-jau::math::Quaternion& jau::math::Mat4f::getRotation(jau::math::Quaternion& res) const noexcept {
-    return res.setFromMat3(*this);
-}
-
-jau::math::geom::Frustum& jau::math::Mat4f::getFrustum(jau::math::geom::Frustum& frustum) noexcept {
-    return frustum.setFromMat4(*this);
-}
-
-std::string& jau::row_to_string(std::string& sb, const std::string& f,
-                                const float a[],
-                                const jau::nsize_t rows, const jau::nsize_t columns,
-                                const bool rowMajorOrder, const jau::nsize_t row) noexcept {
-  if(rowMajorOrder) {
-      for(jau::nsize_t c=0; c<columns; ++c) {
-          sb.append( jau::format_string(f.c_str(), a[ row*columns + c ] ) );
-          sb.append(", ");
-      }
-  } else {
-      for(jau::nsize_t r=0; r<columns; ++r) {
-          sb.append( jau::format_string(f.c_str(), a[ row + r*rows ] ) );
-          sb.append(", ");
-      }
-  }
-  return sb;
-}
-
-std::string& jau::mat_to_string(std::string& sb, const std::string& rowPrefix, const std::string& f,
-                                const float a[], const jau::nsize_t rows, const jau::nsize_t columns,
-                                const bool rowMajorOrder) noexcept {
-    sb.append(rowPrefix).append("{ ");
-    for(jau::nsize_t i=0; i<rows; ++i) {
-        if( 0 < i ) {
-            sb.append(rowPrefix).append("  ");
-        }
-        row_to_string(sb, f, a, rows, columns, rowMajorOrder, i);
-        sb.append("\n");
-    }
-    sb.append(rowPrefix).append("}").append("\n");
-    return sb;
-}
-
-std::string jau::math::Mat4f::toString(const std::string& rowPrefix, const std::string& f) const noexcept {
-    std::string sb;
-    float tmp[16];
-    get(tmp);
-    return jau::mat_to_string(sb, rowPrefix, f, tmp, 4, 4, false /* rowMajorOrder */); // creates a copy-out!
-}
-
-jau::math::AABBox3f& jau::math::AABBox3f::resize(const jau::math::AABBox3f& newBox, const jau::math::geom::plane::AffineTransform& t, Vec3f& tmpV3) noexcept {
-    /** test low */
-    {
-        const jau::math::Vec3f& newBL = t.transform(newBox.low(), tmpV3);
-        if (newBL.x < m_lo.x) {
-            m_lo.x = newBL.x;
-        }
-        if (newBL.y < m_lo.y) {
-            m_lo.y = newBL.y;
-        }
-        if (newBL.z < m_lo.z) {
-            m_lo.z = newBL.z;
-        }
-    }
-
-    /** test high */
-    {
-        const jau::math::Vec3f& newTR = t.transform(newBox.high(), tmpV3);
-        if (newTR.x > m_hi.x) {
-            m_hi.x = newTR.x;
-        }
-        if (newTR.y > m_hi.y) {
-            m_hi.y = newTR.y;
-        }
-        if (newTR.z > m_hi.z) {
-            m_hi.z = newTR.z;
-        }
-    }
-    computeCenter();
-    return *this;
-}
diff --git a/test/test_math_quaternion.cpp b/test/test_math_quaternion.cpp
index f9410f7..1a57305 100644
--- a/test/test_math_quaternion.cpp
+++ b/test/test_math_quaternion.cpp
@@ -48,7 +48,7 @@ static const float HALF_PI = (float)M_PI_2;
 static const float QUARTER_PI = (float)M_PI_4;
 static const float EPSILON = std::numeric_limits<float>::epsilon();
 
-static const Quaternion QUAT_IDENT = Quaternion(0, 0, 0, 1);
+static const Quat4f QUAT_IDENT = Quat4f(0, 0, 0, 1);
 
 static const Vec3f ZERO       = Vec3f (  0,  0,  0 );
 static const Vec3f ONE        = Vec3f (  1,  1,  1 );
@@ -68,14 +68,14 @@ static const Vec4f ONE_v4     = Vec4f (  1,  1,  1, 0 );
 //
 
 TEST_CASE( "Test 01 Normalize", "[quaternion][linear_algebra][math]" ) {
-    const Quaternion quat(0, 1, 2, 3);
-    const Quaternion quat2 = Quaternion(quat).normalize();
+    const Quat4f quat(0, 1, 2, 3);
+    const Quat4f quat2 = Quat4f(quat).normalize();
     // Assert.assertTrue(Math.abs(1 - quat2.magnitude()) <= MACH_EPSILON);
     REQUIRE( true == jau::equals( 0.0f, std::abs( 1 - quat2.magnitude() ) ) ) ;
 }
 
 TEST_CASE( "Test 02 Rotate Zero Vector", "[quaternion][linear_algebra][math]" ) {
-    Quaternion quat;
+    Quat4f quat;
     const Vec3f rotVec0 = quat.rotateVector(ZERO);
     REQUIRE( ZERO == rotVec0 );
 }
@@ -83,22 +83,22 @@ TEST_CASE( "Test 02 Rotate Zero Vector", "[quaternion][linear_algebra][math]" )
 TEST_CASE( "Test 03 Invert and Conugate", "[quaternion][linear_algebra][math]" ) {
     // inversion check
     {
-        const Quaternion quat0(0, 1, 2, 3);
-        Quaternion quat0Inv = Quaternion(quat0).invert();
+        const Quat4f quat0(0, 1, 2, 3);
+        Quat4f quat0Inv = Quat4f(quat0).invert();
         REQUIRE( quat0 == quat0Inv.invert() );
     }
     // conjugate check
     {
-        const Quaternion quat0(-1, -2, -3, 4);
-        const Quaternion quat0Conj = Quaternion( 1,  2,  3, 4).conjugate();
+        const Quat4f quat0(-1, -2, -3, 4);
+        const Quat4f quat0Conj = Quat4f( 1,  2,  3, 4).conjugate();
         REQUIRE( quat0 == quat0Conj );
     }
 }
 
 TEST_CASE( "Test 04 Dot", "[quaternion][linear_algebra][math]" ) {
-    const Quaternion quat(7, 2, 5, -1);
+    const Quat4f quat(7, 2, 5, -1);
     REQUIRE( 35.0f == quat.dot(3, 1, 2, -2));
-    REQUIRE( -11.0f == quat.dot(Quaternion(-1, 1, -1, 1)));
+    REQUIRE( -11.0f == quat.dot(Quat4f(-1, 1, -1, 1)));
 }
 
 //
@@ -106,8 +106,8 @@ TEST_CASE( "Test 04 Dot", "[quaternion][linear_algebra][math]" ) {
 //
 
 TEST_CASE( "Test 10 Angle Axis", "[quaternion][linear_algebra][math]" ) {
-    Quaternion quat1 = Quaternion().setFromAngleAxis(HALF_PI, Vec3f ( 2, 0, 0 ));
-    Quaternion quat2 = Quaternion().setFromAngleNormalAxis(HALF_PI, Vec3f ( 1, 0, 0 ) );
+    Quat4f quat1 = Quat4f().setFromAngleAxis(HALF_PI, Vec3f ( 2, 0, 0 ));
+    Quat4f quat2 = Quat4f().setFromAngleNormalAxis(HALF_PI, Vec3f ( 1, 0, 0 ) );
 
     REQUIRE( quat2 == quat1 );
     REQUIRE( true == jau::equals(0.0f, 1 - quat2.magnitude()));
@@ -136,10 +136,10 @@ TEST_CASE( "Test 10 Angle Axis", "[quaternion][linear_algebra][math]" ) {
 
 TEST_CASE( "Test 11 From Vec to Vec", "[quaternion][linear_algebra][math]" ) {
     Vec3f vecOut;
-    Quaternion quat;
+    Quat4f quat;
     quat.setFromVectors(UNIT_Z, UNIT_X);
 
-    Quaternion quat2;
+    Quat4f quat2;
     quat2.setFromNormalVectors(UNIT_Z, UNIT_X);
     REQUIRE( quat == quat2 );
 
@@ -148,19 +148,19 @@ TEST_CASE( "Test 11 From Vec to Vec", "[quaternion][linear_algebra][math]" ) {
 
     quat.setFromVectors(UNIT_Z, UNIT_Z_NEG);
     vecOut = quat.rotateVector(UNIT_Z);
-    REQUIRE_THAT( std::abs( UNIT_Z_NEG.dist(vecOut) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+    REQUIRE_THAT( std::abs( UNIT_Z_NEG.dist(vecOut) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
 
     quat.setFromVectors(UNIT_X, UNIT_X_NEG);
     vecOut = quat.rotateVector(UNIT_X);
-    REQUIRE_THAT( std::abs( UNIT_X_NEG.dist(vecOut) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+    REQUIRE_THAT( std::abs( UNIT_X_NEG.dist(vecOut) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
 
     quat.setFromVectors(UNIT_Y, UNIT_Y_NEG);
     vecOut = quat.rotateVector(UNIT_Y);
-    REQUIRE_THAT( std::abs( UNIT_Y_NEG.dist(vecOut) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+    REQUIRE_THAT( std::abs( UNIT_Y_NEG.dist(vecOut) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
 
     quat.setFromVectors(ONE, ONE_NEG);
     vecOut = quat.rotateVector(ONE);
-    REQUIRE_THAT( std::abs( ONE_NEG.dist(vecOut) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+    REQUIRE_THAT( std::abs( ONE_NEG.dist(vecOut) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
 
     quat.setFromVectors(ZERO, ZERO);
     REQUIRE( QUAT_IDENT == quat );
@@ -169,7 +169,7 @@ TEST_CASE( "Test 11 From Vec to Vec", "[quaternion][linear_algebra][math]" ) {
 
 TEST_CASE( "Test 12 From and to Euler Angles", "[quaternion][linear_algebra][math]" ) {
     // Y.Z.X -> X.Y.Z
-    Quaternion quat;
+    Quat4f quat;
     Vec3f angles0Exp( 0, HALF_PI, 0 );
     quat.setFromEuler(angles0Exp);
     REQUIRE_THAT( quat.magnitude(), Catch::Matchers::WithinAbs(1.0f, EPSILON) );
@@ -177,7 +177,7 @@ TEST_CASE( "Test 12 From and to Euler Angles", "[quaternion][linear_algebra][mat
     Vec3f angles0Has = quat.toEuler();
     REQUIRE( angles0Exp == angles0Has );
 
-    Quaternion quat2;
+    Quat4f quat2;
     quat2.setFromEuler(angles0Has);
     REQUIRE( quat == quat2 );
 
@@ -207,7 +207,7 @@ TEST_CASE( "Test 12 From and to Euler Angles", "[quaternion][linear_algebra][mat
 }
 
 TEST_CASE( "Test 13 From Euler Angles and Rotate Vec", "[quaternion][linear_algebra][math]" ) {
-    Quaternion quat;
+    Quat4f quat;
     quat.setFromEuler(0, HALF_PI, 0); // 90 degrees y-axis
     REQUIRE_THAT( quat.magnitude(), Catch::Matchers::WithinAbs(1.0f, EPSILON) );
 
@@ -231,14 +231,15 @@ TEST_CASE( "Test 14 Matrix", "[matrix][quaternion][linear_algebra][math]" ) {
     // Vec4f vecOut4;
     Mat4f mat1;
     Mat4f mat2;
-    Quaternion quat;
+    Quat4f quat1;
+    Quat4f quat2;
 
     //
     // IDENTITY CHECK
     //
     mat1.loadIdentity();
-    quat.set(0, 0, 0, 0);
-    quat.toMatrix(mat2);
+    quat1.set(0, 0, 0, 0);
+    quat1.toMatrix(mat2);
     // mat2 = quat.toMatrix();
     REQUIRE( mat1 == mat2 );
 
@@ -261,11 +262,13 @@ TEST_CASE( "Test 14 Matrix", "[matrix][quaternion][linear_algebra][math]" ) {
     }
     {
         // Validate Matrix via Euler rotation on Quaternion!
-        quat.setFromEuler(a, 0, 0);
+        quat1.setFromEuler(a, 0, 0);
         {
-            quat.toMatrix(mat2);
+            quat1.toMatrix(mat2);
             // mat2 = quat.toMatrix();
             REQUIRE( mat1 == mat2 );
+            quat2.setFromMat(mat1);
+            REQUIRE( quat1 == quat2 );
 
             float mat2_0[16];
             mat2.get(mat2_0);
@@ -274,32 +277,32 @@ TEST_CASE( "Test 14 Matrix", "[matrix][quaternion][linear_algebra][math]" ) {
             REQUIRE( mat2 == mat2c );
             REQUIRE( mat1 == mat2c );
         }
-        vecHas = quat.rotateVector(UNIT_Y);
-        REQUIRE_THAT( std::abs( UNIT_Z.dist(vecHas) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+        vecHas = quat1.rotateVector(UNIT_Y);
+        REQUIRE_THAT( std::abs( UNIT_Z.dist(vecHas) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
+    }
+    {
+        quat1.toMatrix(mat1);
+        quat2.setFromMat(mat1);
+        REQUIRE( quat1 == quat2 );
     }
-    mat1.getRotation(quat);
-    quat.setFromMat3(mat1);
-    vecHas = quat.rotateVector(UNIT_Y);
-    REQUIRE_THAT( std::abs( UNIT_Z.dist(vecHas) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+    vecHas = quat1.rotateVector(UNIT_Y);
+    REQUIRE_THAT( std::abs( UNIT_Z.dist(vecHas) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
 
-    quat.toMatrix(mat2);
+    quat1.toMatrix(mat2);
     REQUIRE( mat1 == mat2 );
 
-    vecHas = quat.rotateVector(ONE_NEG);
+    vecHas = quat1.rotateVector(ONE_NEG);
     {
         // use Vec3f math
         mat2.mulVec3f(ONE_NEG, vecOut3);
-        REQUIRE_THAT( std::abs( vecHas.dist(vecOut3) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+        REQUIRE_THAT( std::abs( vecHas.dist(vecOut3) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
         REQUIRE( vecHas == vecOut3);
     }
     {
         // use Vec4f math
+        (mat2 * ONE_NEG_v4).getVec3(vecOut3);
 
-        // mat2.mulVec4f(ONE_NEG_v4, vecOut4);
-        // vecOut3.set(vecOut4);
-        vecOut3.set( mat2 * ONE_NEG_v4 ); // simpler
-
-        REQUIRE_THAT( std::abs( vecHas.dist(vecOut3) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+        REQUIRE_THAT( std::abs( vecHas.dist(vecOut3) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
         REQUIRE( vecHas == vecOut3);
     }
 
@@ -317,24 +320,22 @@ TEST_CASE( "Test 14 Matrix", "[matrix][quaternion][linear_algebra][math]" ) {
     }
     {
         // Validate Matrix via Euler rotation on Quaternion!
-        quat.setFromEuler(a, 0, 0);
-        quat.toMatrix(mat2);
+        quat1.setFromEuler(a, 0, 0);
+        quat1.toMatrix(mat2);
         REQUIRE(mat1 == mat2);
-        vecHas = quat.rotateVector(UNIT_Y);
-        REQUIRE_THAT( std::abs( UNIT_Y_NEG.dist(vecHas) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+        vecHas = quat1.rotateVector(UNIT_Y);
+        REQUIRE_THAT( std::abs( UNIT_Y_NEG.dist(vecHas) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
     }
-    quat.setFromMat3(mat1);
-    vecHas = quat.rotateVector(UNIT_Y);
-    REQUIRE_THAT( std::abs( UNIT_Y_NEG.dist(vecHas) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+    quat1.setFromMat(mat1);
+    vecHas = quat1.rotateVector(UNIT_Y);
+    REQUIRE_THAT( std::abs( UNIT_Y_NEG.dist(vecHas) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
 
-    quat.toMatrix(mat2);
+    quat1.toMatrix(mat2);
     REQUIRE(mat1 == mat2);
 
-    vecHas = quat.rotateVector(ONE);
-    // mat2.mulVec4f(ONE_v4, vecOut4);
-    // vecOut3.set(vecOut4);
-    vecOut3.set( mat2 * ONE_v4 ); // simpler
-    REQUIRE_THAT( std::abs( vecHas.dist(vecOut3) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+    vecHas = quat1.rotateVector(ONE);
+    (mat2 * ONE_v4).getVec3(vecOut3);
+    REQUIRE_THAT( std::abs( vecHas.dist(vecOut3) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
 
     //
     // 180 degrees rotation on Y
@@ -350,25 +351,23 @@ TEST_CASE( "Test 14 Matrix", "[matrix][quaternion][linear_algebra][math]" ) {
     }
     {
         // Validate Matrix via Euler rotation on Quaternion!
-        quat.setFromEuler(0, a, 0);
-        quat.toMatrix(mat2);
+        quat1.setFromEuler(0, a, 0);
+        quat1.toMatrix(mat2);
         REQUIRE(mat1 == mat2);
 
-        vecHas = quat.rotateVector(UNIT_X);
-        REQUIRE_THAT( std::abs( UNIT_X_NEG.dist(vecHas) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+        vecHas = quat1.rotateVector(UNIT_X);
+        REQUIRE_THAT( std::abs( UNIT_X_NEG.dist(vecHas) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
     }
-    quat.setFromMat3(mat1);
-    vecHas = quat.rotateVector(UNIT_X);
-    REQUIRE_THAT( std::abs( UNIT_X_NEG.dist(vecHas) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+    quat1.setFromMat(mat1);
+    vecHas = quat1.rotateVector(UNIT_X);
+    REQUIRE_THAT( std::abs( UNIT_X_NEG.dist(vecHas) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
 
-    quat.toMatrix(mat2);
+    quat1.toMatrix(mat2);
     REQUIRE(mat1 == mat2);
 
-    vecHas = quat.rotateVector(ONE_NEG);
-    // mat2.mulVec4f(ONE_NEG_v4, vecOut4);
-    // vecOut3.set(vecOut4);
-    vecOut3.set( mat2 * ONE_NEG_v4 ); // simpler
-    REQUIRE_THAT( std::abs( vecHas.dist(vecOut3) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+    vecHas = quat1.rotateVector(ONE_NEG);
+    (mat2 * ONE_NEG_v4).getVec3(vecOut3);
+    REQUIRE_THAT( std::abs( vecHas.dist(vecOut3) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
 
     //
     // 180 degrees rotation on Z
@@ -384,24 +383,23 @@ TEST_CASE( "Test 14 Matrix", "[matrix][quaternion][linear_algebra][math]" ) {
     }
     {
         // Validate Matrix via Euler rotation on Quaternion!
-        quat.setFromEuler(0, 0, a);
-        quat.toMatrix(mat2);
+        quat1.setFromEuler(0, 0, a);
+        quat1.toMatrix(mat2);
         REQUIRE(mat1 == mat2);
-        vecHas = quat.rotateVector(UNIT_X);
-        REQUIRE_THAT( std::abs( UNIT_X_NEG.dist(vecHas) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+        vecHas = quat1.rotateVector(UNIT_X);
+        REQUIRE_THAT( std::abs( UNIT_X_NEG.dist(vecHas) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
     }
-    quat.setFromMat3(mat1);
-    vecHas = quat.rotateVector(UNIT_X);
-    REQUIRE_THAT( std::abs( UNIT_X_NEG.dist(vecHas) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+    quat1.setFromMat(mat1);
+    vecHas = quat1.rotateVector(UNIT_X);
+    REQUIRE_THAT( std::abs( UNIT_X_NEG.dist(vecHas) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
 
-    quat.toMatrix(mat2);
+    quat1.toMatrix(mat2);
     REQUIRE(mat1 == mat2);
 
-    vecHas = quat.rotateVector(ONE);
-    // mat2.mulVec4f(ONE_v4, vecOut4);
-    // vecOut3.set(vecOut4);
-    vecOut3.set( mat2 * ONE_v4 ); // simpler
-    REQUIRE_THAT( std::abs( vecHas.dist(vecOut3) ), Catch::Matchers::WithinAbs(0.0f, Quaternion::ALLOWED_DEVIANCE) );
+    vecHas = quat1.rotateVector(ONE);
+    vecOut3 = to_vec3( mat2 * ONE_v4 );
+    // (mat2 * ONE_v4).getVec3(vecOut3);
+    REQUIRE_THAT( std::abs( vecHas.dist(vecOut3) ), Catch::Matchers::WithinAbs(0.0f, Quat4f::allowed_deviation) );
 
     //
     // Test Matrix-Columns
diff --git a/test/test_math_vec.cpp b/test/test_math_vec.cpp
index 02f4852..641cae6 100644
--- a/test/test_math_vec.cpp
+++ b/test/test_math_vec.cpp
@@ -28,6 +28,8 @@
 
 #include <jau/test/catch2_ext.hpp>
 
+#include <jau/int_math.hpp>
+#include <jau/float_math.hpp>
 #include <jau/math/vec2f.hpp>
 #include <jau/math/vec2i.hpp>
 #include <jau/math/vec3f.hpp>
@@ -39,16 +41,59 @@
 #include <jau/math/mat4f.hpp>
 #include <jau/math/recti.hpp>
 #include <jau/math/math_error.hpp>
+#include <jau/math/util/sstack.hpp>
 
 using namespace jau;
 using namespace jau::math;
+using namespace jau::math::util;
+
+template<class T, class U>
+void test_vec(std::ostream& out, const char* prefix) {
+    out << "Test: " << std::string(prefix) << ", sizeof(U) = " << sizeof(U) << std::endl;
+    T a( U(1) ), b( U(2) );
+    out << "- a: " << a << ", len = " << std::to_string(a.length()) << ", len(normal(a)) = " << std::to_string(T(a).normalize().length()) << std::endl;
+    out << "- b: " << b << ", len = " << std::to_string(b.length()) << ", len(normal(b)) = " << std::to_string(T(b).normalize().length()) << std::endl;
+
+    REQUIRE( T(U(2)) == a * U(2) );
+    REQUIRE( T(U(1)) == b / U(2) );
+
+    a.normalize();
+    b.normalize();
+
+    REQUIRE( jau::equals(U(1), a.length() ) );
+    REQUIRE( jau::equals(U(1), b.length() ) );
+}
+
+template<class T>
+void dump_align_props(std::ostream& out, const char* prefix) {
+    out << std::string(prefix) << "{size " << std::to_string( sizeof(T) ) << ", alignment " << std::to_string( alignof(T) )
+               << " }" << std::endl;
+}
 
 TEST_CASE( "Math Vec Test 00", "[vec][linear_algebra][math]" ) {
+    static_assert(alignof(int) == alignof(Vec2i));
+    static_assert(alignof(float) == alignof(Vec2f));
+    static_assert(alignof(float) == alignof(Vec3f));
+    static_assert(alignof(float) == alignof(Vec4f));
+    static_assert(alignof(float) == alignof(Mat4f));
+
+    dump_align_props<int>(std::cout, "int");
+    dump_align_props<float>(std::cout, "float");
+    dump_align_props<Vec2i>(std::cout, "Vec2i");
+    dump_align_props<Vec2f>(std::cout, "Vec2f");
+    dump_align_props<Vec3f>(std::cout, "Vec3f");
+    dump_align_props<Vec4f>(std::cout, "Vec4f");
+    dump_align_props<Mat4f>(std::cout, "Mat4f");
+
     std::cout << "A v2 " << Vec2f(1, 2) << std::endl;
     std::cout << "A v3 " << Vec3f(1, 2, 3) << std::endl;
     std::cout << "A v4 " << Vec4f(1, 2, 3, 4) << std::endl;
     {
-        const float mf[] = { 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 16.0f };
+        const float mf[] = {  1.0f,  2.0f,  3.0f,  4.0f, // column 0
+                              5.0f,  6.0f,  7.0f,  8.0f, // column 1
+                              9.0f, 10.0f, 11.0f, 12.0f, // column 2
+                             13.0f, 14.0f, 15.0f, 16.0f  // column 3
+                           };
         std::cout << "A mat4 " << Mat4f(mf) << std::endl;
     }
 
@@ -56,6 +101,18 @@ TEST_CASE( "Math Vec Test 00", "[vec][linear_algebra][math]" ) {
     REQUIRE( Vec3f() == Vec3f(0, 0, 0) );
     REQUIRE( 0.0f == Vec2f().length() );
     REQUIRE( 0.0f == Vec3f().length() );
+
+    test_vec<Vector2I<int>, int>(std::cout, "Vector2I<int>");
+    test_vec<Vector2I<long>, long>(std::cout, "Vector2I<long>");
+
+    test_vec<Vector2F<float>, float>(std::cout, "Vector2F<float>");
+    test_vec<Vector2F<double>, double>(std::cout, "Vector2F<double>");
+
+    test_vec<Vector3F<float>, float>(std::cout, "Vector3F<float>");
+    test_vec<Vector3F<double>, double>(std::cout, "Vector3F<double>");
+
+    test_vec<Vector4F<float>, float>(std::cout, "Vector4F<float>");
+    test_vec<Vector4F<double>, double>(std::cout, "Vector4F<double>");
 }
 
 TEST_CASE( "Math Vec Normalize Test 01", "[vec][linear_algebra][math]" ) {