Vecs and tests

13 files changed, 1849 insertions, 85 deletions

diff --git a/README.md b/README.md
index 728a637..46f7a15 100644
--- a/README.md
+++ b/README.md
@@ -26,15 +26,19 @@ Only core features of C99 were used:
 
 # Ungrateful
 
-Note: we need to manipulate register of SSE so we disable subnormals (Flush to zero).
+Todo: we need to manipulate register of SSE so we disable subnormals (Flush to zero).
 
+You should compile library and program with `UN_DOUBLE_PRECISION` if you want to make `real` type -- 64 bit float (aka don't define it in code).
 Types, Strings, basic memory allocators (platform specific), and so on. Features:
 
+
 + Memory allocation;
 + Strings;
 + easing;
-- splines (qubic);
-- vecs/ivec;
++ vecs;
++ splines (qubic);
+
+Todo:
 - matrix;
 - Quaternions;
 - raycasting;
diff --git a/build.sh b/build.sh
index a2f91ec..82eb864 100755
--- a/build.sh
+++ b/build.sh
@@ -51,6 +51,7 @@ if [[ $1 == "no_tests" ]]; then
 fi
 
 echo 
+echo "---- Ungrateful tests ----"
 
 for test in tests/un/*.c; do
     fname=$(basename -- "$test")
@@ -65,7 +66,8 @@ for test in tests/un/*.c; do
     fi
 done
 
-wait
+echo
+echo "---- Cynic tests ----"
 
 for test in tests/cyn/*.c; do
     fname=$(basename -- "$test")
@@ -82,7 +84,12 @@ done
 
 wait
 
+echo "---- Compiled ----"
 echo
+echo "---- Running tests ----"
+
+passed=0
+failed=0
 
 for case in bin/*; do
     if [[ $1 == "quiet" ]]; then
@@ -92,8 +99,14 @@ for case in bin/*; do
     fi
 
     if [[ $? -eq 0 ]]; then
+        let "passed++"
         echo "[DONE] $case"
     else
+        let "failed++"
         echo "[FAIL] $case"
     fi
 done
+
+echo
+echo "---- Done ----"
+echo "Passed: $passed, Failed: $failed"
diff --git a/src/un_math.c b/src/un_math.c
index 772ef70..42a1c4b 100644
--- a/src/un_math.c
+++ b/src/un_math.c
@@ -6,7 +6,7 @@ real un_m_lerpr(real a, real b, real t) {
 
 real un_m_ease_isiner(real t) {
 #if defined(UN_DOUBLE_PRECISION)
-    return 1.0 - cos((t * PI) / 2.0);
+    return 1.0 - cos((t * PI_D) / 2.0);
 #else
     return 1.0 - cosf((t * PI) / 2.0);
 #endif
@@ -14,7 +14,7 @@ real un_m_ease_isiner(real t) {
 
 real un_m_ease_iosiner(real t) {
 #if defined(UN_DOUBLE_PRECISION)
-    return sin((t * PI) / 2.0f);
+    return sin((t * PI_D) / 2.0f);
 #else
     return sinf((t * PI) / 2.0f);
 #endif
@@ -22,7 +22,7 @@ real un_m_ease_iosiner(real t) {
 
 real un_m_ease_osiner(real t) {
 #if defined(UN_DOUBLE_PRECISION)
-    return -(cos(t * PI) - 1.0) / 2.0;
+    return -(cos(t * PI_D) - 1.0) / 2.0;
 #else
     return -(cosf(t * PI) - 1.0) / 2.0;
 #endif
diff --git a/src/un_splines.c b/src/un_splines.c
index 7b27ee0..c6d2844 100644
--- a/src/un_splines.c
+++ b/src/un_splines.c
@@ -1,40 +1,93 @@
-real un_m_bezierr(real a, real q0, real q1, real b, real t) {
-    real i0, i1, i2, p1, p2;
+void un_m_bezierr(real *v, real a, real q0, real q1, real b, real t) {
+    real i0, i1, i2, p0, p1;
 
     i0 = un_m_lerpr(a,  q0, t);
     i1 = un_m_lerpr(q0, q1, t);
     i2 = un_m_lerpr(q1, b,  t);
 
-    p1 = un_m_lerpr(i0, i1, t);
-    p2 = un_m_lerpr(i1, i2, t);
+    p0 = un_m_lerpr(i0, i1, t);
+    p1 = un_m_lerpr(i1, i2, t);
 
-    return un_m_lerpr(p1, p2, t);
+    *v = un_m_lerpr(p0, p1, t);
 }
 
-f32 un_m_bezierf(f32 a, f32 q0, f32 q1, f32 b, f32 t) {
-    f32 i0, i1, i2, p1, p2;
+void un_m_bezierf(f32 *v, f32 a, f32 q0, f32 q1, f32 b, f32 t) {
+    f32 i0, i1, i2, p0, p1;
 
     i0 = un_m_lerpf(a,  q0, t);
     i1 = un_m_lerpf(q0, q1, t);
     i2 = un_m_lerpf(q1, b,  t);
 
-    p1 = un_m_lerpf(i0, i1, t);
-    p2 = un_m_lerpf(i1, i2, t);
+    p0 = un_m_lerpf(i0, i1, t);
+    p1 = un_m_lerpf(i1, i2, t);
 
-    return un_m_lerpf(p1, p2, t);
+    *v = un_m_lerpf(p0, p1, t);
 }
 
-f64 un_m_bezierd(f64 a, f64 q0, f64 q1, f64 b, f64 t) {
-    f64 i0, i1, i2, p1, p2;
+void un_m_bezierd(f64 *v, f64 a, f64 q0, f64 q1, f64 b, f64 t) {
+    f64 i0, i1, i2, p0, p1;
 
     i0 = un_m_lerpd(a,  q0, t);
     i1 = un_m_lerpd(q0, q1, t);
     i2 = un_m_lerpd(q1, b,  t);
 
-    p1 = un_m_lerpd(i0, i1, t);
-    p2 = un_m_lerpd(i1, i2, t);
+    p0 = un_m_lerpd(i0, i1, t);
+    p1 = un_m_lerpd(i1, i2, t);
 
-    return un_m_lerpd(p1, p2, t);
+    *v = un_m_lerpd(p0, p1, t);
 }
 
-// add v2-4 variants
+void un_m_bezier2r(real *v, real *a, real *q0, real *q1, real *b, real t) {
+    un_m_bezierr(v + 0, a[0], q0[0], q1[0], b[0], t);
+    un_m_bezierr(v + 1, a[1], q0[1], q1[1], b[1], t);
+}
+
+void un_m_bezier3r(real *v, real *a, real *q0, real *q1, real *b, real t) {
+    un_m_bezierr(v + 0, a[0], q0[0], q1[0], b[0], t);
+    un_m_bezierr(v + 1, a[1], q0[1], q1[1], b[1], t);
+    un_m_bezierr(v + 2, a[2], q0[2], q1[2], b[2], t);
+}
+
+void un_m_bezier4r(real *v, real *a, real *q0, real *q1, real *b, real t) {
+    un_m_bezierr(v + 0, a[0], q0[0], q1[0], b[0], t);
+    un_m_bezierr(v + 1, a[1], q0[1], q1[1], b[1], t);
+    un_m_bezierr(v + 2, a[2], q0[2], q1[2], b[2], t);
+    un_m_bezierr(v + 3, a[3], q0[3], q1[3], b[3], t);
+}
+
+void un_m_bezier2f(f32 *v, f32 *a, f32 *q0, f32 *q1, f32 *b, f32 t) {
+    un_m_bezierf(v + 0, a[0], q0[0], q1[0], b[0], t);
+    un_m_bezierf(v + 1, a[1], q0[1], q1[1], b[1], t);
+}
+
+void un_m_bezier3f(f32 *v, f32 *a, f32 *q0, f32 *q1, f32 *b, f32 t) {
+    un_m_bezierf(v + 0, a[0], q0[0], q1[0], b[0], t);
+    un_m_bezierf(v + 1, a[1], q0[1], q1[1], b[1], t);
+    un_m_bezierf(v + 2, a[2], q0[2], q1[2], b[2], t);
+}
+
+void un_m_bezier4f(f32 *v, f32 *a, f32 *q0, f32 *q1, f32 *b, f32 t) {
+    un_m_bezierf(v + 0, a[0], q0[0], q1[0], b[0], t);
+    un_m_bezierf(v + 1, a[1], q0[1], q1[1], b[1], t);
+    un_m_bezierf(v + 2, a[2], q0[2], q1[2], b[2], t);
+    un_m_bezierf(v + 3, a[3], q0[3], q1[3], b[3], t);
+}
+
+
+void un_m_bezier2d(f64 *v, f64 *a, f64 *q0, f64 *q1, f64 *b, f64 t) {
+    un_m_bezierd(v + 0, a[0], q0[0], q1[0], b[0], t);
+    un_m_bezierd(v + 1, a[1], q0[1], q1[1], b[1], t);
+}
+
+void un_m_bezier3d(f64 *v, f64 *a, f64 *q0, f64 *q1, f64 *b, f64 t) {
+    un_m_bezierd(v + 0, a[0], q0[0], q1[0], b[0], t);
+    un_m_bezierd(v + 1, a[1], q0[1], q1[1], b[1], t);
+    un_m_bezierd(v + 2, a[2], q0[2], q1[2], b[2], t);
+}
+
+void un_m_bezier4d(f64 *v, f64 *a, f64 *q0, f64 *q1, f64 *b, f64 t) {
+    un_m_bezierd(v + 0, a[0], q0[0], q1[0], b[0], t);
+    un_m_bezierd(v + 1, a[1], q0[1], q1[1], b[1], t);
+    un_m_bezierd(v + 2, a[2], q0[2], q1[2], b[2], t);
+    un_m_bezierd(v + 3, a[3], q0[3], q1[3], b[3], t);
+}
diff --git a/src/un_vec.c b/src/un_vec.c
index 1e1ccda..d6c984e 100644
--- a/src/un_vec.c
+++ b/src/un_vec.c
@@ -1,31 +1,3 @@
-
-
-
-/*
- *
- * vector from two vectors
- *
- * add
- * sub
- *
- * add scalar
- * sub scalar
- *
- * mul scalar
- * div scalar
- *
- * distance
- * discance sqr
- *
- * magnitude
- * magnitude sqr
- *
- * dot, cross, surface_cross, normalize, reflect, flatten, project
- *
- * 
- *
- */
-
 void un_m_add2f(f32 *v, f32 *a, f32 *b) {
     v[0] = a[0] + b[0];
     v[1] = a[1] + b[1];
@@ -199,6 +171,23 @@ void un_m_magnitude_sqr4f(f32 *v, f32 *a) {
     *v = a[0] * a[0] + a[1] * a[1] + a[2] * a[2] + a[3] * a[3];
 }
 
+void un_m_distance2f(f32 *v, f32 *a, f32 *b) {
+    f32 f[2];
+    un_m_sub2f(f, b, a);
+    *v = sqrtf(f[0] * f[0] + f[1] * f[1]);
+}
+
+void un_m_distance3f(f32 *v, f32 *a, f32 *b) {
+    f32 f[3];
+    un_m_sub3f(f, b, a);
+    *v = sqrtf(f[0] * f[0] + f[1] * f[1] + f[2] * f[2]);
+}
+
+void un_m_distance4f(f32 *v, f32 *a, f32 *b) {
+    f32 f[4];
+    un_m_sub4f(f, b, a);
+    *v = sqrtf(f[0] * f[0] + f[1] * f[1] + f[2] * f[2] + f[3] * f[3]);
+}
 
 void un_m_distance_sqr2f(f32 *v, f32 *a, f32 *b) {
     f32 f[2];
@@ -221,25 +210,25 @@ void un_m_distance_sqr4f(f32 *v, f32 *a, f32 *b) {
 void un_m_normalize2f(f32 *v, f32 *a) {
     f32 f;
     un_m_magnitude2f(&f, a);
-    v[0] /= f;
-    v[1] /= f;
+    v[0] = a[0] / f;
+    v[1] = a[1] / f;
 }
 
 void un_m_normalize3f(f32 *v, f32 *a) {
     f32 f;
     un_m_magnitude3f(&f, a);
-    v[0] /= f;
-    v[1] /= f;
-    v[2] /= f;
+    v[0] = a[0] / f;
+    v[1] = a[1] / f;
+    v[2] = a[2] / f;
 }
 
 void un_m_normalize4f(f32 *v, f32 *a) {
     f32 f;
     un_m_magnitude4f(&f, a);
-    v[0] /= f;
-    v[1] /= f;
-    v[2] /= f;
-    v[3] /= f;
+    v[0] = a[0] / f;
+    v[1] = a[1] / f;
+    v[2] = a[2] / f;
+    v[3] = a[3] / f;
 }
 
 void un_m_project2f(f32 *v, f32 *a, f32 *onto) {
@@ -292,3 +281,27 @@ void un_m_project4f(f32 *v, f32 *a, f32 *onto) {
     v[2] = onto[2] * scale;
     v[3] = onto[3] * scale;
 }
+
+void un_m_reflect2f(f32 *v, f32 *a, f32 *normal) {
+    f32 dot;
+
+    un_m_dot2f(&dot, a, normal);
+    un_m_mul_scalar2f(v, normal, 2.0f * dot);
+    un_m_sub2f(v, v, a);
+}
+
+void un_m_reflect3f(f32 *v, f32 *a, f32 *normal) {
+    f32 dot;
+
+    un_m_dot3f(&dot, a, normal);
+    un_m_mul_scalar3f(v, normal, 2.0f * dot);
+    un_m_sub3f(v, v, a);
+}
+
+void un_m_reflect4f(f32 *v, f32 *a, f32 *normal) {
+    f32 dot;
+
+    un_m_dot4f(&dot, a, normal);
+    un_m_mul_scalar4f(v, normal, 2.0f * dot);
+    un_m_sub4f(v, v, a);
+}
diff --git a/src/un_vecd.c b/src/un_vecd.c
new file mode 100644
index 0000000..0d138a4
--- /dev/null
+++ b/src/un_vecd.c
@@ -0,0 +1,307 @@
+void un_m_add2d(f64 *v, f64 *a, f64 *b) {
+    v[0] = a[0] + b[0];
+    v[1] = a[1] + b[1];
+}
+
+void un_m_add3d(f64 *v, f64 *a, f64 *b) {
+    v[0] = a[0] + b[0];
+    v[1] = a[1] + b[1];
+    v[2] = a[2] + b[2];
+}
+
+void un_m_add4d(f64 *v, f64 *a, f64 *b) {
+    v[0] = a[0] + b[0];
+    v[1] = a[1] + b[1];
+    v[2] = a[2] + b[2];
+    v[3] = a[3] + b[3];
+}
+
+void un_m_sub2d(f64 *v, f64 *a, f64 *b) {
+    v[0] = a[0] - b[0];
+    v[1] = a[1] - b[1];
+}
+
+void un_m_sub3d(f64 *v, f64 *a, f64 *b) {
+    v[0] = a[0] - b[0];
+    v[1] = a[1] - b[1];
+    v[2] = a[2] - b[2];
+}
+
+void un_m_sub4d(f64 *v, f64 *a, f64 *b) {
+    v[0] = a[0] - b[0];
+    v[1] = a[1] - b[1];
+    v[2] = a[2] - b[2];
+    v[3] = a[3] - b[3];
+}
+
+void un_m_add_scalar2d(f64 *v, f64 *a, f64 scalar) {
+    v[0] = a[0] + scalar;
+    v[1] = a[1] + scalar;
+}
+
+void un_m_add_scalar3d(f64 *v, f64 *a, f64 scalar) {
+    v[0] = a[0] + scalar;
+    v[1] = a[1] + scalar;
+    v[2] = a[2] + scalar;
+}
+
+void un_m_add_scalar4d(f64 *v, f64 *a, f64 scalar) {
+    v[0] = a[0] + scalar;
+    v[1] = a[1] + scalar;
+    v[2] = a[2] + scalar;
+    v[3] = a[3] + scalar;
+}
+
+void un_m_sub_scalar2d(f64 *v, f64 *a, f64 scalar) {
+    v[0] = a[0] - scalar;
+    v[1] = a[1] - scalar;
+}
+
+void un_m_sub_scalar3d(f64 *v, f64 *a, f64 scalar) {
+    v[0] = a[0] - scalar;
+    v[1] = a[1] - scalar;
+    v[2] = a[2] - scalar;
+}
+
+void un_m_sub_scalar4d(f64 *v, f64 *a, f64 scalar) {
+    v[0] = a[0] - scalar;
+    v[1] = a[1] - scalar;
+    v[2] = a[2] - scalar;
+    v[3] = a[3] - scalar;
+}
+
+void un_m_mul_scalar2d(f64 *v, f64 *a, f64 scalar) {
+    v[0] = a[0] * scalar;
+    v[1] = a[1] * scalar;
+}
+
+void un_m_mul_scalar3d(f64 *v, f64 *a, f64 scalar) {
+    v[0] = a[0] * scalar;
+    v[1] = a[1] * scalar;
+    v[2] = a[2] * scalar;
+}
+
+void un_m_mul_scalar4d(f64 *v, f64 *a, f64 scalar) {
+    v[0] = a[0] * scalar;
+    v[1] = a[1] * scalar;
+    v[2] = a[2] * scalar;
+    v[3] = a[3] * scalar;
+}
+
+void un_m_div_scalar2d(f64 *v, f64 *a, f64 scalar) {
+    v[0] = a[0] / scalar;
+    v[1] = a[1] / scalar;
+}
+
+void un_m_div_scalar3d(f64 *v, f64 *a, f64 scalar) {
+    v[0] = a[0] / scalar;
+    v[1] = a[1] / scalar;
+    v[2] = a[2] / scalar;
+}
+
+void un_m_div_scalar4d(f64 *v, f64 *a, f64 scalar) {
+    v[0] = a[0] / scalar;
+    v[1] = a[1] / scalar;
+    v[2] = a[2] / scalar;
+    v[3] = a[3] / scalar;
+}
+
+void un_m_dot2d(f64 *v, f64 *a, f64 *b) {
+    *v = a[0] * b[0] + a[1] * b[1];
+}
+
+void un_m_dot3d(f64 *v, f64 *a, f64 *b) {
+    *v = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
+}
+
+void un_m_dot4d(f64 *v, f64 *a, f64 *b) {
+    *v = a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[3];
+}
+
+void un_m_hadamard2d(f64 *v, f64 *a, f64 *b) {
+    v[0] = a[0] * b[0];
+    v[1] = a[1] * b[1];
+}
+
+void un_m_hadamard3d(f64 *v, f64 *a, f64 *b) {
+    v[0] = a[0] * b[0];
+    v[1] = a[1] * b[1];
+    v[2] = a[2] * b[2];
+}
+
+void un_m_hadamard4d(f64 *v, f64 *a, f64 *b) {
+    v[0] = a[0] * b[0];
+    v[1] = a[1] * b[1];
+    v[2] = a[2] * b[2];
+    v[3] = a[3] * b[3];
+}
+
+void un_m_cross2d(f64 *v, f64 *a) {
+    v[0] = a[1];
+    v[1] =-a[0];
+}
+
+void un_m_cross3d(f64 *v, f64 *a, f64 *b) {
+	v[0] = a[1] * b[2] - a[2] * b[1];
+	v[1] = a[2] * b[0] - a[0] * b[2];
+	v[2] = a[0] * b[1] - a[1] * b[0];
+}
+
+void un_m_magnitude2d(f64 *v, f64 *a) {
+    *v = sqrt(a[0] * a[0] + a[1] * a[1]);
+}
+
+void un_m_magnitude3d(f64 *v, f64 *a) {
+    *v = sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);
+}
+
+void un_m_magnitude4d(f64 *v, f64 *a) {
+    *v = sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2] + a[3] * a[3]);
+}
+
+void un_m_magnitude_sqr2d(f64 *v, f64 *a) {
+    *v = a[0] * a[0] + a[1] * a[1];
+}
+
+void un_m_magnitude_sqr3d(f64 *v, f64 *a) {
+    *v = a[0] * a[0] + a[1] * a[1] + a[2] * a[2];
+}
+
+void un_m_magnitude_sqr4d(f64 *v, f64 *a) {
+    *v = a[0] * a[0] + a[1] * a[1] + a[2] * a[2] + a[3] * a[3];
+}
+
+void un_m_distance2d(f64 *v, f64 *a, f64 *b) {
+    f64 f[2];
+    un_m_sub2d(f, b, a);
+    *v = sqrt(f[0] * f[0] + f[1] * f[1]);
+}
+
+void un_m_distance3d(f64 *v, f64 *a, f64 *b) {
+    f64 f[3];
+    un_m_sub3d(f, b, a);
+    *v = sqrt(f[0] * f[0] + f[1] * f[1] + f[2] * f[2]);
+}
+
+void un_m_distance4d(f64 *v, f64 *a, f64 *b) {
+    f64 f[4];
+    un_m_sub4d(f, b, a);
+    *v = sqrt(f[0] * f[0] + f[1] * f[1] + f[2] * f[2] + f[3] * f[3]);
+}
+
+void un_m_distance_sqr2d(f64 *v, f64 *a, f64 *b) {
+    f64 f[2];
+    un_m_sub2d(f, b, a);
+    *v = f[0] * f[0] + f[1] * f[1];
+}
+
+void un_m_distance_sqr3d(f64 *v, f64 *a, f64 *b) {
+    f64 f[3];
+    un_m_sub3d(f, b, a);
+    *v = f[0] * f[0] + f[1] * f[1] + f[2] * f[2];
+}
+
+void un_m_distance_sqr4d(f64 *v, f64 *a, f64 *b) {
+    f64 f[4];
+    un_m_sub4d(f, b, a);
+    *v = f[0] * f[0] + f[1] * f[1] + f[2] * f[2] + f[3] * f[3];
+}
+
+void un_m_normalize2d(f64 *v, f64 *a) {
+    f64 f;
+    un_m_magnitude2d(&f, a);
+    v[0] = a[0] / f;
+    v[1] = a[1] / f;
+}
+
+void un_m_normalize3d(f64 *v, f64 *a) {
+    f64 f;
+    un_m_magnitude3d(&f, a);
+    v[0] = a[0] / f;
+    v[1] = a[1] / f;
+    v[2] = a[2] / f;
+}
+
+void un_m_normalize4d(f64 *v, f64 *a) {
+    f64 f;
+    un_m_magnitude4d(&f, a);
+    v[0] = a[0] / f;
+    v[1] = a[1] / f;
+    v[2] = a[2] / f;
+    v[3] = a[3] / f;
+}
+
+void un_m_project2d(f64 *v, f64 *a, f64 *onto) {
+    f64 dot, magn, scale;
+
+    un_m_dot2d(&dot, a, onto);
+    un_m_magnitude_sqr2d(&magn, onto);
+
+    if (magn < EPSILON) {
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+        return;
+    }
+
+    scale = dot / magn;
+    v[0] = onto[0] * scale;
+    v[1] = onto[1] * scale;
+}
+
+void un_m_project3d(f64 *v, f64 *a, f64 *onto) {
+    f64 dot, magn, scale;
+
+    un_m_dot3d(&dot, a, onto);
+    un_m_magnitude_sqr3d(&magn, onto);
+
+    if (magn < EPSILON_D) {
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+        return;
+    }
+
+    scale = dot / magn;
+    v[0] = onto[0] * scale;
+    v[1] = onto[1] * scale;
+    v[2] = onto[2] * scale;
+}
+
+void un_m_project4d(f64 *v, f64 *a, f64 *onto) {
+    f64 dot, magn, scale;
+
+    un_m_dot4d(&dot, a, onto);
+    un_m_magnitude_sqr4d(&magn, onto);
+
+    if (magn < EPSILON_D) {
+        un_memory_set((void*)v, 0, sizeof(*v) * 4.0);
+        return;
+    }
+
+    scale = dot / magn;
+    v[0] = onto[0] * scale;
+    v[1] = onto[1] * scale;
+    v[2] = onto[2] * scale;
+    v[3] = onto[3] * scale;
+}
+
+void un_m_reflect2d(f64 *v, f64 *a, f64 *normal) {
+    f64 dot;
+
+    un_m_dot2d(&dot, a, normal);
+    un_m_mul_scalar2d(v, normal, 2.0 * dot);
+    un_m_sub2d(v, v, a);
+}
+
+void un_m_reflect3d(f64 *v, f64 *a, f64 *normal) {
+    f64 dot;
+
+    un_m_dot3d(&dot, a, normal);
+    un_m_mul_scalar3d(v, normal, 2.0 * dot);
+    un_m_sub3d(v, v, a);
+}
+
+void un_m_reflect4d(f64 *v, f64 *a, f64 *normal) {
+    f64 dot;
+
+    un_m_dot4d(&dot, a, normal);
+    un_m_mul_scalar4d(v, normal, 2.0 * dot);
+    un_m_sub4d(v, v, a);
+}
diff --git a/src/un_vecr.c b/src/un_vecr.c
new file mode 100644
index 0000000..0fcf672
--- /dev/null
+++ b/src/un_vecr.c
@@ -0,0 +1,353 @@
+void un_m_add2r(real *v, real *a, real *b) {
+    v[0] = a[0] + b[0];
+    v[1] = a[1] + b[1];
+}
+
+void un_m_add3r(real *v, real *a, real *b) {
+    v[0] = a[0] + b[0];
+    v[1] = a[1] + b[1];
+    v[2] = a[2] + b[2];
+}
+
+void un_m_add4r(real *v, real *a, real *b) {
+    v[0] = a[0] + b[0];
+    v[1] = a[1] + b[1];
+    v[2] = a[2] + b[2];
+    v[3] = a[3] + b[3];
+}
+
+void un_m_sub2r(real *v, real *a, real *b) {
+    v[0] = a[0] - b[0];
+    v[1] = a[1] - b[1];
+}
+
+void un_m_sub3r(real *v, real *a, real *b) {
+    v[0] = a[0] - b[0];
+    v[1] = a[1] - b[1];
+    v[2] = a[2] - b[2];
+}
+
+void un_m_sub4r(real *v, real *a, real *b) {
+    v[0] = a[0] - b[0];
+    v[1] = a[1] - b[1];
+    v[2] = a[2] - b[2];
+    v[3] = a[3] - b[3];
+}
+
+void un_m_add_scalar2r(real *v, real *a, real scalar) {
+    v[0] = a[0] + scalar;
+    v[1] = a[1] + scalar;
+}
+
+void un_m_add_scalar3r(real *v, real *a, real scalar) {
+    v[0] = a[0] + scalar;
+    v[1] = a[1] + scalar;
+    v[2] = a[2] + scalar;
+}
+
+void un_m_add_scalar4r(real *v, real *a, real scalar) {
+    v[0] = a[0] + scalar;
+    v[1] = a[1] + scalar;
+    v[2] = a[2] + scalar;
+    v[3] = a[3] + scalar;
+}
+
+void un_m_sub_scalar2r(real *v, real *a, real scalar) {
+    v[0] = a[0] - scalar;
+    v[1] = a[1] - scalar;
+}
+
+void un_m_sub_scalar3r(real *v, real *a, real scalar) {
+    v[0] = a[0] - scalar;
+    v[1] = a[1] - scalar;
+    v[2] = a[2] - scalar;
+}
+
+void un_m_sub_scalar4r(real *v, real *a, real scalar) {
+    v[0] = a[0] - scalar;
+    v[1] = a[1] - scalar;
+    v[2] = a[2] - scalar;
+    v[3] = a[3] - scalar;
+}
+
+void un_m_mul_scalar2r(real *v, real *a, real scalar) {
+    v[0] = a[0] * scalar;
+    v[1] = a[1] * scalar;
+}
+
+void un_m_mul_scalar3r(real *v, real *a, real scalar) {
+    v[0] = a[0] * scalar;
+    v[1] = a[1] * scalar;
+    v[2] = a[2] * scalar;
+}
+
+void un_m_mul_scalar4r(real *v, real *a, real scalar) {
+    v[0] = a[0] * scalar;
+    v[1] = a[1] * scalar;
+    v[2] = a[2] * scalar;
+    v[3] = a[3] * scalar;
+}
+
+void un_m_div_scalar2r(real *v, real *a, real scalar) {
+    v[0] = a[0] / scalar;
+    v[1] = a[1] / scalar;
+}
+
+void un_m_div_scalar3r(real *v, real *a, real scalar) {
+    v[0] = a[0] / scalar;
+    v[1] = a[1] / scalar;
+    v[2] = a[2] / scalar;
+}
+
+void un_m_div_scalar4r(real *v, real *a, real scalar) {
+    v[0] = a[0] / scalar;
+    v[1] = a[1] / scalar;
+    v[2] = a[2] / scalar;
+    v[3] = a[3] / scalar;
+}
+
+void un_m_dot2r(real *v, real *a, real *b) {
+    *v = a[0] * b[0] + a[1] * b[1];
+}
+
+void un_m_dot3r(real *v, real *a, real *b) {
+    *v = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
+}
+
+void un_m_dot4r(real *v, real *a, real *b) {
+    *v = a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[3];
+}
+
+void un_m_hadamard2r(real *v, real *a, real *b) {
+    v[0] = a[0] * b[0];
+    v[1] = a[1] * b[1];
+}
+
+void un_m_hadamard3r(real *v, real *a, real *b) {
+    v[0] = a[0] * b[0];
+    v[1] = a[1] * b[1];
+    v[2] = a[2] * b[2];
+}
+
+void un_m_hadamard4r(real *v, real *a, real *b) {
+    v[0] = a[0] * b[0];
+    v[1] = a[1] * b[1];
+    v[2] = a[2] * b[2];
+    v[3] = a[3] * b[3];
+}
+
+void un_m_cross2r(real *v, real *a) {
+    v[0] = a[1];
+    v[1] =-a[0];
+}
+
+void un_m_cross3r(real *v, real *a, real *b) {
+	v[0] = a[1] * b[2] - a[2] * b[1];
+	v[1] = a[2] * b[0] - a[0] * b[2];
+	v[2] = a[0] * b[1] - a[1] * b[0];
+}
+
+void un_m_magnitude2r(real *v, real *a) {
+#if defined(UN_DOUBLE_PRECISION)
+    *v = sqrt(a[0] * a[0] + a[1] * a[1]);
+#else
+    *v = sqrtf(a[0] * a[0] + a[1] * a[1]);
+#endif
+}
+
+void un_m_magnitude3r(real *v, real *a) {
+#if defined(UN_DOUBLE_PRECISION)
+    *v = sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);
+#else
+    *v = sqrtf(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);
+#endif
+}
+
+void un_m_magnitude4r(real *v, real *a) {
+#if defined(UN_DOUBLE_PRECISION)
+    *v = sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2] + a[3] * a[3]);
+#else
+    *v = sqrtf(a[0] * a[0] + a[1] * a[1] + a[2] * a[2] + a[3] * a[3]);
+#endif
+}
+
+void un_m_magnitude_sqr2r(real *v, real *a) {
+    *v = a[0] * a[0] + a[1] * a[1];
+}
+
+void un_m_magnitude_sqr3r(real *v, real *a) {
+    *v = a[0] * a[0] + a[1] * a[1] + a[2] * a[2];
+}
+
+void un_m_magnitude_sqr4r(real *v, real *a) {
+    *v = a[0] * a[0] + a[1] * a[1] + a[2] * a[2] + a[3] * a[3];
+}
+
+void un_m_distance2r(real *v, real *a, real *b) {
+    real f[2];
+    un_m_sub2r(f, b, a);
+#if defined(UN_DOUBLE_PRECISION)
+    *v = sqrt(f[0] * f[0] + f[1] * f[1]);
+#else
+    *v = sqrtf(f[0] * f[0] + f[1] * f[1]);
+#endif
+}
+
+void un_m_distance3r(real *v, real *a, real *b) {
+    real f[3];
+    un_m_sub3r(f, b, a);
+#if defined(UN_DOUBLE_PRECISION)
+    *v = sqrt(f[0] * f[0] + f[1] * f[1] + f[2] * f[2]);
+#else
+    *v = sqrtf(f[0] * f[0] + f[1] * f[1] + f[2] * f[2]);
+#endif
+}
+
+void un_m_distance4r(real *v, real *a, real *b) {
+    real f[4];
+    un_m_sub4r(f, b, a);
+#if defined(UN_DOUBLE_PRECISION)
+    *v = sqrt(f[0] * f[0] + f[1] * f[1] + f[2] * f[2] + f[3] * f[3]);
+#else
+    *v = sqrtf(f[0] * f[0] + f[1] * f[1] + f[2] * f[2] + f[3] * f[3]);
+#endif
+}
+
+void un_m_distance_sqr2r(real *v, real *a, real *b) {
+    real f[2];
+    un_m_sub2r(f, b, a);
+    *v = f[0] * f[0] + f[1] * f[1];
+}
+
+void un_m_distance_sqr3r(real *v, real *a, real *b) {
+    real f[3];
+    un_m_sub3r(f, b, a);
+    *v = f[0] * f[0] + f[1] * f[1] + f[2] * f[2];
+}
+
+void un_m_distance_sqr4r(real *v, real *a, real *b) {
+    real f[4];
+    un_m_sub4r(f, b, a);
+    *v = f[0] * f[0] + f[1] * f[1] + f[2] * f[2] + f[3] * f[3];
+}
+
+void un_m_normalize2r(real *v, real *a) {
+    real f;
+    un_m_magnitude2r(&f, a);
+    v[0] = a[0] / f;
+    v[1] = a[1] / f;
+}
+
+void un_m_normalize3r(real *v, real *a) {
+    real f;
+    un_m_magnitude3r(&f, a);
+    v[0] = a[0] / f;
+    v[1] = a[1] / f;
+    v[2] = a[2] / f;
+}
+
+void un_m_normalize4r(real *v, real *a) {
+    real f;
+    un_m_magnitude4r(&f, a);
+    v[0] = a[0] / f;
+    v[1] = a[1] / f;
+    v[2] = a[2] / f;
+    v[3] = a[3] / f;
+}
+
+void un_m_project2r(real *v, real *a, real *onto) {
+    real dot, magn, scale;
+
+    un_m_dot2r(&dot, a, onto);
+    un_m_magnitude_sqr2r(&magn, onto);
+
+
+#if defined(UN_DOUBLE_PRECISION)
+    if (magn < EPSILON_D) {
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+        return;
+    }
+#else
+    if (magn < EPSILON) {
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+        return;
+    }
+#endif
+
+    scale = dot / magn;
+    v[0] = onto[0] * scale;
+    v[1] = onto[1] * scale;
+}
+
+void un_m_project3r(real *v, real *a, real *onto) {
+    real dot, magn, scale;
+
+    un_m_dot3r(&dot, a, onto);
+    un_m_magnitude_sqr3r(&magn, onto);
+
+#if defined(UN_DOUBLE_PRECISION)
+    if (magn < EPSILON_D) {
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+        return;
+    }
+#else
+    if (magn < EPSILON) {
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+        return;
+    }
+#endif
+
+    scale = dot / magn;
+    v[0] = onto[0] * scale;
+    v[1] = onto[1] * scale;
+    v[2] = onto[2] * scale;
+}
+
+void un_m_project4r(real *v, real *a, real *onto) {
+    real dot, magn, scale;
+
+    un_m_dot4r(&dot, a, onto);
+    un_m_magnitude_sqr4r(&magn, onto);
+
+#if defined(UN_DOUBLE_PRECISION)
+    if (magn < EPSILON_D) {
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+        return;
+    }
+#else
+    if (magn < EPSILON) {
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+        return;
+    }
+#endif
+
+    scale = dot / magn;
+    v[0] = onto[0] * scale;
+    v[1] = onto[1] * scale;
+    v[2] = onto[2] * scale;
+    v[3] = onto[3] * scale;
+}
+
+void un_m_reflect2r(real *v, real *a, real *normal) {
+    real dot;
+
+    un_m_dot2r(&dot, a, normal);
+    un_m_mul_scalar2r(v, normal, 2.0 * dot);
+    un_m_sub2r(v, v, a);
+}
+
+void un_m_reflect3r(real *v, real *a, real *normal) {
+    real dot;
+
+    un_m_dot3r(&dot, a, normal);
+    un_m_mul_scalar3r(v, normal, 2.0 * dot);
+    un_m_sub3r(v, v, a);
+}
+
+void un_m_reflect4r(real *v, real *a, real *normal) {
+    real dot;
+
+    un_m_dot4r(&dot, a, normal);
+    un_m_mul_scalar4r(v, normal, 2.0 * dot);
+    un_m_sub4r(v, v, a);
+}
diff --git a/src/ungrateful.c b/src/ungrateful.c
index 6bfd381..61d8185 100644
--- a/src/ungrateful.c
+++ b/src/ungrateful.c
@@ -3,6 +3,11 @@
 #include "un_memory.c"
 #include "un_strings.c"
 #include "un_list.c"
