diff --git a/CMakeLists.txt b/CMakeLists.txt
index c09bfd9..ffc11f4 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -18,10 +18,10 @@ if (NOT DEFINED CMAKE_RUNTIME_OUTPUT_DIRECTORY)
     set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/bin) 
 endif ()
 
-add_executable(libguf_example src/test/example.c src/test/guf_str_impl.c src/test/guf_dict_impl.c)
+add_executable(libguf_example src/test/example.c src/test/guf_str_impl.c src/test/guf_dict_impl.c src/test/guf_linalg_impl.c)
 target_include_directories(libguf_example PRIVATE src src/test)
 
-add_executable(libguf_test src/test/test.cpp src/test/guf_init_impl.c src/test/guf_dbuf_impl.c src/test/guf_str_impl.c src/test/guf_dict_impl.c src/test/guf_rand_impl.c src/test/guf_sort_impl.c)
+add_executable(libguf_test src/test/test.cpp src/test/guf_init_impl.c src/test/guf_dbuf_impl.c src/test/guf_str_impl.c src/test/guf_dict_impl.c src/test/guf_rand_impl.c src/test/guf_sort_impl.c src/test/guf_linalg_impl.c)
 target_include_directories(libguf_test PRIVATE src src/test)
 
 if (NOT DEFINED MSVC)
diff --git a/src/guf_math_linalg.h b/src/guf_linalg.h
similarity index 72%
rename from src/guf_math_linalg.h
rename to src/guf_linalg.h
index e47ee24..41d9952 100644
--- a/src/guf_math_linalg.h
+++ b/src/guf_linalg.h
@@ -1,11 +1,11 @@
-#if defined(GUF_MATH_LINALG_IMPL_STATIC)
-    #define GUF_MATH_LINALG_KWRDS static
+#if defined(GUF_LINALG_IMPL_STATIC)
+    #define GUF_LINALG_KWRDS static
 #else
-    #define GUF_MATH_LINALG_KWRDS
+    #define GUF_LINALG_KWRDS
 #endif
 
-#ifndef GUF_MATH_LINALG_H
-    #define GUF_MATH_LINALG_H
+#ifndef GUF_LINALG_H
+    #define GUF_LINALG_H
     #include "guf_assert.h"
     #include "guf_math.h"
 
@@ -51,44 +51,44 @@
         float angle;
     } guf_axis_angle;
 
-    GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_zero(guf_mat4x4 *m);
-    GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_identity(guf_mat4x4 *m);
-    GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_with_basis(guf_mat4x4 *m, guf_vec3 x_basis, guf_vec3 y_basis, guf_vec3 z_basis);
-    GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_from_quaternion(guf_mat4x4 *rotmat, guf_quaternion q);
+    GUF_LINALG_KWRDS void guf_mat4x4_init_zero(guf_mat4x4 *m);
+    GUF_LINALG_KWRDS void guf_mat4x4_init_identity(guf_mat4x4 *m);
+    GUF_LINALG_KWRDS void guf_mat4x4_init_with_basis(guf_mat4x4 *m, guf_vec3 x_basis, guf_vec3 y_basis, guf_vec3 z_basis);
+    GUF_LINALG_KWRDS void guf_mat4x4_init_from_quaternion(guf_mat4x4 *rotmat, guf_quaternion q);
 
-    GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_zero(guf_mat3x3 *m);
-    GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_identity(guf_mat3x3 *m);
-    GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_with_basis(guf_mat3x3 *m, guf_vec3 x_basis, guf_vec3 y_basis, guf_vec3 z_basis);
-    GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_from_quaternion(guf_mat3x3 *rotmat, guf_quaternion q);
+    GUF_LINALG_KWRDS void guf_mat3x3_init_zero(guf_mat3x3 *m);
+    GUF_LINALG_KWRDS void guf_mat3x3_init_identity(guf_mat3x3 *m);
+    GUF_LINALG_KWRDS void guf_mat3x3_init_with_basis(guf_mat3x3 *m, guf_vec3 x_basis, guf_vec3 y_basis, guf_vec3 z_basis);
+    GUF_LINALG_KWRDS void guf_mat3x3_init_from_quaternion(guf_mat3x3 *rotmat, guf_quaternion q);
 
     #ifdef __cplusplus
         // C++ has no restrict keyword like C99...
