diff --git a/ZEngine/ZEngine/Core/Maths/Matrix.h b/ZEngine/ZEngine/Core/Maths/Matrix.h
index 5661722f..ec67cf8d 100644
--- a/ZEngine/ZEngine/Core/Maths/Matrix.h
+++ b/ZEngine/ZEngine/Core/Maths/Matrix.h
@@ -132,7 +132,7 @@ namespace ZEngine::Core::Maths
 
         Mat2()
         {
-            Helpers::secure_memset(m_data, 0, sizeof(m_data));
+            Helpers::secure_memset(m_data, 0, sizeof(m_data), sizeof(m_data));
         }
 
         Mat2(T m00, T m01, T m10, T m11)
@@ -192,7 +192,7 @@ namespace ZEngine::Core::Maths
 
         Mat3()
         {
-            Helpers::secure_memset(m_data, 0, sizeof(m_data));
+            Helpers::secure_memset(m_data, 0, sizeof(m_data), sizeof(m_data));
         }
 
         Mat3(T m00, T m01, T m02, T m10, T m11, T m12, T m20, T m21, T m22)
@@ -270,7 +270,7 @@ namespace ZEngine::Core::Maths
 
             float          invDet = 1.0f / dot(r2, c);
 
-            return Mat<T>(r0.x * invDet, r0.y * invDet, r0.z * invDet, r1.x * invDet, r1.y * invDet, r1.z * invDet, r2.x * invDet, r2.y * invDet, r2.z * invDet);
+            return Mat3<T>(r0.x * invDet, r0.y * invDet, r0.z * invDet, r1.x * invDet, r1.y * invDet, r1.z * invDet, r2.x * invDet, r2.y * invDet, r2.z * invDet);
         };
     };
 
@@ -280,29 +280,29 @@ namespace ZEngine::Core::Maths
         using Matrix<T, 4, 4>::m_data;
         Mat4()
         {
-            Helpers::secure_memset(m_data, 0, sizeof(m_data));
+            Helpers::secure_memset(m_data, 0, sizeof(m_data), sizeof(m_data));
         }
 
         Mat4(T m00, T m01, T m02, T m03, T m10, T m11, T m12, T m13, T m20, T m21, T m22, T m23, T m30, T m31, T m32, T m33)
         {
             m_data[0][0] = m00;
-            m_data[1][0] = m10;
-            m_data[2][0] = m20;
-            m_data[3][0] = m30;
+            m_data[1][0] = m01;
+            m_data[2][0] = m02;
+            m_data[3][0] = m03;
 
-            m_data[0][1] = m01;
+            m_data[0][1] = m10;
             m_data[1][1] = m11;
-            m_data[2][1] = m21;
-            m_data[3][1] = m31;
+            m_data[2][1] = m12;
+            m_data[3][1] = m13;
 
-            m_data[0][2] = m02;
-            m_data[1][2] = m12;
+            m_data[0][2] = m20;
+            m_data[1][2] = m21;
             m_data[2][2] = m22;
-            m_data[3][2] = m32;
+            m_data[3][2] = m23;
 
-            m_data[0][3] = m03;
-            m_data[1][3] = m13;
-            m_data[2][3] = m23;
+            m_data[0][3] = m30;
+            m_data[1][3] = m31;
+            m_data[2][3] = m32;
             m_data[3][3] = m33;
         }
 