+
 #include "un_math.c"
+
 #include "un_vec.c"
+#include "un_vecd.c"
+#include "un_vecr.c"
+
 #include "un_splines.c"
diff --git a/src/ungrateful.h b/src/ungrateful.h
index ff8c24b..bbbea47 100644
--- a/src/ungrateful.h
+++ b/src/ungrateful.h
@@ -97,45 +97,47 @@ typedef int8_t b8;
 typedef float  f32;
 typedef double f64;
 
-#if defined(UN_DOUBLE_PRECISION)
-typedef double real;
-
 #ifndef PI
-#   define PI 3.14159265358979323846
+#   define PI 3.14159265358979323846f
 #endif
 
 #ifndef EPSILON
-#   define EPSILON 0.000001
+#   define EPSILON 0.000001f
 #endif
 
 #ifndef DEG2RAD
-#   define DEG2RAD (PI/180.0)
+#   define DEG2RAD (PI/180.0f)
 #endif
 
 #ifndef RAD2DEG
-#   define RAD2DEG (180.0/PI)
+#   define RAD2DEG (180.0f/PI)
 #endif
 
+#if defined(UN_DOUBLE_PRECISION)
+    typedef double real;
 #else
-typedef float real;
+    typedef float real;
+#endif // UN_DOUBLE_PRECISION
 