-        GUF_MATH_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4 *transposed);
+        GUF_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4 *transposed);
     #else
-        GUF_MATH_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4 * restrict transposed);
+        GUF_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4 * restrict transposed);
     #endif
 
     #ifdef __cplusplus
         // C++ has no restrict keyword like C99...
-        GUF_MATH_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3 *transposed);
+        GUF_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3 *transposed);
     #else
-        GUF_MATH_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3 * restrict transposed);
+        GUF_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3 * restrict transposed);
     #endif
 
     #ifdef __cplusplus
         // C++ has no restrict keyword like C99...
-        GUF_MATH_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4 *a, const guf_mat4x4 *b, guf_mat4x4 *result);
+        GUF_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4 *a, const guf_mat4x4 *b, guf_mat4x4 *result);
     #else
-        GUF_MATH_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4f *a, const guf_mat4x4f *b, guf_mat4x4f * restrict result);
+        GUF_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4 *a, const guf_mat4x4 *b, guf_mat4x4 * restrict result);
     #endif
 
-    GUF_MATH_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *inverted);
+    GUF_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *inverted);
 
     #ifdef __cplusplus
         // C++ has no restrict keyword like C99...
-        GUF_MATH_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3 *b, guf_mat3x3 *result);
+        GUF_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3 *b, guf_mat3x3 *result);
     #else
-        GUF_MATH_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3 *b, guf_mat3x3 * restrict result);
+        GUF_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3 *b, guf_mat3x3 * restrict result);
     #endif
 
     static inline void guf_mat4x4_set_trans(guf_mat4x4 *m, guf_vec3 trans)
@@ -99,12 +99,12 @@
         m->data[2][3] = trans.z;
     }
 
-    GUF_MATH_LINALG_KWRDS void guf_mat4x4_set_rotmat(guf_mat4x4 *m, const guf_mat3x3 *rotmat_3x3);
+    GUF_LINALG_KWRDS void guf_mat4x4_set_rotmat(guf_mat4x4 *m, const guf_mat3x3 *rotmat_3x3);
 
-    GUF_MATH_LINALG_KWRDS bool guf_mat4x4_is_identity(const guf_mat4x4 *m, float eps);
+    GUF_LINALG_KWRDS bool guf_mat4x4_is_identity(const guf_mat4x4 *m, float eps);
 
-    GUF_MATH_LINALG_KWRDS guf_vec3 guf_vec3_transformed(const guf_mat4x4 *mat, const guf_vec3 *v);
-    GUF_MATH_LINALG_KWRDS guf_vec4 guf_vec4_transformed(const guf_mat4x4 *mat, const guf_vec4 *v);
+    GUF_LINALG_KWRDS guf_vec3 guf_vec3_transformed(const guf_mat4x4 *mat, const guf_vec3 *v);
+    GUF_LINALG_KWRDS guf_vec4 guf_vec4_transformed(const guf_mat4x4 *mat, const guf_vec4 *v);
 
     static inline guf_vec2 guf_vec2_add(guf_vec2 a, guf_vec2 b) 
     {
@@ -174,13 +174,13 @@
         return a.x * b.y - a.y * b.x; // The z-component of a 3d cross-product of two vectors in a plane.
     }
 
-    GUF_MATH_LINALG_KWRDS guf_vec2 *guf_vec2_normalise(guf_vec2 *v);    
-    GUF_MATH_LINALG_KWRDS guf_vec3 *guf_vec3_normalise(guf_vec3 *v);
+    GUF_LINALG_KWRDS guf_vec2 *guf_vec2_normalise(guf_vec2 *v);    
+    GUF_LINALG_KWRDS guf_vec3 *guf_vec3_normalise(guf_vec3 *v);
 
-    static inline guf_vec2 guf_vec2_normalised(guf_vec2 v) {return *guf_vec2_normalise(&v);};
-    static inline guf_vec3 guf_vec3_normalised(guf_vec3 v) {return *guf_vec3_normalise(&v);};
+    static inline guf_vec2 guf_vec2_normalised(guf_vec2 v) {return *guf_vec2_normalise(&v);}
+    static inline guf_vec3 guf_vec3_normalised(guf_vec3 v) {return *guf_vec3_normalise(&v);}
 