@@ -331,35 +331,38 @@ namespace ZEngine::Core::Maths
 
         Mat4<T> Inverse() const
         {
-            auto&        M = *this;
+            auto&    M = *this;
 
-            Vec3<T>      a(M(0, 0), M(1, 0), M(2, 0));
-            Vec3<T>      b(M(0, 1), M(1, 1), M(2, 1));
-            Vec3<T>      c(M(0, 2), M(1, 2), M(2, 2));
-            Vec3<T>      d(M(0, 3), M(1, 3), M(2, 3));
+            Vec3<T>  a(M(0, 0), M(0, 1), M(0, 2));
+            Vec3<T>  b(M(1, 0), M(1, 1), M(1, 2));
+            Vec3<T>  c(M(2, 0), M(2, 1), M(2, 2));
+            Vec3<T>  d(M(3, 0), M(3, 1), M(3, 2));
 
-            const float& x       = M(0, 3);
-            const float& y       = M(1, 3);
-            const float& z       = M(2, 3);
-            const float& w       = M(3, 3);
+            const T& x   = M(0, 3);
+            const T& y   = M(1, 3);
+            const T& z   = M(2, 3);
+            const T& w   = M(3, 3);
 
-            Vec3<T>      s       = cross3d(a, b);
-            Vec3<T>      t       = cross3d(c, d);
-            Vec3<T>      u       = a * y - b * x;
-            Vec3<T>      v       = c * w - d * z;
+            Vec3<T>  s   = cross3d(a, b);
+            Vec3<T>  t   = cross3d(c, d);
+            Vec3<T>  u   = a * y - b * x;
+            Vec3<T>  v   = c * w - d * z;
 
-            float        invDet  = T(1) / (dot(s, v) + dot(t, u));
-            s                   *= invDet;
-            t                   *= invDet;
-            u                   *= invDet;
-            v                   *= invDet;
+            T        det = dot(s, v) + dot(t, u);
+            ZENGINE_VALIDATE_ASSERT(det != T(0), "Matrix is singular and cannot be inverted");
+
+            T invDet    = T(1) / det;
+            s          *= invDet;
+            t          *= invDet;
+            u          *= invDet;
+            v          *= invDet;
 
-            Vec3<T> r0           = cross3d(b, v) + (t * y);
-            Vec3<T> r1           = cross3d(v, a) - (t * x);
-            Vec3<T> r2           = cross3d(d, u) + (s * w);
-            Vec3<T> r3           = cross3d(u, c) - (s * z);
+            Vec3<T> r0  = cross3d(b, v) + (t * y);
+            Vec3<T> r1  = cross3d(v, a) - (t * x);
+            Vec3<T> r2  = cross3d(d, u) + (s * w);
+            Vec3<T> r3  = cross3d(u, c) - (s * z);
 
-            return Mat4<T>(r0.x, r0.y, r0.z, -dot(b, t), r1.x, r1.y, r1.z, dot(a, t), r2.x, r2.y, r2.z, -dot(d, s), r3.x, r3.y, r3.z, dot(c, s));
+            return Mat4<T>(r0.x, r1.x, r2.x, r3.x, r0.y, r1.y, r2.y, r3.y, r0.z, r1.z, r2.z, r3.z, -dot(r0, d), -dot(r1, d), -dot(r2, a), -dot(r3, a));
         }
     };
 
@@ -434,4 +437,29 @@ namespace ZEngine::Core::Maths
 
         return m00 * (m11 * (m22 * m33 - m23 * m32) - m12 * (m21 * m33 - m23 * m31) + m13 * (m21 * m32 - m22 * m31)) - m01 * (m10 * (m22 * m33 - m23 * m32) - m12 * (m20 * m33 - m23 * m30) + m13 * (m20 * m32 - m22 * m30)) + m02 * (m10 * (m21 * m33 - m23 * m31) - m11 * (m20 * m33 - m23 * m30) + m13 * (m20 * m31 - m21 * m30)) - m03 * (m10 * (m21 * m32 - m22 * m31) - m11 * (m20 * m32 - m22 * m30) + m12 * (m20 * m31 - m21 * m30));
     }
+
+    template <typename T>
+    inline Mat4<T> operator*(const Mat4<T>& a, const Mat4<T>& b)
+    {
+        return (Mat4<T>(
+            a(0, 0) * b(0, 0) + a(0, 1) * b(1, 0) + a(0, 2) * b(2, 0) + a(0, 3) * b(3, 0),
+            a(0, 0) * b(0, 1) + a(0, 1) * b(1, 1) + a(0, 2) * b(2, 1) + a(0, 3) * b(3, 1),
+            a(0, 0) * b(0, 2) + a(0, 1) * b(1, 2) + a(0, 2) * b(2, 2) + a(0, 3) * b(3, 2),
+            a(0, 0) * b(0, 3) + a(0, 1) * b(1, 3) + a(0, 2) * b(2, 3) + a(0, 3) * b(3, 3),
+
+            a(1, 0) * b(0, 0) + a(1, 1) * b(1, 0) + a(1, 2) * b(2, 0) + a(1, 3) * b(3, 0),
+            a(1, 0) * b(0, 1) + a(1, 1) * b(1, 1) + a(1, 2) * b(2, 1) + a(1, 3) * b(3, 1),
+            a(1, 0) * b(0, 2) + a(1, 1) * b(1, 2) + a(1, 2) * b(2, 2) + a(1, 3) * b(3, 2),
+            a(1, 0) * b(0, 3) + a(1, 1) * b(1, 3) + a(1, 2) * b(2, 3) + a(1, 3) * b(3, 3),
+
+            a(2, 0) * b(0, 0) + a(2, 1) * b(1, 0) + a(2, 2) * b(2, 0) + a(2, 3) * b(3, 0),
+            a(2, 0) * b(0, 1) + a(2, 1) * b(1, 1) + a(2, 2) * b(2, 1) + a(2, 3) * b(3, 1),
+            a(2, 0) * b(0, 2) + a(2, 1) * b(1, 2) + a(2, 2) * b(2, 2) + a(2, 3) * b(3, 2),
+            a(2, 0) * b(0, 3) + a(2, 1) * b(1, 3) + a(2, 2) * b(2, 3) + a(2, 3) * b(3, 3),
+
+            a(3, 0) * b(0, 0) + a(3, 1) * b(1, 0) + a(3, 2) * b(2, 0) + a(3, 3) * b(3, 0),
+            a(3, 0) * b(0, 1) + a(3, 1) * b(1, 1) + a(3, 2) * b(2, 1) + a(3, 3) * b(3, 1),
+            a(3, 0) * b(0, 2) + a(3, 1) * b(1, 2) + a(3, 2) * b(2, 2) + a(3, 3) * b(3, 2),
+            a(3, 0) * b(0, 3) + a(3, 1) * b(1, 3) + a(3, 2) * b(2, 3) + a(3, 3) * b(3, 3)));
+    }
 } // namespace ZEngine::Core::Maths