-#ifndef PI
-#   define PI 3.14159265358979323846f
+#ifndef PI_D
+#   define PI_D 3.14159265358979323846
 #endif
 
-#ifndef EPSILON
-#   define EPSILON 0.000001f
+#ifndef EPSILON_D
+#   define EPSILON_D 0.000000000000001
 #endif
 
-#ifndef DEG2RAD
-#   define DEG2RAD (PI/180.0f)
+#ifndef DEG2RAD_D
+#   define DEG2RAD_D (PI/180.0)
 #endif
-
-#ifndef RAD2DEG
-#   define RAD2DEG (180.0f/PI)
+#ifndef RAD2DEG_D
+#   define RAD2DEG_D (180.0/PI)
 #endif
 
-#endif // UN_DOUBLE_PRECISION
+
+#define un_vec_unwrap2(vec) vec[0], vec[1]
+#define un_vec_unwrap3(vec) vec[0], vec[1], vec[2]
+#define un_vec_unwrap4(vec) vec[0], vec[1], vec[2], vec[3]
 
 /* ---- Memory Allocators API ---- */
 
@@ -232,7 +234,7 @@ real un_m_lerpr(real a, real b, real t);
 f32  un_m_lerpf(f32 a, f32 b, f32 t);
 f64  un_m_lerpd(f64 a, f64 b, f64 t);
 