-    GUF_MATH_LINALG_KWRDS guf_vec3 guf_vec4_to_vec3(guf_vec4 v);
+    GUF_LINALG_KWRDS guf_vec3 guf_vec4_to_vec3(guf_vec4 v);
 
     static inline void guf_quaternion_init_identity(guf_quaternion *q)
     {
@@ -199,7 +199,7 @@
         return q.w * q.w + q.x * q.x + q.y * q.y + q.z * q.z;
     }
 
-    GUF_MATH_LINALG_KWRDS guf_quaternion *guf_quaternion_normalise(guf_quaternion *q);
+    GUF_LINALG_KWRDS guf_quaternion *guf_quaternion_normalise(guf_quaternion *q);
 
     static inline guf_quaternion guf_quaternion_normalised(guf_quaternion q)
     {
@@ -216,28 +216,29 @@
 
     #ifdef __cplusplus
         // C++ has no restrict keyword like C99...
-        GUF_MATH_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf_quaternion *b, guf_quaternion *result);
+        GUF_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf_quaternion *b, guf_quaternion *result);
     #else
-        GUF_MATH_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf_quaternion *b, guf_quaternion * restrict result);
+        GUF_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf_quaternion *b, guf_quaternion * restrict result);
     #endif
 
-    GUF_MATH_LINALG_KWRDS guf_vec3 *guf_vec3_rotate_by_quat(guf_vec3 *p, const guf_quaternion *q);
+    GUF_LINALG_KWRDS guf_vec3 *guf_vec3_rotate_by_quat(guf_vec3 *p, const guf_quaternion *q);
     static inline guf_vec3 guf_vec3_rotated_by_quat(guf_vec3 p, const guf_quaternion *q)
     {
         return *guf_vec3_rotate_by_quat(&p, q);
     }
 
-    GUF_MATH_LINALG_KWRDS guf_quaternion guf_quaternion_from_euler(guf_euler_angles euler);
-    GUF_MATH_LINALG_KWRDS guf_euler_angles guf_quaternion_to_euler(guf_quaternion q);
+    GUF_LINALG_KWRDS guf_quaternion guf_quaternion_from_euler(guf_euler_angles euler);
+    GUF_LINALG_KWRDS guf_euler_angles guf_quaternion_to_euler(guf_quaternion q);
 
 #endif 
 
+#define GUF_LINALG_IMPL_STATIC /* debug */
 
-#if defined(GUF_MATH_LINALG_IMPL) || defined(GUF_MATH_LINALG_IMPL_STATIC)
+#if defined(GUF_LINALG_IMPL) || defined(GUF_LINALG_IMPL_STATIC)
 
 #include <string.h>
 
-GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_zero(guf_mat4x4 *m)
+GUF_LINALG_KWRDS void guf_mat4x4_init_zero(guf_mat4x4 *m)
 {
     GUF_ASSERT(m);
     for (int row = 0; row < 4; ++row) {
@@ -247,7 +248,7 @@ GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_zero(guf_mat4x4 *m)
     }
 }
 
-GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_zero(guf_mat3x3 *m)
+GUF_LINALG_KWRDS void guf_mat3x3_init_zero(guf_mat3x3 *m)
 {
     GUF_ASSERT(m);
     for (int row = 0; row < 3; ++row) {
@@ -257,7 +258,7 @@ GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_zero(guf_mat3x3 *m)
     }
 }
 
-GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_identity(guf_mat4x4 *m)
+GUF_LINALG_KWRDS void guf_mat4x4_init_identity(guf_mat4x4 *m)
 {
     GUF_ASSERT(m);
     guf_mat4x4_init_zero(m);
@@ -266,7 +267,7 @@ GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_identity(guf_mat4x4 *m)
     }
 }
 
-GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_identity(guf_mat3x3 *m)
+GUF_LINALG_KWRDS void guf_mat3x3_init_identity(guf_mat3x3 *m)
 {
     GUF_ASSERT(m);
     guf_mat3x3_init_zero(m);
@@ -275,7 +276,7 @@ GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_identity(guf_mat3x3 *m)
     }
 }
 
-GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_with_basis(guf_mat4x4 *m, guf_vec3 x_basis, guf_vec3 y_basis, guf_vec3 z_basis)
+GUF_LINALG_KWRDS void guf_mat4x4_init_with_basis(guf_mat4x4 *m, guf_vec3 x_basis, guf_vec3 y_basis, guf_vec3 z_basis)
 {
     GUF_ASSERT(m);
 
@@ -298,7 +299,7 @@ GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_with_basis(guf_mat4x4 *m, guf_vec3 x_
     m->data[3][3] = 1;
 }
 
-GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_with_basis(guf_mat3x3 *m, guf_vec3 x_basis, guf_vec3 y_basis, guf_vec3 z_basis)
+GUF_LINALG_KWRDS void guf_mat3x3_init_with_basis(guf_mat3x3 *m, guf_vec3 x_basis, guf_vec3 y_basis, guf_vec3 z_basis)
 {
     GUF_ASSERT(m);
 
@@ -315,10 +316,10 @@ GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_with_basis(guf_mat3x3 *m, guf_vec3 x_
     m->data[2][2] = z_basis.z;
 }
 
-#if defined(__cplusplus) && defined(GUF_MATH_LINALG_IMPL_STATIC)
-GUF_MATH_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4 *a, const guf_mat4x4 *b, guf_mat4x4 *result)
+#if defined(__cplusplus) && defined(GUF_LINALG_IMPL_STATIC)
+GUF_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4 *a, const guf_mat4x4 *b, guf_mat4x4 *result)
 #else
-GUF_MATH_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4 *a, const guf_mat4x4 *b, guf_mat4x4 * restrict result)
+GUF_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4 *a, const guf_mat4x4 *b, guf_mat4x4 * restrict result)
 #endif
 {
     GUF_ASSERT(a && b && result);
@@ -333,10 +334,10 @@ GUF_MATH_LINALG_KWRDS void guf_mat4x4_mul(const guf_mat4x4 *a, const guf_mat4x4
     }
 }
 
-#if defined(__cplusplus) && defined(GUF_MATH_LINALG_IMPL_STATIC)
-GUF_MATH_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3 *b, guf_mat3x3 *result)
+#if defined(__cplusplus) && defined(GUF_LINALG_IMPL_STATIC)
+GUF_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3 *b, guf_mat3x3 *result)
 #else
-GUF_MATH_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3 *b, guf_mat3x3 * restrict result)
+GUF_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3 *b, guf_mat3x3 * restrict result)
 #endif
 {
     GUF_ASSERT(a && b && result);
@@ -351,10 +352,10 @@ GUF_MATH_LINALG_KWRDS void guf_mat3x3_mul(const guf_mat3x3 *a, const guf_mat3x3
     }
 }
 
-#if defined(__cplusplus) && defined(GUF_MATH_LINALG_IMPL_STATIC)
-GUF_MATH_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4 *transposed)
+#if defined(__cplusplus) && defined(GUF_LINALG_IMPL_STATIC)
+GUF_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4 *transposed)
 #else
-GUF_MATH_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4 * restrict transposed)
+GUF_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4 * restrict transposed)
 #endif
 {
     GUF_ASSERT(m && transposed);
@@ -366,10 +367,10 @@ GUF_MATH_LINALG_KWRDS void guf_mat4x4_transposed(const guf_mat4x4 *m, guf_mat4x4
     }
 }
 
-#if defined(__cplusplus) && defined(GUF_MATH_LINALG_IMPL_STATIC)
-GUF_MATH_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3 *transposed)
+#if defined(__cplusplus) && defined(GUF_LINALG_IMPL_STATIC)
+GUF_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3 *transposed)
 #else
-GUF_MATH_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3 * restrict transposed)
+GUF_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3 * restrict transposed)
 #endif
 {
     GUF_ASSERT(m && transposed);
@@ -381,7 +382,7 @@ GUF_MATH_LINALG_KWRDS void guf_mat3x3_transposed(const guf_mat3x3 *m, guf_mat3x3
     }
 }
 