diff --git a/ZEngine/tests/Maths/Matrix_test.cpp b/ZEngine/tests/Maths/Matrix_test.cpp
index 6a9d46ef..e57e6fde 100644
--- a/ZEngine/tests/Maths/Matrix_test.cpp
+++ b/ZEngine/tests/Maths/Matrix_test.cpp
@@ -5,6 +5,381 @@ using namespace ZEngine::Core::Maths;
 
 constexpr float EPSILON = 1e-5f;
 
+// ========== CONSTRUCTION TESTS ==========
+
+TEST(MatrixTest, DefaultConstruction)
+{
+    Mat2f m2;
+    Mat3f m3;
+    Mat4f m4;
+
+    // Default construction should zero-initialize
+    for (size_t i = 0; i < 2; ++i)
+    {
+        for (size_t j = 0; j < 2; ++j)
+        {
+            EXPECT_NEAR(m2(i, j), 0.0f, EPSILON);
+        }
+    }
+
+    for (size_t i = 0; i < 3; ++i)
+    {
+        for (size_t j = 0; j < 3; ++j)
+        {
+            EXPECT_NEAR(m3(i, j), 0.0f, EPSILON);
+        }
+    }
+
+    for (size_t i = 0; i < 4; ++i)
+    {
+        for (size_t j = 0; j < 4; ++j)
+        {
+            EXPECT_NEAR(m4(i, j), 0.0f, EPSILON);
+        }
+    }
+}
+
+TEST(MatrixTest, ParameterConstruction)
+{
+    // Mat2 construction
+    Mat2f m2(1.1f, 2.2f, 3.3f, 4.4f);
+    EXPECT_NEAR(m2(0, 0), 1.1f, EPSILON);
+    EXPECT_NEAR(m2(0, 1), 2.2f, EPSILON);
+    EXPECT_NEAR(m2(1, 0), 3.3f, EPSILON);
+    EXPECT_NEAR(m2(1, 1), 4.4f, EPSILON);
+
+    // Mat3 construction
+    Mat3f m3(1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f);
+    EXPECT_NEAR(m3(0, 0), 1.0f, EPSILON);
+    EXPECT_NEAR(m3(0, 1), 2.0f, EPSILON);
+    EXPECT_NEAR(m3(0, 2), 3.0f, EPSILON);
+    EXPECT_NEAR(m3(1, 0), 4.0f, EPSILON);
+    EXPECT_NEAR(m3(1, 1), 5.0f, EPSILON);
+    EXPECT_NEAR(m3(1, 2), 6.0f, EPSILON);
+    EXPECT_NEAR(m3(2, 0), 7.0f, EPSILON);
+    EXPECT_NEAR(m3(2, 1), 8.0f, EPSILON);
+    EXPECT_NEAR(m3(2, 2), 9.0f, EPSILON);
+
+    // Mat4 construction
+    Mat4f m4(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16);
+    EXPECT_NEAR(m4(0, 0), 1.0f, EPSILON);
+    EXPECT_NEAR(m4(1, 0), 5.0f, EPSILON);
+    EXPECT_NEAR(m4(3, 3), 16.0f, EPSILON);
+}
+
+TEST(MatrixTest, VectorConstruction)
+{
+    Vec2<float> v1(1.0f, 2.0f);
+    Vec2<float> v2(3.0f, 4.0f);
+    Mat2f       m2(v1, v2);
+
+    EXPECT_NEAR(m2(0, 0), 1.0f, EPSILON);
+    EXPECT_NEAR(m2(1, 0), 2.0f, EPSILON);
+    EXPECT_NEAR(m2(0, 1), 3.0f, EPSILON);
+    EXPECT_NEAR(m2(1, 1), 4.0f, EPSILON);
+
+    Vec3<float> v3a(1.0f, 2.0f, 3.0f);
+    Vec3<float> v3b(4.0f, 5.0f, 6.0f);
+    Vec3<float> v3c(7.0f, 8.0f, 9.0f);
+    Mat3f       m3(v3a, v3b, v3c);
+
+    EXPECT_NEAR(m3(0, 0), 1.0f, EPSILON);
+    EXPECT_NEAR(m3(1, 0), 2.0f, EPSILON);
+    EXPECT_NEAR(m3(2, 0), 3.0f, EPSILON);
+}
+
+TEST(MatrixTest, CopyConstruction)
+{
+    Mat2f original(1.0f, 2.0f, 3.0f, 4.0f);
+    Mat2f copy(original);
+
+    for (size_t i = 0; i < 2; ++i)
+    {
+        for (size_t j = 0; j < 2; ++j)
+        {
+            EXPECT_NEAR(copy(i, j), original(i, j), EPSILON);
+        }
+    }
+}
+
+// ========== INDEXING TESTS ==========
+
+TEST(MatrixTest, IndexingOperators)
+{
+    Mat3f        m3(1, 2, 3, 4, 5, 6, 7, 8, 9);
+
+    // Test column access operator for Mat3
+    Vec3<float>& col0 = m3[0];
+    EXPECT_NEAR(col0.x, 1.0f, EPSILON);
+    EXPECT_NEAR(col0.y, 4.0f, EPSILON);
+}
+
+// ========== ARITHMETIC TESTS ==========
+
+TEST(MatrixTest, Addition)
+{
+    Mat2f m1(1, 2, 3, 4);
+    Mat2f m2(5, 6, 7, 8);
+    Mat2f result = m1 + m2;
+
+    EXPECT_NEAR(result(0, 0), 6.0f, EPSILON);
+    EXPECT_NEAR(result(0, 1), 8.0f, EPSILON);
+    EXPECT_NEAR(result(1, 0), 10.0f, EPSILON);
+    EXPECT_NEAR(result(1, 1), 12.0f, EPSILON);
+
+    // Test compound assignment
+    Mat2f m3(1, 2, 3, 4);
+    m3 += m2;
+    EXPECT_NEAR(m3(0, 0), 6.0f, EPSILON);
+    EXPECT_NEAR(m3(1, 1), 12.0f, EPSILON);
+}
+
+TEST(MatrixTest, Subtraction)
+{
+    Mat2f m1(5, 6, 7, 8);
+    Mat2f m2(1, 2, 3, 4);
+    Mat2f result = m1 - m2;
+
+    EXPECT_NEAR(result(0, 0), 4.0f, EPSILON);
+    EXPECT_NEAR(result(0, 1), 4.0f, EPSILON);
+    EXPECT_NEAR(result(1, 0), 4.0f, EPSILON);
+    EXPECT_NEAR(result(1, 1), 4.0f, EPSILON);
+
+    Mat2f m3(5, 6, 7, 8);
+    m3 -= m2;
+    EXPECT_NEAR(m3(0, 0), 4.0f, EPSILON);
+    EXPECT_NEAR(m3(1, 1), 4.0f, EPSILON);
+}
+
+TEST(MatrixTest, ScalarMultiplication)
+{
+    Mat2f m(1, 2, 3, 4);
+    Mat2f result = m * 2.0f;
+
+    EXPECT_NEAR(result(0, 0), 2.0f, EPSILON);
+    EXPECT_NEAR(result(0, 1), 4.0f, EPSILON);
+    EXPECT_NEAR(result(1, 0), 6.0f, EPSILON);
+    EXPECT_NEAR(result(1, 1), 8.0f, EPSILON);
+
+    Mat2f m2(1, 2, 3, 4);
+    m2 *= 3.0f;
+    EXPECT_NEAR(m2(0, 0), 3.0f, EPSILON);
+    EXPECT_NEAR(m2(1, 1), 12.0f, EPSILON);
+}
+
+TEST(MatrixTest, ScalarDivision)
+{
+    Mat2f m(2, 4, 6, 8);
+    Mat2f result = m / 2.0f;
+
+    EXPECT_NEAR(result(0, 0), 1.0f, EPSILON);
+    EXPECT_NEAR(result(0, 1), 2.0f, EPSILON);
+    EXPECT_NEAR(result(1, 0), 3.0f, EPSILON);
+    EXPECT_NEAR(result(1, 1), 4.0f, EPSILON);
+
+    Mat2f m2(4, 8, 12, 16);
+    m2 /= 4.0f;
+    EXPECT_NEAR(m2(0, 0), 1.0f, EPSILON);
+    EXPECT_NEAR(m2(1, 1), 4.0f, EPSILON);
+}
+
+TEST(MatrixTest, MatrixMultiplication)
+{
+    Mat3f m1(1, 2, 3, 0, 1, 4, 5, 6, 0);
+    Mat3f m2     = Identity<Mat3f>();
+
+    Mat3f result = m1 * m2;
+
+    // Multiplying by identity should result in the original matrix
+    for (size_t i = 0; i < 3; ++i)
+    {
+        for (size_t j = 0; j < 3; ++j)
+        {
+            EXPECT_NEAR(result(i, j), m1(i, j), EPSILON);
+        }
+    }
+
+    // Test specific multiplication
+    Mat3f a = Identity<Mat3f>();
+    Mat3f b(2, 0, 0, 0, 2, 0, 0, 0, 2);
+    Mat3f product = a * b;
+
+    EXPECT_NEAR(product(0, 0), 2.0f, EPSILON);
+    EXPECT_NEAR(product(1, 1), 2.0f, EPSILON);
+    EXPECT_NEAR(product(2, 2), 2.0f, EPSILON);
+}
+
+// ========== DETERMINANT TESTS ==========
+
+TEST(MatrixTest, Determinant2x2)
+{
+    Mat2f m1(1, 2, 3, 4);
+    EXPECT_NEAR(Determinant(m1), -2.0f, EPSILON);
+
+    Mat2f m2(2, 0, 0, 2);
+    EXPECT_NEAR(Determinant(m2), 4.0f, EPSILON);
+
+    Mat2f m3 = Identity<Mat2f>();
+    EXPECT_NEAR(Determinant(m3), 1.0f, EPSILON);
+}
+
+TEST(MatrixTest, Determinant3x3)
+{
+    Mat3f m1(1, 2, 3, 0, 1, 4, 5, 6, 0);
+    EXPECT_NEAR(Determinant(m1), 1.0f, EPSILON);
+
+    Mat3f identity = Identity<Mat3f>();
+    EXPECT_NEAR(Determinant(identity), 1.0f, EPSILON);
+
+    Mat3f singular(1, 2, 3, 2, 4, 6, 3, 6, 9); // Singular matrix
+    EXPECT_NEAR(Determinant(singular), 0.0f, EPSILON);
+}
+
+TEST(MatrixTest, Determinant4x4)
+{
+    Mat4f identity = Identity<Mat4f>();
+    EXPECT_NEAR(Determinant(identity), 1.0f, EPSILON);
+
+    Mat4f m(2, 0, 0, 0, 0, 2, 0, 0, 0, 0, 2, 0, 0, 0, 0, 2);
+    EXPECT_NEAR(Determinant(m), 16.0f, EPSILON); // 2^4 = 16
+}
+
+// ========== IDENTITY MATRIX TESTS ==========
+
+TEST(MatrixTest, IdentityMatrices)
+{
+    Mat2f id2 = Identity<Mat2f>();
+    EXPECT_NEAR(id2(0, 0), 1.0f, EPSILON);
+    EXPECT_NEAR(id2(0, 1), 0.0f, EPSILON);
+    EXPECT_NEAR(id2(1, 0), 0.0f, EPSILON);
+    EXPECT_NEAR(id2(1, 1), 1.0f, EPSILON);
+
+    Mat3f id3 = Identity<Mat3f>();
+    for (size_t i = 0; i < 3; ++i)
+    {
+        for (size_t j = 0; j < 3; ++j)
+        {
+            if (i == j)
+            {
+                EXPECT_NEAR(id3(i, j), 1.0f, EPSILON);
+            }
+            else
+            {
+                EXPECT_NEAR(id3(i, j), 0.0f, EPSILON);
+            }
+        }
+    }
+
+    Mat4f id4 = Identity<Mat4f>();
+    for (size_t i = 0; i < 4; ++i)
+    {
+        for (size_t j = 0; j < 4; ++j)
+        {
+            if (i == j)
+            {
+                EXPECT_NEAR(id4(i, j), 1.0f, EPSILON);
+            }
+            else
+            {
+                EXPECT_NEAR(id4(i, j), 0.0f, EPSILON);
+            }
+        }
+    }
+}
+
+// ========== INVERSE TESTS ==========
+
+TEST(MatrixTest, Inverse3x3)
+{
+    Mat3f m(1, 0, 2, 2, 1, 0, 1, 1, 1);
+    Mat3f inv      = m.Inverse();
+    Mat3f product  = m * inv;
+    Mat3f identity = Identity<Mat3f>();
+
+    // Check if (m * m^-1) ≈ I
+    for (size_t i = 0; i < 3; ++i)
+    {
+        for (size_t j = 0; j < 3; ++j)
+        {
+            EXPECT_NEAR(product(i, j), identity(i, j), EPSILON);
+        }
+    }
+}
+
+TEST(MatrixTest, Inverse4x4)
+{
+    Mat4f m(1, 0, 0, 1, 0, 1, 0, 2, 0, 0, 1, 3, 0, 0, 0, 1);
+    Mat4f inv      = m.Inverse();
+    Mat4f product  = m * inv;
+    Mat4f identity = Identity<Mat4f>();
+
+    // Check if m * m^-1 ≈ I
+    for (size_t i = 0; i < 4; ++i)
+    {
+        for (size_t j = 0; j < 4; ++j)
+        {
+            EXPECT_NEAR(product(i, j), identity(i, j), EPSILON);
+        }
+    }
+}
+
+// ========== EDGE CASES AND ERROR CONDITIONS ==========
+
+TEST(MatrixTest, ZeroDivision)
+{
+    Mat2f m(1, 2, 3, 4);
+    Mat2f result = m / 0.0f;
+
+    // Check for infinity or very large values
+    EXPECT_TRUE(std::isinf(result(0, 0)) || std::abs(result(0, 0)) > 1e30f);
+}
+
+TEST(MatrixTest, IntegerMatrices)
+{
+    IMat2 m(1, 2, 3, 4);
+    EXPECT_EQ(m(0, 0), 1);
+    EXPECT_EQ(m(0, 1), 2);
+    EXPECT_EQ(m(1, 0), 3);
+    EXPECT_EQ(m(1, 1), 4);
+
+    IMat2 m2(5, 6, 7, 8);
+    IMat2 sum = m + m2;
+    EXPECT_EQ(sum(0, 0), 6);
+    EXPECT_EQ(sum(1, 1), 12);
+}
+
+TEST(MatrixTest, LargeValueHandling)
+{
+    Mat2f m(1e20f, 1e20f, 1e20f, 1e20f);
+    Mat2f result = m * 2.0f;
+
+    EXPECT_NEAR(result(0, 0), 2e20f, 1e15f);
+    EXPECT_NEAR(result(1, 1), 2e20f, 1e15f);
+}
+
+TEST(MatrixTest, SmallValueHandling)
+{
+    Mat2f m(1e-20f, 1e-20f, 1e-20f, 1e-20f);
+    Mat2f result = m * 2.0f;
+
+    EXPECT_NEAR(result(0, 0), 2e-20f, 1e-25f);
+    EXPECT_NEAR(result(1, 1), 2e-20f, 1e-25f);
+}
+
+TEST(MatrixTest, ChainedOperations)
+{
+    Mat3f a = Identity<Mat3f>();
+    Mat3f b(2, 0, 0, 0, 2, 0, 0, 0, 2);
+    Mat3f c(1, 1, 1, 1, 1, 1, 1, 1, 1);
+
+    // Test: (a * b) + c - a
+    Mat3f result = ((a * b) + c) - a;
+
+    EXPECT_NEAR(result(0, 0), 2.0f, EPSILON); // (1*2) + 1 - 1 = 2
+    EXPECT_NEAR(result(0, 1), 1.0f, EPSILON); // (0*0) + 1 - 0 = 1
+    EXPECT_NEAR(result(1, 1), 2.0f, EPSILON); // (1*2) + 1 - 1 = 2
+}
+
 TEST(MatrixTest, ConstructionAndIndexingFloat)
 {
     Mat2<float> m2(1.1f, 2.2f, 3.3f, 4.4f);