-/* 2d */
+/* 2d f32 */
 void un_m_add2f(f32 *v, f32 *a, f32 *b);
 void un_m_sub2f(f32 *v, f32 *a, f32 *b);
 void un_m_add_scalar2f(f32 *v, f32 *a, f32 scalar);
@@ -243,14 +245,16 @@ void un_m_div_scalar2f(f32 *v, f32 *a, f32 scalar);
 void un_m_dot2f(f32 *v, f32 *a, f32 *b);
 void un_m_hadamard2f(f32 *v, f32 *a, f32 *b);
 void un_m_project2f(f32 *v, f32 *a, f32 *onto);
+void un_m_reflect2f(f32 *v, f32 *a, f32 *normal);
 void un_m_cross2f(f32 *v, f32 *a);
 
 void un_m_normalize2f(f32 *v, f32 *a);
 void un_m_magnitude2f(f32 *v, f32 *a);
 void un_m_magnitude_sqr2f(f32 *v, f32 *a);
+void un_m_distance2f(f32 *v, f32 *a, f32 *b);
 void un_m_distance_sqr2f(f32 *v, f32 *a, f32 *b);
 
-/* 3d */
+/* 3d f32 */
 void un_m_add3f(f32 *v, f32 *a, f32 *b);
 void un_m_sub3f(f32 *v, f32 *a, f32 *b);
 void un_m_add_scalar3f(f32 *v, f32 *a, f32 scalar);