-GUF_MATH_LINALG_KWRDS void guf_mat4x4_set_rotmat(guf_mat4x4 *m, const guf_mat3x3 *rotmat_3x3)
+GUF_LINALG_KWRDS void guf_mat4x4_set_rotmat(guf_mat4x4 *m, const guf_mat3x3 *rotmat_3x3)
 {
     GUF_ASSERT(m && rotmat_3x3)
     for (int row = 0; row < 3; ++row) {
@@ -391,16 +392,16 @@ GUF_MATH_LINALG_KWRDS void guf_mat4x4_set_rotmat(guf_mat4x4 *m, const guf_mat3x3
     }
 }
 
-GUF_MATH_LINALG_KWRDS bool guf_mat4x4_is_identity(const guf_mat4x4 *m, float eps)
+GUF_LINALG_KWRDS bool guf_mat4x4_is_identity(const guf_mat4x4 *m, float eps)
 {
     for (int row = 0; row < 4; ++row) {
         for (int col = 0; col < 4; ++col) {
             if (row == col) { // Check 1 entries.
-                if (fabs(m->data[row][col] - 1) > eps) { 
+                if (!guf_nearly_one_f32(m->data[row][col], eps)) { 
                     return false;
                 }
             } else { // Check 0 entries.
-                if (fabs(m->data[row][col]) > eps) {
+                if (!guf_nearly_one_f32(m->data[row][col], eps)) {
                     return false;
                 }
             }
@@ -431,7 +432,7 @@ static void guf_mat4x4_swap_rows(guf_mat4x4 *m, int row_a_idx, int row_b_idx)
     Note: In the special case of square orthonormal matrices, the matrix inverse can be calculated more efficiently
           and accurately using matrix transposition (e.g. valid 3x3 rotation matrices can be inverted by transposing them).
 */
-GUF_MATH_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *inverted)
+GUF_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *inverted)
 {
     GUF_ASSERT(m && inverted);
     
@@ -441,9 +442,6 @@ GUF_MATH_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *
     const float EPS = 1e-6f;
     const float EPS_ASSERT = 0.05f;
 
-    #define guf_nearly_zero(x, eps) (fabsf((x)) < (eps))
-    #define guf_nearly_one(x, eps)  (fabsf((x) - 1) < (eps))
-
     for (int row_start = 0; row_start < 3; ++row_start) { // 1.) Try to bring into row echelon form.
         // Find pivot, i.e. the element with row in range [row_start, 3] which has
         // the the largest value in the current_column (current_column = row_start)
@@ -451,7 +449,7 @@ GUF_MATH_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *
         float largest = -INFINITY;
         for (int row = row_start; row < 4; ++row) {
             const float candidate = fabsf(rref.data[row][row_start]);
-            if (!guf_nearly_zero(candidate, EPS) && candidate > largest) {
+            if (!guf_nearly_zero_f32(candidate, EPS) && candidate > largest) {
                 largest = candidate;
                 first_nonzero_row = row;
             }   
@@ -469,11 +467,11 @@ GUF_MATH_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *
         const int pivot_row = row_start;
         const int pivot_col = row_start;
         const float pivot_val = rref.data[pivot_row][pivot_col];
-        GUF_ASSERT(!guf_nearly_zero(pivot_val, EPS));
+        GUF_ASSERT(!guf_nearly_zero_f32(pivot_val, EPS));
 
         for (int row = pivot_row + 1; row < 4; ++row) { // Make all values below the pivot element zero. 
             const float val = rref.data[row][pivot_col]; // val: the current value we want to turn into zero.
-            if (guf_nearly_zero(val, EPS)) {
+            if (guf_nearly_zero_f32(val, EPS)) {
                 rref.data[row][pivot_col] = 0;
                 continue;
             }
@@ -492,21 +490,21 @@ GUF_MATH_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *
     }
     
     for (int i = 3; i >= 0; --i) { // 2.) Try to bring into *reduced* row echelon form.
-        if (guf_nearly_zero(rref.data[i][i], EPS)) {
+        if (guf_nearly_zero_f32(rref.data[i][i], EPS)) {
             return false;
         } 
         for (int col = 0; col < 4; ++col) { // Scale row_i to make elem[i][i] 1
             inverted->data[i][col] *= 1.f / rref.data[i][i];
             GUF_ASSERT((col != i) ? rref.data[i][col] == 0 : rref.data[i][col] != 0);
         }
-        GUF_ASSERT(guf_nearly_one(rref.data[i][i], EPS_ASSERT));
+        GUF_ASSERT(guf_nearly_one_f32(rref.data[i][i], EPS_ASSERT));
         rref.data[i][i] = 1; // *= 1.f / rref.data[i][i] (== 1)
 
         for (int row = i - 1; row >= 0; --row) { // Make all elements in column above rref.data[i][i] zero (by adding a scaled row_i to row).
             for (int col = 0; col < 4; ++col) {
                 inverted->data[row][col] += inverted->data[i][col] * -rref.data[row][i];
             }
-            GUF_ASSERT(guf_nearly_zero(rref.data[row][i], EPS_ASSERT));
+            GUF_ASSERT(guf_nearly_zero_f32(rref.data[row][i], EPS_ASSERT));
             rref.data[row][i] = 0; // += rref.data[i][i] * -rref.data[row][i] (== 0)
         }
     }
@@ -514,11 +512,10 @@ GUF_MATH_LINALG_KWRDS bool guf_mat4x4_inverted(const guf_mat4x4 *m, guf_mat4x4 *
     GUF_ASSERT(guf_mat4x4_is_identity(&rref, EPS));
     (void)EPS_ASSERT;
     return true;
-    #undef guf_nearly_zero
 }
 
 // Assumes the last row of mat is [0, 0, 0, 1] and v->w implicitly is 1 (affine transformation)
-GUF_MATH_LINALG_KWRDS guf_vec3 guf_vec3_transformed(const guf_mat4x4 *mat, const guf_vec3 *v) 
+GUF_LINALG_KWRDS guf_vec3 guf_vec3_transformed(const guf_mat4x4 *mat, const guf_vec3 *v) 
 {
     GUF_ASSERT(mat && v);
     guf_vec3 v_transformed = {
@@ -529,7 +526,7 @@ GUF_MATH_LINALG_KWRDS guf_vec3 guf_vec3_transformed(const guf_mat4x4 *mat, const
     return v_transformed;
 }
 
-GUF_MATH_LINALG_KWRDS guf_vec4 guf_vec4_transformed(const guf_mat4x4 *mat, const guf_vec4 *v) 
+GUF_LINALG_KWRDS guf_vec4 guf_vec4_transformed(const guf_mat4x4 *mat, const guf_vec4 *v) 
 {
     GUF_ASSERT(mat && v);
     guf_vec4 v_transformed = {
@@ -541,7 +538,7 @@ GUF_MATH_LINALG_KWRDS guf_vec4 guf_vec4_transformed(const guf_mat4x4 *mat, const
     return v_transformed;
 }
 
-GUF_MATH_LINALG_KWRDS guf_vec2 *guf_vec2_normalise(guf_vec2 *v)
+GUF_LINALG_KWRDS guf_vec2 *guf_vec2_normalise(guf_vec2 *v)
 {
     GUF_ASSERT(v);
     const float mag = guf_vec2_mag(*v);
@@ -554,7 +551,7 @@ GUF_MATH_LINALG_KWRDS guf_vec2 *guf_vec2_normalise(guf_vec2 *v)
     return v;
 }
 
-GUF_MATH_LINALG_KWRDS guf_vec3 *guf_vec3_normalise(guf_vec3 *v)
+GUF_LINALG_KWRDS guf_vec3 *guf_vec3_normalise(guf_vec3 *v)
 { 
     GUF_ASSERT(v);
     const float mag = guf_vec3_mag(*v);
@@ -568,10 +565,10 @@ GUF_MATH_LINALG_KWRDS guf_vec3 *guf_vec3_normalise(guf_vec3 *v)
     return v;
 }
 
-GUF_MATH_LINALG_KWRDS guf_vec3 guf_vec4_to_vec3(guf_vec4 v) 
+GUF_LINALG_KWRDS guf_vec3 guf_vec4_to_vec3(guf_vec4 v) 
 {
     if (guf_isclose_abstol_f32(v.w, 0, 1e-21f)) { // If w is extremely close 0, we treat it as zero.
-        guf_vec3 inf_vec = {(float)INFINITY, (float)INFINITY, float(INFINITY)};
+        guf_vec3 inf_vec = {(float)INFINITY, (float)INFINITY, (float)INFINITY};
         return inf_vec;
     } else if (guf_isclose_reltol_f32(v.w, 1, 1e-9f)) { // If w is almost 1, we treat it as exactly 1. 
         guf_vec3 vec = {.x = v.x, .y = v.y, .z = v.z};
@@ -580,14 +577,14 @@ GUF_MATH_LINALG_KWRDS guf_vec3 guf_vec4_to_vec3(guf_vec4 v)
         guf_vec3 vec = {.x = v.x / v.w, .y = v.y / v.w, .z = v.z / v.w};
         return vec;
     }
-};
+}
 
 /* 
     (Unit)quaternion operations: 
     cf. https://danceswithcode.net/engineeringnotes/quaternions/quaternions.html (last-retrieved 2025-03-06)
 */
 
-GUF_MATH_LINALG_KWRDS guf_quaternion guf_quaternion_from_axis_angle(float angle, guf_vec3 unit_axis)
+GUF_LINALG_KWRDS guf_quaternion guf_quaternion_from_axis_angle(float angle, guf_vec3 unit_axis)
 {
     const float sin_halfangle = sinf(angle / 2.f);
     guf_quaternion q = {
@@ -599,7 +596,7 @@ GUF_MATH_LINALG_KWRDS guf_quaternion guf_quaternion_from_axis_angle(float angle,
     return q;
 }
 
-GUF_MATH_LINALG_KWRDS guf_quaternion guf_quaternion_from_euler(guf_euler_angles euler)
+GUF_LINALG_KWRDS guf_quaternion guf_quaternion_from_euler(guf_euler_angles euler)
 {
     const float r = euler.roll / 2;
     const float p = euler.pitch / 2;
@@ -613,7 +610,7 @@ GUF_MATH_LINALG_KWRDS guf_quaternion guf_quaternion_from_euler(guf_euler_angles
     return q;
 }
 
-GUF_MATH_LINALG_KWRDS guf_euler_angles guf_quaternion_to_euler(guf_quaternion q)
+GUF_LINALG_KWRDS guf_euler_angles guf_quaternion_to_euler(guf_quaternion q)
 {
     guf_euler_angles euler = {0};
     euler.pitch = asinf(2 * (q.w * q.y - q.x * q.z));
@@ -632,7 +629,7 @@ GUF_MATH_LINALG_KWRDS guf_euler_angles guf_quaternion_to_euler(guf_quaternion q)
     return euler;
 }
 
-GUF_MATH_LINALG_KWRDS float guf_quaternion_to_axis_angle(const guf_quaternion *q, guf_vec3 *axis_result)
+GUF_LINALG_KWRDS float guf_quaternion_to_axis_angle(const guf_quaternion *q, guf_vec3 *axis_result)
 {
     GUF_ASSERT(q);
     if (fabsf(q->w - 1) <= GUF_QUAT_TO_AXIS_ANGLE_EPS) { // Singularity: q->w is very close to 1 (identity-quaternions produce no rotation).
@@ -653,7 +650,7 @@ GUF_MATH_LINALG_KWRDS float guf_quaternion_to_axis_angle(const guf_quaternion *q
     return angle;
 }
 
-GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_from_quaternion(guf_mat4x4 *rotmat, guf_quaternion q)
+GUF_LINALG_KWRDS void guf_mat4x4_init_from_quaternion(guf_mat4x4 *rotmat, guf_quaternion q)
 {
     GUF_ASSERT(rotmat);
     
@@ -676,7 +673,7 @@ GUF_MATH_LINALG_KWRDS void guf_mat4x4_init_from_quaternion(guf_mat4x4 *rotmat, g
     rotmat->data[3][3] = 1;
 }
 
-GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_from_quaternion(guf_mat3x3 *rotmat, guf_quaternion q)
+GUF_LINALG_KWRDS void guf_mat3x3_init_from_quaternion(guf_mat3x3 *rotmat, guf_quaternion q)
 {
     GUF_ASSERT(rotmat);
     
@@ -693,7 +690,7 @@ GUF_MATH_LINALG_KWRDS void guf_mat3x3_init_from_quaternion(guf_mat3x3 *rotmat, g
     rotmat->data[2][2] = 1 - 2 * q.x * q.x - 2 * q.y * q.y;
 }
 
-GUF_MATH_LINALG_KWRDS guf_quaternion *guf_quaternion_normalise(guf_quaternion *q)
+GUF_LINALG_KWRDS guf_quaternion *guf_quaternion_normalise(guf_quaternion *q)
 {
     GUF_ASSERT(q);
     const float sq_mag = guf_quaternion_mag_squared(*q);
@@ -712,10 +709,10 @@ GUF_MATH_LINALG_KWRDS guf_quaternion *guf_quaternion_normalise(guf_quaternion *q
     return q;
 }
 
-#if defined(__cplusplus) && defined(GUF_MATH_LINALG_IMPL_STATIC)
-GUF_MATH_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf_quaternion *b, guf_quaternion *result)
+#if defined(__cplusplus) && defined(GUF_LINALG_IMPL_STATIC)
+GUF_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf_quaternion *b, guf_quaternion *result)
 #else
-GUF_MATH_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf_quaternion *b, guf_quaternion * restrict result)
+GUF_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf_quaternion *b, guf_quaternion * restrict result)
 #endif
 {
     GUF_ASSERT(a && b && result);
@@ -726,7 +723,7 @@ GUF_MATH_LINALG_KWRDS void guf_quaternion_mul(const guf_quaternion *a, const guf
     result->z = a->w * b->z - a->x * b->y + a->y * b->x + a->z * b->w; // (r0s3 − r1s2 + r2s1 + r3s0)
 }
 
-GUF_MATH_LINALG_KWRDS guf_vec3 *guf_vec3_rotate_by_quat(guf_vec3 *p, const guf_quaternion *q)
+GUF_LINALG_KWRDS guf_vec3 *guf_vec3_rotate_by_quat(guf_vec3 *p, const guf_quaternion *q)
 {
     // Active rotation of the point/vector p by the quaternion q. 
     GUF_ASSERT(p && q);
@@ -745,9 +742,9 @@ GUF_MATH_LINALG_KWRDS guf_vec3 *guf_vec3_rotate_by_quat(guf_vec3 *p, const guf_q
     return p;
 }
 
-#undef GUF_MATH_LINALG_IMPL
-#undef GUF_MATH_LINALG_IMPL_STATIC
+#undef GUF_LINALG_IMPL
+#undef GUF_LINALG_IMPL_STATIC
 
 #endif /* end impl */
 
-#undef GUF_MATH_LINALG_KWRDS
+#undef GUF_LINALG_KWRDS
diff --git a/src/guf_math.h b/src/guf_math.h
index d80b43a..868cf9c 100644
--- a/src/guf_math.h
+++ b/src/guf_math.h
@@ -72,8 +72,14 @@ static inline uint16_t  guf_absdiff_i16(int16_t a, int16_t b) {return a > b ? (u
 static inline uint32_t  guf_absdiff_i32(int32_t a, int32_t b) {return a > b ? (uint32_t)a - (uint32_t)b : (uint32_t)b - (uint32_t)a;}
 static inline uint64_t  guf_absdiff_i64(int64_t a, int64_t b) {return a > b ? (uint64_t)a - (uint64_t)b : (uint64_t)b - (uint64_t)a;}
 
+static bool guf_nearly_zero_f32(float x, float eps) {return fabsf(x) <= eps;}
+static bool guf_nearly_one_f32(float x, float eps)  {return fabsf(x - 1) <= eps;}
+
+static bool guf_nearly_zero_f64(double x, double eps) {return fabs(x) <= eps;}
+static bool guf_nearly_one_f64(double x, double eps)  {return fabs(x - 1) <= eps;}
+
 /*
-    Floating-point comparison:
+    General floating-point comparison:
 
     cf. https://peps.python.org/pep-0485/ (last-retrieved 2025-03-04)
         https://github.com/PythonCHB/close_pep/blob/master/is_close.py (last-retrieved 2025-03-04)
diff --git a/src/test/guf_linalg_impl.c b/src/test/guf_linalg_impl.c
new file mode 100644
index 0000000..abd171f
--- /dev/null
+++ b/src/test/guf_linalg_impl.c
@@ -0,0 +1,2 @@
+#define GUF_MATH_LINALG_IMPL
+#include "guf_linalg.h"