diff --git a/headers/vmath.hpp/vmath_mat_ext.hpp b/headers/vmath.hpp/vmath_mat_ext.hpp
new file mode 100644
index 0000000..9771b31
--- /dev/null
+++ b/headers/vmath.hpp/vmath_mat_ext.hpp
@@ -0,0 +1,120 @@
+/*******************************************************************************
+ * This file is part of the "https://github.com/blackmatov/vmath.hpp"
+ * For conditions of distribution and use, see copyright notice in LICENSE.md
+ * Copyright (C) 2020, by Matvey Cherevko (blackmatov@gmail.com)
+ ******************************************************************************/
+
+#pragma once
+
+#include "vmath_fwd.hpp"
+
+#include "vmath_fun.hpp"
+#include "vmath_mat.hpp"
+#include "vmath_mat_fun.hpp"
+
+namespace vmath_hpp
+{
+    // identity
+
+    template < typename T >
+    mat<T, 4> identity() {
+        return {
+            {1, 0, 0, 0},
+            {0, 1, 0, 0},
+            {0, 0, 1, 0},
+            {0, 0, 0, 1}};
+    }
+
+    // translate
+
+    template < typename T >
+    mat<T, 4> translate(T x, T y, T z) {
+        return {
+            {1, 0, 0, 0},
+            {0, 1, 0, 0},
+            {0, 0, 1, 0},
+            {x, y, z, 1}};
+    }
+
+    template < typename T >
+    mat<T, 4> translate(const vec<T, 3>& xyz) {
+        return translate(xyz.x, xyz.y, xyz.z);
+    }
+
+    template < typename T >
+    mat<T, 4> translate(const mat<T, 4>& m, T x, T y, T z) {
+        return m * translate(x, y, z);
+    }
+
+    template < typename T >
+    mat<T, 4> translate(const mat<T, 4>& m, const vec<T, 3>& xyz) {
+        return m * translate(xyz);
+    }
+
+    // rotate
+
+    template < typename T >
+    mat<T, 4> rotate(T angle, T axis_x, T axis_y, T axis_z) {
+        const T x = axis_x;
+        const T y = axis_y;
+        const T z = axis_z;
+        const T px = x * x;
+        const T py = y * y;
+        const T pz = z * z;
+        const T cs = cos(angle);
+        const T sn = sin(angle);
+        const T ics = T(1) - cs;
+        const T xym = x * y * ics;
+        const T xzm = x * z * ics;
+        const T yzm = y * z * ics;
+        const T xsn = x * sn;
+        const T ysn = y * sn;
+        const T zsn = z * sn;
+        return {
+            px * ics + cs, xym + zsn,     xzm - ysn,     0,
+            xym - zsn,     py * ics + cs, yzm + xsn,     0,
+            xzm + ysn,     yzm - xsn,     pz * ics + cs, 0,
+            0,             0,             0,             1};
+    }
+
+    template < typename T >
+    mat<T, 4> rotate(T angle, const vec<T, 3>& axis) {
+        return rotate(angle, axis.x, axis.y, axis.z);
+    }
+
+    template < typename T >
+    mat<T, 4> rotate(const mat<T, 4>& m, T angle, T axis_x, T axis_y, T axis_z) {
+        return m * rotate(angle, axis_x, axis_y, axis_z);
+    }
+
+    template < typename T >
+    mat<T, 4> rotate(const mat<T, 4>& m, T angle, const vec<T, 3>& axis) {
+        return m * rotate(angle, axis);
+    }
+
+    // scale
+
+    template < typename T >
+    mat<T, 4> scale(T x, T y, T z) {
+        return {
+            {x, 0, 0, 0},
+            {0, y, 0, 0},
+            {0, 0, z, 0},
+            {0, 0, 0, 1}};
+    }
+
+    template < typename T >
+    mat<T, 4> scale(const vec<T, 3>& xyz) {
+        return scale(xyz.x, xyz.y, xyz.z);
+    }
+
+    template < typename T >
+    mat<T, 4> scale(const mat<T, 4>& m, T x, T y, T z) {
+        return m * scale(x, y, z);
+    }
+
+    template < typename T >
+    mat<T, 4> scale(const mat<T, 4>& m, const vec<T, 3>& xyz) {
+        return m * scale(xyz);
+    }
+}
diff --git a/untests/vmath_mat_ext_tests.cpp b/untests/vmath_mat_ext_tests.cpp
new file mode 100644
index 0000000..3a3ceb0
--- /dev/null
+++ b/untests/vmath_mat_ext_tests.cpp
@@ -0,0 +1,103 @@
+/*******************************************************************************
+ * This file is part of the "https://github.com/blackmatov/vmath.hpp"
+ * For conditions of distribution and use, see copyright notice in LICENSE.md
+ * Copyright (C) 2020, by Matvey Cherevko (blackmatov@gmail.com)
+ ******************************************************************************/
+
+#include <vmath.hpp/vmath_mat_ext.hpp>
+
+#define CATCH_CONFIG_FAST_COMPILE
+#include <catch2/catch.hpp>
+
+namespace
+{
+    using namespace vmath_hpp;
+
+    template < typename T >
+    class approx2 {
+    public:
+        constexpr explicit approx2(T v) : value_(v) {}
+        constexpr explicit approx2(T x, T y) : value_(x, y) {}
+
+        friend constexpr bool operator==(const vec<T, 2>& l, const approx2& r) {
+            return (r.value_.x < l.x + epsilon)
+                && (l.x < r.value_.x + epsilon)
+                && (r.value_.y < l.y + epsilon)
+                && (l.y < r.value_.y + epsilon);
+        }
+    private:
+        vec<T, 2> value_;
+        static constexpr T epsilon = std::numeric_limits<T>::epsilon() * 100;
+    };
+
+    template < typename T >
+    class approx3 {
+    public:
+        constexpr explicit approx3(T v) : value_(v) {}
+        constexpr explicit approx3(T x, T y, T z) : value_(x, y, z) {}
+
+        friend constexpr bool operator==(const vec<T, 3>& l, const approx3& r) {
+            return (r.value_.x < l.x + epsilon)
+                && (l.x < r.value_.x + epsilon)
+                && (r.value_.y < l.y + epsilon)
+                && (l.y < r.value_.y + epsilon)
+                && (r.value_.z < l.z + epsilon)
+                && (l.z < r.value_.z + epsilon);
+        }
+    private:
+        vec<T, 3> value_;
+        static constexpr T epsilon = std::numeric_limits<T>::epsilon() * 100;
+    };
+
+    template < typename T >
+    class approx4 {
+    public:
+        constexpr explicit approx4(T v) : value_(v) {}
+        constexpr explicit approx4(T x, T y, T z, T w) : value_(x, y, z, w) {}
+
+        friend constexpr bool operator==(const vec<T, 4>& l, const approx4& r) {
+            return (r.value_.x < l.x + epsilon)
+                && (l.x < r.value_.x + epsilon)
+                && (r.value_.y < l.y + epsilon)
+                && (l.y < r.value_.y + epsilon)
+                && (r.value_.z < l.z + epsilon)
+                && (l.z < r.value_.z + epsilon)
+                && (r.value_.w < l.w + epsilon)
+                && (l.w < r.value_.w + epsilon);
+        }
+    private:
+        vec<T, 4> value_;
+        static constexpr T epsilon = std::numeric_limits<T>::epsilon() * 100;
+    };
+}
+
+TEST_CASE("vmath/mat_ext") {
+    SECTION("identity") {
+        REQUIRE(vec4f(2.f,3.f,4.f,1.f) * identity<float>() == approx4(2.f,3.f,4.f,1.f));
+        REQUIRE(vec4f(2.f,3.f,4.f,1.f) * identity<float>() == approx4(2.f,3.f,4.f,1.f));
+    }
+
+    SECTION("translate") {
+        REQUIRE(vec4f(2.f,3.f,4.f,1.f) * translate(1.f,2.f,3.f) == approx4(3.f,5.f,7.f,1.f));
+        REQUIRE(vec4f(2.f,3.f,4.f,1.f) * translate(vec3f{1.f,2.f,3.f}) == approx4(3.f,5.f,7.f,1.f));
+
+        REQUIRE(vec4f(2.f,3.f,4.f,1.f) * translate(translate(1.f,2.f,3.f), 1.f,2.f,3.f) == approx4(4.f,7.f,10.f,1.f));
+        REQUIRE(vec4f(2.f,3.f,4.f,1.f) * translate(translate(1.f,2.f,3.f), vec3f{1.f,2.f,3.f}) == approx4(4.f,7.f,10.f,1.f));
+    }
+
+    SECTION("rotate") {
+        REQUIRE(vec4f(2.f,3.f,4.f,1.f) * rotate(float(M_PI),0.f,0.f,1.f) == approx4(-2.f,-3.f,4.f,1.f));
+        REQUIRE(vec4f(2.f,3.f,4.f,1.f) * rotate(float(M_PI),vec3f{0.f,0.f,1.f}) == approx4(-2.f,-3.f,4.f,1.f));
+
+        REQUIRE(vec4f(2.f,3.f,4.f,1.f) * rotate(rotate(float(M_PI_2),0.f,0.f,1.f),float(M_PI_2),0.f,0.f,1.f) == approx4(-2.f,-3.f,4.f,1.f));
+        REQUIRE(vec4f(2.f,3.f,4.f,1.f) * rotate(rotate(float(M_PI_2),0.f,0.f,1.f),float(M_PI_2),vec3f{0.f,0.f,1.f}) == approx4(-2.f,-3.f,4.f,1.f));
+    }
+
+    SECTION("scale") {
+        REQUIRE(vec4f(2.f,3.f,4.f,1.f) * scale(2.f,3.f,4.f) == approx4(4.f,9.f,16.f,1.f));
+        REQUIRE(vec4f(2.f,3.f,4.f,1.f) * scale(vec3f{2.f,3.f,4.f}) == approx4(4.f,9.f,16.f,1.f));
+
+        REQUIRE(vec4f(2.f,3.f,4.f,1.f) * scale(scale(2.f,2.f,2.f), 2.f,3.f,4.f) == approx4(8.f,18.f,32.f,1.f));
+        REQUIRE(vec4f(2.f,3.f,4.f,1.f) * scale(scale(2.f,2.f,2.f), vec3f{2.f,3.f,4.f}) == approx4(8.f,18.f,32.f,1.f));
+    }
+}