@@ -261,14 +265,16 @@ void un_m_div_scalar3f(f32 *v, f32 *a, f32 scalar);
 void un_m_dot3f(f32 *v, f32 *a, f32 *b);
 void un_m_hadamard3f(f32 *v, f32 *a, f32 *b);
 void un_m_project3f(f32 *v, f32 *a, f32 *onto);
+void un_m_reflect3f(f32 *v, f32 *a, f32 *normal);
 void un_m_cross3f(f32 *v, f32 *a, f32 *b);
 
 void un_m_normalize3f(f32 *v, f32 *a);
 void un_m_magnitude3f(f32 *v, f32 *a);
 void un_m_magnitude_sqr3f(f32 *v, f32 *a);
+void un_m_distance3f(f32 *v, f32 *a, f32 *b);
 void un_m_distance_sqr3f(f32 *v, f32 *a, f32 *b);
 
-/* 4d */
+/* 4d f32 */
 void un_m_add4f(f32 *v, f32 *a, f32 *b);
 void un_m_sub4f(f32 *v, f32 *a, f32 *b);
 void un_m_add_scalar4f(f32 *v, f32 *a, f32 scalar);
@@ -279,17 +285,155 @@ void un_m_div_scalar4f(f32 *v, f32 *a, f32 scalar);
 void un_m_dot4f(f32 *v, f32 *a, f32 *b);
 void un_m_hadamard4f(f32 *v, f32 *a, f32 *b);
 void un_m_project4f(f32 *v, f32 *a, f32 *onto);
+void un_m_reflect4f(f32 *v, f32 *a, f32 *normal);
 
 void un_m_normalize4f(f32 *v, f32 *a);
 void un_m_magnitude4f(f32 *v, f32 *a);
 void un_m_magnitude_sqr4f(f32 *v, f32 *a);
+void un_m_distance4f(f32 *v, f32 *a, f32 *b);
 void un_m_distance_sqr4f(f32 *v, f32 *a, f32 *b);
 
+
+/* 2d f64 */
+void un_m_add2d(f64 *v, f64 *a, f64 *b);
+void un_m_sub2d(f64 *v, f64 *a, f64 *b);
+
+void un_m_add_scalar2d(f64 *v, f64 *a, f64 scalar);
+void un_m_sub_scalar2d(f64 *v, f64 *a, f64 scalar);
+void un_m_mul_scalar2d(f64 *v, f64 *a, f64 scalar);
+void un_m_div_scalar2d(f64 *v, f64 *a, f64 scalar);
+
+void un_m_dot2d(f64 *v, f64 *a, f64 *b);
+void un_m_hadamard2d(f64 *v, f64 *a, f64 *b);
+void un_m_project2d(f64 *v, f64 *a, f64 *onto);
+void un_m_reflect2d(f64 *v, f64 *a, f64 *normal);
+void un_m_cross2d(f64 *v, f64 *a);
+
+void un_m_normalize2d(f64 *v, f64 *a);
+void un_m_magnitude2d(f64 *v, f64 *a);
+void un_m_magnitude_sqr2d(f64 *v, f64 *a);
+void un_m_distance2d(f64 *v, f64 *a, f64 *b);
+void un_m_distance_sqr2d(f64 *v, f64 *a, f64 *b);
+
+/* 3d f64 */
+void un_m_add3d(f64 *v, f64 *a, f64 *b);
+void un_m_sub3d(f64 *v, f64 *a, f64 *b);
+
+void un_m_add_scalar3d(f64 *v, f64 *a, f64 scalar);
+void un_m_sub_scalar3d(f64 *v, f64 *a, f64 scalar);
+void un_m_mul_scalar3d(f64 *v, f64 *a, f64 scalar);
+void un_m_div_scalar3d(f64 *v, f64 *a, f64 scalar);
+
+void un_m_dot3d(f64 *v, f64 *a, f64 *b);
+void un_m_hadamard3d(f64 *v, f64 *a, f64 *b);
+void un_m_project3d(f64 *v, f64 *a, f64 *onto);
+void un_m_reflect3d(f64 *v, f64 *a, f64 *normal);
+void un_m_cross3d(f64 *v, f64 *a, f64 *b);
+
+void un_m_normalize3d(f64 *v, f64 *a);
+void un_m_magnitude3d(f64 *v, f64 *a);
+void un_m_magnitude_sqr3d(f64 *v, f64 *a);
+void un_m_distance3d(f64 *v, f64 *a, f64 *b);
+void un_m_distance_sqr3d(f64 *v, f64 *a, f64 *b);
+
+/* 4d f64 */
+void un_m_add4d(f64 *v, f64 *a, f64 *b);
+void un_m_sub4d(f64 *v, f64 *a, f64 *b);
+
+void un_m_add_scalar4d(f64 *v, f64 *a, f64 scalar);
+void un_m_sub_scalar4d(f64 *v, f64 *a, f64 scalar);
+void un_m_mul_scalar4d(f64 *v, f64 *a, f64 scalar);
+void un_m_div_scalar4d(f64 *v, f64 *a, f64 scalar);
+
+void un_m_dot4d(f64 *v, f64 *a, f64 *b);
+void un_m_hadamard4d(f64 *v, f64 *a, f64 *b);
+void un_m_project4d(f64 *v, f64 *a, f64 *onto);
+void un_m_reflect4d(f64 *v, f64 *a, f64 *normal);
+
+void un_m_normalize4d(f64 *v, f64 *a);
+void un_m_magnitude4d(f64 *v, f64 *a);
+void un_m_magnitude_sqr4d(f64 *v, f64 *a);
+void un_m_distance4d(f64 *v, f64 *a, f64 *b);
+void un_m_distance_sqr4d(f64 *v, f64 *a, f64 *b);
+
+/* 2d real */
+void un_m_add2r(real *v, real *a, real *b);
+void un_m_sub2r(real *v, real *a, real *b);
+
+void un_m_add_scalar2r(real *v, real *a, real scalar);
+void un_m_sub_scalar2r(real *v, real *a, real scalar);
+void un_m_mul_scalar2r(real *v, real *a, real scalar);
+void un_m_div_scalar2r(real *v, real *a, real scalar);
+
+void un_m_dot2r(real *v, real *a, real *b);
+void un_m_hadamard2r(real *v, real *a, real *b);
+void un_m_project2r(real *v, real *a, real *onto);
+void un_m_reflect2r(real *v, real *a, real *normal);
+void un_m_cross2r(real *v, real *a);
+
+void un_m_normalize2r(real *v, real *a);
+void un_m_magnitude2r(real *v, real *a);
+void un_m_magnitude_sqr2r(real *v, real *a);
+void un_m_distance2r(real *v, real *a, real *b);
+void un_m_distance_sqr2r(real *v, real *a, real *b);
+
+/* 3d real */
+void un_m_add3r(real *v, real *a, real *b);
+void un_m_sub3r(real *v, real *a, real *b);
+
+void un_m_add_scalar3r(real *v, real *a, real scalar);
+void un_m_sub_scalar3r(real *v, real *a, real scalar);
+void un_m_mul_scalar3r(real *v, real *a, real scalar);
+void un_m_div_scalar3r(real *v, real *a, real scalar);
+
+void un_m_dot3r(real *v, real *a, real *b);
+void un_m_hadamard3r(real *v, real *a, real *b);
+void un_m_project3r(real *v, real *a, real *onto);
+void un_m_reflect3r(real *v, real *a, real *normal);
+void un_m_cross3r(real *v, real *a, real *b);
+
+void un_m_normalize3r(real *v, real *a);
+void un_m_magnitude3r(real *v, real *a);
+void un_m_magnitude_sqr3r(real *v, real *a);
+void un_m_distance3r(real *v, real *a, real *b);
+void un_m_distance_sqr3r(real *v, real *a, real *b);
+
+/* 4d real */
+void un_m_add4r(real *v, real *a, real *b);
+void un_m_sub4r(real *v, real *a, real *b);
+
+void un_m_add_scalar4r(real *v, real *a, real scalar);
+void un_m_sub_scalar4r(real *v, real *a, real scalar);
+void un_m_mul_scalar4r(real *v, real *a, real scalar);
+void un_m_div_scalar4r(real *v, real *a, real scalar);
+
+void un_m_dot4r(real *v, real *a, real *b);
+void un_m_hadamard4r(real *v, real *a, real *b);
+void un_m_project4r(real *v, real *a, real *onto);
+void un_m_reflect4r(real *v, real *a, real *normal);
+
+void un_m_normalize4r(real *v, real *a);
+void un_m_magnitude4r(real *v, real *a);
+void un_m_magnitude_sqr4r(real *v, real *a);
+void un_m_distance4r(real *v, real *a, real *b);
+void un_m_distance_sqr4r(real *v, real *a, real *b);
+
 /* ---- splines ---- */
 
-real un_m_bezierr(real a, real q0, real q1, real b, real t);
-f32  un_m_bezierf(f32 a, f32 q0, f32 q1, f32 b, f32 t);
-f64  un_m_bezierd(f64 a, f64 q0, f64 q1, f64 b, f64 t);
+void un_m_bezierr(real *v, real a, real q0, real q1, real b, real t);
+void un_m_bezier2r(real *v, real *a, real *q0, real *q1, real *b, real t);
+void un_m_bezier3r(real *v, real *a, real *q0, real *q1, real *b, real t);
+void un_m_bezier4r(real *v, real *a, real *q0, real *q1, real *b, real t);
+
+void un_m_bezierf(f32 *v, f32 a, f32 q0, f32 q1, f32 b, f32 t);
+void un_m_bezier2f(f32 *v, f32 *a, f32 *q0, f32 *q1, f32 *b, f32 t);
+void un_m_bezier3f(f32 *v, f32 *a, f32 *q0, f32 *q1, f32 *b, f32 t);
+void un_m_bezier4f(f32 *v, f32 *a, f32 *q0, f32 *q1, f32 *b, f32 t);
+
+void un_m_bezierd(f64 *v, f64 a, f64 q0, f64 q1, f64 b, f64 t);
+void un_m_bezier2d(f64 *v, f64 *a, f64 *q0, f64 *q1, f64 *b, f64 t);
+void un_m_bezier3d(f64 *v, f64 *a, f64 *q0, f64 *q1, f64 *b, f64 t);
+void un_m_bezier4d(f64 *v, f64 *a, f64 *q0, f64 *q1, f64 *b, f64 t);
 
 /* ---- easing ---- */
 
diff --git a/tests/un/math.c b/tests/un/math.c
index 6cc7cd8..eab0dcb 100644
--- a/tests/un/math.c
+++ b/tests/un/math.c
@@ -1,5 +1,37 @@
 #include <ungrateful.h>
 
-int main(void) {
+int main() {
+    // Lerping
+    {
+        f64 d;
+        real r;
+        f32 f;
+        
+        r = un_m_lerpr(1.0f, 2.0f, 0.5f);
+        assert(fabs(r - 1.5f) < EPSILON);
+        
+        f = un_m_lerpf(1.0f, 2.0f, 0.5f);
+        assert(fabsf(f - 1.5f) < EPSILON);
+        
+        d = un_m_lerpd(1.0, 2.0, 0.5);
+        assert(fabs(d - 1.5) < 0.000001);
+    }
 
+    // Test Bezier curve functions
+    {
+        f64 d;
+        real r;
+        f32 f;
+
+        un_m_bezierf(&f, 0.0f, 1.0f, 2.0f, 3.0f, 0.5f);
+        assert(fabsf(f - 1.5f) < EPSILON);
+
+        un_m_bezierd(&d, 0.0, 1.0, 2.0, 3.0, 0.5);
+        assert(fabs(d - 1.5) < EPSILON);
+
+        un_m_bezierr(&r, 0.0, 1.0, 2.0, 3.0, 0.5);
+        assert(fabs(r - 1.5) < EPSILON);
+    }
+
+    return 0;
 }
diff --git a/tests/un/mathd.c b/tests/un/mathd.c
new file mode 100644
index 0000000..3381dfb
--- /dev/null
+++ b/tests/un/mathd.c
@@ -0,0 +1,280 @@
+#include <ungrateful.h>
+
+int main() {
+    /* 2d */
+    {
+        f64 a[2] = {1.0, 2.0};
+        f64 b[2] = {3.0, 4.0};
+        f64 n[2] = {0.0, 1.0};
+        f64 v[2];
+
+        un_m_add2d(v, a, b);
+        assert(fabs(v[0] - 4.0) < EPSILON);
+        assert(fabs(v[1] - 6.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_sub2d(v, a, b);
+        assert(fabs(v[0] + 2.0) < EPSILON);
+        assert(fabs(v[1] + 2.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_add_scalar2d(v, a, 1);
+        assert(fabs(v[0] - 2.0) < EPSILON);
+        assert(fabs(v[1] - 3.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_sub_scalar2d(v, a, 1);
+        assert(fabs(v[0])        < EPSILON);
+        assert(fabs(v[1] - 1.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_mul_scalar2d(v, a, 2);
+        assert(fabs(v[0] - 2.0) < EPSILON);
+        assert(fabs(v[1] - 4.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_div_scalar2d(v, a, 2);
+        assert(fabs(v[0] - 0.5) < EPSILON);
+        assert(fabs(v[1] - 1.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_dot2d(v, a, b);
+        assert(fabs(v[0] - 11.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_hadamard2d(v, a, b);
+        assert(fabs(v[0] - 3.0) < EPSILON);
+        assert(fabs(v[1] - 8.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_project2d(v, a, b);
+        assert(fabs(v[0] - 1.32) < EPSILON);
+        assert(fabs(v[1] - 1.76) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_cross2d(v, a);
+        assert(fabs(v[0] - 2.0) < EPSILON);
+        assert(fabs(v[1] + 1.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_magnitude2d(v, b);
+        assert(fabs(v[0] - 5.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_magnitude_sqr2d(v, b);
+        assert(fabs(v[0] - 25.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_distance2d(v, a, b);
+        assert(fabs(v[0] - sqrt(8.0)) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_distance_sqr2d(v, a, b);
+        assert(fabs(v[0] - 8.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_normalize2d(v, b);
+        assert(fabs(v[0] - 0.6) < EPSILON);
+        assert(fabs(v[1] - 0.8) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_reflect2d(v, a, n);
+        assert(fabs(v[0] + a[0]) < EPSILON);
+        assert(fabs(v[1] - a[1]) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+    }
+
+    /* 3d */
+    {
+        f64 a[3] = {1.0, 2.0, 3.0};
+        f64 b[3] = {3.0, 4.0, 5.0};
+        f64 n[3] = {0.0, 1.0, 0.0};
+        f64 v[3];
+
+        un_m_add3d(v, a, b);
+        assert(fabs(v[0] - 4.0) < EPSILON);
+        assert(fabs(v[1] - 6.0) < EPSILON);
+        assert(fabs(v[2] - 8.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_sub3d(v, a, b);
+        assert(fabs(v[0] + 2.0) < EPSILON);
+        assert(fabs(v[1] + 2.0) < EPSILON);
+        assert(fabs(v[2] + 2.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_add_scalar3d(v, a, 1);
+        assert(fabs(v[0] - 2.0) < EPSILON);
+        assert(fabs(v[1] - 3.0) < EPSILON);
+        assert(fabs(v[2] - 4.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_sub_scalar3d(v, a, 1);
+        assert(fabs(v[0])        < EPSILON);
+        assert(fabs(v[1] - 1.0) < EPSILON);
+        assert(fabs(v[2] - 2.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_mul_scalar3d(v, a, 2);
+        assert(fabs(v[0] - 2.0)  < EPSILON);
+        assert(fabs(v[1] - 4.0)  < EPSILON);
+        assert(fabs(v[2] - 6.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_div_scalar3d(v, a, 2);
+        assert(fabs(v[0] - 0.5) < EPSILON);
+        assert(fabs(v[1] - 1.0) < EPSILON);
+        assert(fabs(v[2] - 1.5) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_dot3d(v, a, b);
+        assert(fabs(v[0] - 26.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_hadamard3d(v, a, b);
+        assert(fabs(v[0] - 3.0)  < EPSILON);
+        assert(fabs(v[1] - 8.0)  < EPSILON);
+        assert(fabs(v[2] - 15.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_project3d(v, a, b);
+        assert(fabs(v[0] - 1.56) < EPSILON);
+        assert(fabs(v[1] - 2.08) < EPSILON);
+        assert(fabs(v[2] - 2.6) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_cross3d(v, a, b);
+        assert(fabs(v[0] + 2.0) < EPSILON);
+        assert(fabs(v[1] - 4.0) < EPSILON);
+        assert(fabs(v[2] + 2.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_magnitude3d(v, b);
+        assert(fabs(v[0] - sqrt(50.0)) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_magnitude_sqr3d(v, b);
+        assert(fabs(v[0] - 50.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_distance3d(v, a, b);
+        assert(fabs(v[0] - sqrt(12.0)) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_distance_sqr3d(v, a, b);
+        assert(fabs(v[0] - 12.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_normalize3d(v, b);
+        assert(fabs(v[0] - 0.424264) < EPSILON);
+        assert(fabs(v[1] - 0.565685) < EPSILON);
+        assert(fabs(v[2] - 0.707107) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_reflect3d(v, a, n);
+        assert(fabs(v[0] + a[0]) < EPSILON);
+        assert(fabs(v[1] - a[1]) < EPSILON);
+        assert(fabs(v[2] + a[2]) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+    }
+
+    /* 4d */
+    {
+        f64 a[4] = {1.0, 2.0, 3.0, 4.0};
+        f64 b[4] = {3.0, 4.0, 5.0, 6.0};
+        f64 n[4] = {0.0, 1.0, 0.0, 0.0};
+        f64 v[4];
+
+        un_m_add4d(v, a, b);
+        assert(fabs(v[0] - 4.0) < EPSILON);
+        assert(fabs(v[1] - 6.0) < EPSILON);
+        assert(fabs(v[2] - 8.0) < EPSILON);
+        assert(fabs(v[3] - 10.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_sub4d(v, a, b);
+        assert(fabs(v[0] + 2.0) < EPSILON);
+        assert(fabs(v[1] + 2.0) < EPSILON);
+        assert(fabs(v[2] + 2.0) < EPSILON);
+        assert(fabs(v[3] + 2.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_add_scalar4d(v, a, 1);
+        assert(fabs(v[0] - 2.0) < EPSILON);
+        assert(fabs(v[1] - 3.0) < EPSILON);
+        assert(fabs(v[2] - 4.0) < EPSILON);
+        assert(fabs(v[3] - 5.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_sub_scalar4d(v, a, 1);
+        assert(fabs(v[0])        < EPSILON);
+        assert(fabs(v[1] - 1.0) < EPSILON);
+        assert(fabs(v[2] - 2.0) < EPSILON);
+        assert(fabs(v[3] - 3.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_mul_scalar4d(v, a, 2);
+        assert(fabs(v[0] - 2.0)  < EPSILON);
+        assert(fabs(v[1] - 4.0)  < EPSILON);
+        assert(fabs(v[2] - 6.0) < EPSILON);
+        assert(fabs(v[3] - 8.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_div_scalar4d(v, a, 2);
+        assert(fabs(v[0] - 0.5) < EPSILON);
+        assert(fabs(v[1] - 1.0) < EPSILON);
+        assert(fabs(v[2] - 1.5) < EPSILON);
+        assert(fabs(v[3] - 2.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_dot4d(v, a, b);
+        assert(fabs(v[0] - 50.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_hadamard4d(v, a, b);
+        assert(fabs(v[0] - 3.0)  < EPSILON);
+        assert(fabs(v[1] - 8.0)  < EPSILON);
+        assert(fabs(v[2] - 15.0) < EPSILON);
+        assert(fabs(v[3] - 24.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_project4d(v, a, b);
+        assert(fabs(v[0] - 1.744186) < EPSILON);
+        assert(fabs(v[1] - 2.325581) < EPSILON);
+        assert(fabs(v[2] - 2.906976) < EPSILON);
+        assert(fabs(v[3] - 3.488372) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_magnitude4d(v, b);
+        assert(fabs(v[0] - sqrt(86.0)) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_magnitude_sqr4d(v, b);
+        assert(fabs(v[0] - 86.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_distance4d(v, a, b);
+        assert(fabs(v[0] - 4.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_distance_sqr4d(v, a, b);
+        assert(fabs(v[0] - 16.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_normalize4d(v, b);
+        assert(fabs(v[0] - 0.323498) < EPSILON);
+        assert(fabs(v[1] - 0.431331) < EPSILON);
+        assert(fabs(v[2] - 0.539164) < EPSILON);
+        assert(fabs(v[3] - 0.646997) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_reflect4d(v, a, n);
+        assert(fabs(v[0] + a[0]) < EPSILON);
+        assert(fabs(v[1] - a[1]) < EPSILON);
+        assert(fabs(v[2] + a[2]) < EPSILON);
+        assert(fabs(v[3] + a[3]) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+    }
+
+    return 0;
+}
diff --git a/tests/un/mathf.c b/tests/un/mathf.c
new file mode 100644
index 0000000..cd3a24f
--- /dev/null
+++ b/tests/un/mathf.c
@@ -0,0 +1,280 @@
+#include <ungrateful.h>
+
+int main() {
+    /* 2d */
+    {
+        f32 a[2] = {1.0f, 2.0f};
+        f32 b[2] = {3.0f, 4.0f};
+        f32 n[2] = {0.0f, 1.0f};
+        f32 v[2];
+
+        un_m_add2f(v, a, b);
+        assert(fabsf(v[0] - 4.0f) < EPSILON);
+        assert(fabsf(v[1] - 6.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_sub2f(v, a, b);
+        assert(fabsf(v[0] + 2.0f) < EPSILON);
+        assert(fabsf(v[1] + 2.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_add_scalar2f(v, a, 1);
+        assert(fabsf(v[0] - 2.0f) < EPSILON);
+        assert(fabsf(v[1] - 3.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_sub_scalar2f(v, a, 1);
+        assert(fabsf(v[0])        < EPSILON);
+        assert(fabsf(v[1] - 1.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_mul_scalar2f(v, a, 2);
+        assert(fabsf(v[0] - 2.0f) < EPSILON);
+        assert(fabsf(v[1] - 4.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_div_scalar2f(v, a, 2);
+        assert(fabsf(v[0] - 0.5f) < EPSILON);
+        assert(fabsf(v[1] - 1.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_dot2f(v, a, b);
+        assert(fabsf(v[0] - 11.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_hadamard2f(v, a, b);
+        assert(fabsf(v[0] - 3.0f) < EPSILON);
+        assert(fabsf(v[1] - 8.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_project2f(v, a, b);
+        assert(fabsf(v[0] - 1.32f) < EPSILON);
+        assert(fabsf(v[1] - 1.76f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_cross2f(v, a);
+        assert(fabsf(v[0] - 2.0f) < EPSILON);
+        assert(fabsf(v[1] + 1.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_magnitude2f(v, b);
+        assert(fabsf(v[0] - 5.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_magnitude_sqr2f(v, b);
+        assert(fabsf(v[0] - 25.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_distance2f(v, a, b);
+        assert(fabsf(v[0] - sqrtf(8.0f)) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_distance_sqr2f(v, a, b);
+        assert(fabsf(v[0] - 8.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_normalize2f(v, b);
+        assert(fabsf(v[0] - 0.6f) < EPSILON);
+        assert(fabsf(v[1] - 0.8f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_reflect2f(v, a, n);
+        assert(fabsf(v[0] + a[0]) < EPSILON);
+        assert(fabsf(v[1] - a[1]) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+    }
+
+    /* 3d */
+    {
+        f32 a[3] = {1.0f, 2.0f, 3.0f};
+        f32 b[3] = {3.0f, 4.0f, 5.0f};
+        f32 n[3] = {0.0f, 1.0f, 0.0f};
+        f32 v[3];
+
+        un_m_add3f(v, a, b);
+        assert(fabsf(v[0] - 4.0f) < EPSILON);
+        assert(fabsf(v[1] - 6.0f) < EPSILON);
+        assert(fabsf(v[2] - 8.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_sub3f(v, a, b);
+        assert(fabsf(v[0] + 2.0f) < EPSILON);
+        assert(fabsf(v[1] + 2.0f) < EPSILON);
+        assert(fabsf(v[2] + 2.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_add_scalar3f(v, a, 1);
+        assert(fabsf(v[0] - 2.0f) < EPSILON);
+        assert(fabsf(v[1] - 3.0f) < EPSILON);
+        assert(fabsf(v[2] - 4.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_sub_scalar3f(v, a, 1);
+        assert(fabsf(v[0])        < EPSILON);
+        assert(fabsf(v[1] - 1.0f) < EPSILON);
+        assert(fabsf(v[2] - 2.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_mul_scalar3f(v, a, 2);
+        assert(fabsf(v[0] - 2.0f)  < EPSILON);
+        assert(fabsf(v[1] - 4.0f)  < EPSILON);
+        assert(fabsf(v[2] - 6.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_div_scalar3f(v, a, 2);
+        assert(fabsf(v[0] - 0.5f) < EPSILON);
+        assert(fabsf(v[1] - 1.0f) < EPSILON);
+        assert(fabsf(v[2] - 1.5f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_dot3f(v, a, b);
+        assert(fabsf(v[0] - 26.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_hadamard3f(v, a, b);
+        assert(fabsf(v[0] - 3.0f)  < EPSILON);
+        assert(fabsf(v[1] - 8.0f)  < EPSILON);
+        assert(fabsf(v[2] - 15.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_project3f(v, a, b);
+        assert(fabsf(v[0] - 1.56f) < EPSILON);
+        assert(fabsf(v[1] - 2.08f) < EPSILON);
+        assert(fabsf(v[2] - 2.6f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_cross3f(v, a, b);
+        assert(fabsf(v[0] + 2.0f) < EPSILON);
+        assert(fabsf(v[1] - 4.0f) < EPSILON);
+        assert(fabsf(v[2] + 2.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_magnitude3f(v, b);
+        assert(fabsf(v[0] - sqrtf(50.0f)) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_magnitude_sqr3f(v, b);
+        assert(fabsf(v[0] - 50.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_distance3f(v, a, b);
+        assert(fabsf(v[0] - sqrtf(12.0f)) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_distance_sqr3f(v, a, b);
+        assert(fabsf(v[0] - 12.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_normalize3f(v, b);
+        assert(fabsf(v[0] - 0.424264f) < EPSILON);
+        assert(fabsf(v[1] - 0.565685f) < EPSILON);
+        assert(fabsf(v[2] - 0.707107f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_reflect3f(v, a, n);
+        assert(fabsf(v[0] + a[0]) < EPSILON);
+        assert(fabsf(v[1] - a[1]) < EPSILON);
+        assert(fabsf(v[2] + a[2]) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+    }
+
+    /* 4d */
+    {
+        f32 a[4] = {1.0f, 2.0f, 3.0f, 4.0f};
+        f32 b[4] = {3.0f, 4.0f, 5.0f, 6.0f};
+        f32 n[4] = {0.0f, 1.0f, 0.0f, 0.0f};
+        f32 v[4];
+
+        un_m_add4f(v, a, b);
+        assert(fabsf(v[0] - 4.0f) < EPSILON);
+        assert(fabsf(v[1] - 6.0f) < EPSILON);
+        assert(fabsf(v[2] - 8.0f) < EPSILON);
+        assert(fabsf(v[3] - 10.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_sub4f(v, a, b);
+        assert(fabsf(v[0] + 2.0f) < EPSILON);
+        assert(fabsf(v[1] + 2.0f) < EPSILON);
+        assert(fabsf(v[2] + 2.0f) < EPSILON);
+        assert(fabsf(v[3] + 2.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_add_scalar4f(v, a, 1);
+        assert(fabsf(v[0] - 2.0f) < EPSILON);
+        assert(fabsf(v[1] - 3.0f) < EPSILON);
+        assert(fabsf(v[2] - 4.0f) < EPSILON);
+        assert(fabsf(v[3] - 5.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_sub_scalar4f(v, a, 1);
+        assert(fabsf(v[0])        < EPSILON);
+        assert(fabsf(v[1] - 1.0f) < EPSILON);
+        assert(fabsf(v[2] - 2.0f) < EPSILON);
+        assert(fabsf(v[3] - 3.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_mul_scalar4f(v, a, 2);
+        assert(fabsf(v[0] - 2.0f)  < EPSILON);
+        assert(fabsf(v[1] - 4.0f)  < EPSILON);
+        assert(fabsf(v[2] - 6.0f) < EPSILON);
+        assert(fabsf(v[3] - 8.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_div_scalar4f(v, a, 2);
+        assert(fabsf(v[0] - 0.5f) < EPSILON);
+        assert(fabsf(v[1] - 1.0f) < EPSILON);
+        assert(fabsf(v[2] - 1.5f) < EPSILON);
+        assert(fabsf(v[3] - 2.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_dot4f(v, a, b);
+        assert(fabsf(v[0] - 50.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_hadamard4f(v, a, b);
+        assert(fabsf(v[0] - 3.0f)  < EPSILON);
+        assert(fabsf(v[1] - 8.0f)  < EPSILON);
+        assert(fabsf(v[2] - 15.0f) < EPSILON);
+        assert(fabsf(v[3] - 24.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_project4f(v, a, b);
+        assert(fabsf(v[0] - 1.744186f) < EPSILON);
+        assert(fabsf(v[1] - 2.325581f) < EPSILON);
+        assert(fabsf(v[2] - 2.906976f) < EPSILON);
+        assert(fabsf(v[3] - 3.488372f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_magnitude4f(v, b);
+        assert(fabsf(v[0] - sqrtf(86.0f)) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_magnitude_sqr4f(v, b);
+        assert(fabsf(v[0] - 86.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_distance4f(v, a, b);
+        assert(fabsf(v[0] - 4.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_distance_sqr4f(v, a, b);
+        assert(fabsf(v[0] - 16.0f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_normalize4f(v, b);
+        assert(fabsf(v[0] - 0.323498f) < EPSILON);
+        assert(fabsf(v[1] - 0.431331f) < EPSILON);
+        assert(fabsf(v[2] - 0.539164f) < EPSILON);
+        assert(fabsf(v[3] - 0.646997f) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_reflect4f(v, a, n);
+        assert(fabsf(v[0] + a[0]) < EPSILON);
+        assert(fabsf(v[1] - a[1]) < EPSILON);
+        assert(fabsf(v[2] + a[2]) < EPSILON);
+        assert(fabsf(v[3] + a[3]) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+    }
+
+    return 0;
+}
diff --git a/tests/un/mathr.c b/tests/un/mathr.c
new file mode 100644
index 0000000..05653b4
--- /dev/null
+++ b/tests/un/mathr.c
@@ -0,0 +1,280 @@
+#include <ungrateful.h>
+
+int main() {
+    /* 2d */
+    {
+        real a[2] = {1.0, 2.0};
+        real b[2] = {3.0, 4.0};
+        real n[2] = {0.0, 1.0};
+        real v[2];
+
+        un_m_add2r(v, a, b);
+        assert(fabs(v[0] - 4.0) < EPSILON);
+        assert(fabs(v[1] - 6.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_sub2r(v, a, b);
+        assert(fabs(v[0] + 2.0) < EPSILON);
+        assert(fabs(v[1] + 2.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_add_scalar2r(v, a, 1);
+        assert(fabs(v[0] - 2.0) < EPSILON);
+        assert(fabs(v[1] - 3.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_sub_scalar2r(v, a, 1);
+        assert(fabs(v[0])        < EPSILON);
+        assert(fabs(v[1] - 1.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_mul_scalar2r(v, a, 2);
+        assert(fabs(v[0] - 2.0) < EPSILON);
+        assert(fabs(v[1] - 4.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_div_scalar2r(v, a, 2);
+        assert(fabs(v[0] - 0.5) < EPSILON);
+        assert(fabs(v[1] - 1.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_dot2r(v, a, b);
+        assert(fabs(v[0] - 11.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_hadamard2r(v, a, b);
+        assert(fabs(v[0] - 3.0) < EPSILON);
+        assert(fabs(v[1] - 8.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_project2r(v, a, b);
+        assert(fabs(v[0] - 1.32) < EPSILON);
+        assert(fabs(v[1] - 1.76) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_cross2r(v, a);
+        assert(fabs(v[0] - 2.0) < EPSILON);
+        assert(fabs(v[1] + 1.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_magnitude2r(v, b);
+        assert(fabs(v[0] - 5.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_magnitude_sqr2r(v, b);
+        assert(fabs(v[0] - 25.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_distance2r(v, a, b);
+        assert(fabs(v[0] - sqrt(8.0)) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_distance_sqr2r(v, a, b);
+        assert(fabs(v[0] - 8.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_normalize2r(v, b);
+        assert(fabs(v[0] - 0.6) < EPSILON);
+        assert(fabs(v[1] - 0.8) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+
+        un_m_reflect2r(v, a, n);
+        assert(fabs(v[0] + a[0]) < EPSILON);
+        assert(fabs(v[1] - a[1]) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 2);
+    }
+
+    /* 3d */
+    {
+        real a[3] = {1.0, 2.0, 3.0};
+        real b[3] = {3.0, 4.0, 5.0};
+        real n[3] = {0.0, 1.0, 0.0};
+        real v[3];
+
+        un_m_add3r(v, a, b);
+        assert(fabs(v[0] - 4.0) < EPSILON);
+        assert(fabs(v[1] - 6.0) < EPSILON);
+        assert(fabs(v[2] - 8.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_sub3r(v, a, b);
+        assert(fabs(v[0] + 2.0) < EPSILON);
+        assert(fabs(v[1] + 2.0) < EPSILON);
+        assert(fabs(v[2] + 2.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_add_scalar3r(v, a, 1);
+        assert(fabs(v[0] - 2.0) < EPSILON);
+        assert(fabs(v[1] - 3.0) < EPSILON);
+        assert(fabs(v[2] - 4.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_sub_scalar3r(v, a, 1);
+        assert(fabs(v[0])        < EPSILON);
+        assert(fabs(v[1] - 1.0) < EPSILON);
+        assert(fabs(v[2] - 2.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_mul_scalar3r(v, a, 2);
+        assert(fabs(v[0] - 2.0)  < EPSILON);
+        assert(fabs(v[1] - 4.0)  < EPSILON);
+        assert(fabs(v[2] - 6.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_div_scalar3r(v, a, 2);
+        assert(fabs(v[0] - 0.5) < EPSILON);
+        assert(fabs(v[1] - 1.0) < EPSILON);
+        assert(fabs(v[2] - 1.5) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_dot3r(v, a, b);
+        assert(fabs(v[0] - 26.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_hadamard3r(v, a, b);
+        assert(fabs(v[0] - 3.0)  < EPSILON);
+        assert(fabs(v[1] - 8.0)  < EPSILON);
+        assert(fabs(v[2] - 15.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_project3r(v, a, b);
+        assert(fabs(v[0] - 1.56) < EPSILON);
+        assert(fabs(v[1] - 2.08) < EPSILON);
+        assert(fabs(v[2] - 2.6) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_cross3r(v, a, b);
+        assert(fabs(v[0] + 2.0) < EPSILON);
+        assert(fabs(v[1] - 4.0) < EPSILON);
+        assert(fabs(v[2] + 2.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_magnitude3r(v, b);
+        assert(fabs(v[0] - sqrt(50.0)) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_magnitude_sqr3r(v, b);
+        assert(fabs(v[0] - 50.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_distance3r(v, a, b);
+        assert(fabs(v[0] - sqrt(12.0)) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_distance_sqr3r(v, a, b);
+        assert(fabs(v[0] - 12.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_normalize3r(v, b);
+        assert(fabs(v[0] - 0.424264) < EPSILON);
+        assert(fabs(v[1] - 0.565685) < EPSILON);
+        assert(fabs(v[2] - 0.707107) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+
+        un_m_reflect3r(v, a, n);
+        assert(fabs(v[0] + a[0]) < EPSILON);
+        assert(fabs(v[1] - a[1]) < EPSILON);
+        assert(fabs(v[2] + a[2]) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 3);
+    }
+
+    /* 4d */
+    {
+        real a[4] = {1.0, 2.0, 3.0, 4.0};
+        real b[4] = {3.0, 4.0, 5.0, 6.0};
+        real n[4] = {0.0, 1.0, 0.0, 0.0};
+        real v[4];
+
+        un_m_add4r(v, a, b);
+        assert(fabs(v[0] - 4.0) < EPSILON);
+        assert(fabs(v[1] - 6.0) < EPSILON);
+        assert(fabs(v[2] - 8.0) < EPSILON);
+        assert(fabs(v[3] - 10.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_sub4r(v, a, b);
+        assert(fabs(v[0] + 2.0) < EPSILON);
+        assert(fabs(v[1] + 2.0) < EPSILON);
+        assert(fabs(v[2] + 2.0) < EPSILON);
+        assert(fabs(v[3] + 2.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_add_scalar4r(v, a, 1);
+        assert(fabs(v[0] - 2.0) < EPSILON);
+        assert(fabs(v[1] - 3.0) < EPSILON);
+        assert(fabs(v[2] - 4.0) < EPSILON);
+        assert(fabs(v[3] - 5.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_sub_scalar4r(v, a, 1);
+        assert(fabs(v[0])        < EPSILON);
+        assert(fabs(v[1] - 1.0) < EPSILON);
+        assert(fabs(v[2] - 2.0) < EPSILON);
+        assert(fabs(v[3] - 3.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_mul_scalar4r(v, a, 2);
+        assert(fabs(v[0] - 2.0)  < EPSILON);
+        assert(fabs(v[1] - 4.0)  < EPSILON);
+        assert(fabs(v[2] - 6.0) < EPSILON);
+        assert(fabs(v[3] - 8.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_div_scalar4r(v, a, 2);
+        assert(fabs(v[0] - 0.5) < EPSILON);
+        assert(fabs(v[1] - 1.0) < EPSILON);
+        assert(fabs(v[2] - 1.5) < EPSILON);
+        assert(fabs(v[3] - 2.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_dot4r(v, a, b);
+        assert(fabs(v[0] - 50.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_hadamard4r(v, a, b);
+        assert(fabs(v[0] - 3.0)  < EPSILON);
+        assert(fabs(v[1] - 8.0)  < EPSILON);
+        assert(fabs(v[2] - 15.0) < EPSILON);
+        assert(fabs(v[3] - 24.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_project4r(v, a, b);
+        assert(fabs(v[0] - 1.744186) < EPSILON);
+        assert(fabs(v[1] - 2.325581) < EPSILON);
+        assert(fabs(v[2] - 2.906976) < EPSILON);
+        assert(fabs(v[3] - 3.488372) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_magnitude4r(v, b);
+        assert(fabs(v[0] - sqrt(86.0)) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_magnitude_sqr4r(v, b);
+        assert(fabs(v[0] - 86.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_distance4r(v, a, b);
+        assert(fabs(v[0] - 4.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_distance_sqr4r(v, a, b);
+        assert(fabs(v[0] - 16.0) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_normalize4r(v, b);
+        assert(fabs(v[0] - 0.323498) < EPSILON);
+        assert(fabs(v[1] - 0.431331) < EPSILON);
+        assert(fabs(v[2] - 0.539164) < EPSILON);
+        assert(fabs(v[3] - 0.646997) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+
+        un_m_reflect4r(v, a, n);
+        assert(fabs(v[0] + a[0]) < EPSILON);
+        assert(fabs(v[1] - a[1]) < EPSILON);
+        assert(fabs(v[2] + a[2]) < EPSILON);
+        assert(fabs(v[3] + a[3]) < EPSILON);
+        un_memory_set((void*)v, 0, sizeof(*v) * 4);
+    }
+
+    return 0;
+}