diff --git a/TODO b/TODO
index 4681ef3..148eba9 100644
--- a/TODO
+++ b/TODO
@@ -4,10 +4,8 @@
 	017 [x] re-implement dynamic arrays (should mirror procedures on 'Static_Array')
 	
 [Jesse]
-	011 [math] add more Vec math procedures 
 	034 [math] Rect '#type,distinct Vec(4, T)' type + rectcut
 	002 [math] use simd for Mat4 operations
-	041 [math] small test suite
 
 	
 *** UP NEXT (NO PARTICULAR ORDER) ***
@@ -74,6 +72,8 @@
 
 
 *** DONE ***
+	041 [math] small test suite
+	011 [math] add more Vec math procedures
 	010 [x] make base jai files standalone (ie. [array], [memory], etc... big maybe)
 	042 [x] use absolute imports internally (ie. #import "jc/math"), update symlink/install process in README
 	000 [math] add easing procedures 
diff --git a/math/common.jai b/math/common.jai
new file mode 100644
index 0000000..c2f1f57
--- /dev/null
+++ b/math/common.jai
@@ -0,0 +1,86 @@
+
+PI :: 3.1415926;
+TAU :: 6.2831853;
+
+#if UNITS == .turns {
+   to_rad :: turn_to_rad;
+   from_rad :: rad_to_turn;
+}
+else #if UNITS == .radians {
+   to_rad :: rad_to_rad;
+   from_rad :: rad_to_rad;
+}
+else #if UNITS == .degrees {
+   to_rad :: deg_to_rad;
+   from_rad :: rad_to_deg;
+}
+
+turn_to_rad :: (turns: float64) -> float #expand #no_debug {
+   CONSTANT :float64: PI*2.0;
+   return cast(float)(CONSTANT*turns);
+}
+turn_to_rad :: (turns: float) -> float #expand #no_debug {
+   CONSTANT :float: xx PI*2.0;
+   return CONSTANT*turns;
+}
+
+deg_to_rad :: (deg: float64) -> float #expand #no_debug {
+   CONSTANT :float64: PI/180.0;
+   return cast(float)(CONSTANT*deg);
+}
+deg_to_rad :: (deg: float) -> float #expand #no_debug {
+   CONSTANT :float: xx (PI/180.0);
+   return CONSTANT*deg;
+}
+
+rad_to_rad :: (rad: float64) -> float #expand #no_debug {
+   return cast(float)rad;
+}
+rad_to_rad :: (rad: float) -> float #expand #no_debug {
+   return rad;
+}
+
+rad_to_turn :: (rad: float64) -> float #expand {
+   CONSTANT :float64: (1.0/(PI*2.0));
+   return cast(float)(CONSTANT*rad);
+}
+rad_to_turn :: (rad: float) -> float #expand {
+   CONSTANT :float: xx (1.0/(PI*2.0));
+   return CONSTANT*rad;
+}
+
+rad_to_deg :: (rad: float64) -> float #expand {
+   CONSTANT :float64: (180.0/PI);
+   return cast(float)(CONSTANT*rad);
+}
+rad_to_deg :: (rad: float) -> float #expand {
+   CONSTANT :float: xx (180.0/PI);
+   return CONSTANT*rad;
+}
+
+sin :: (ang: float) -> float #expand {
+   return math.sin(to_rad(ang));
+}
+
+cos :: (ang: float) -> float #expand {
+   return math.cos(to_rad(ang));
+}
+
+atan2 :: (y: float, x: float) -> float #expand {
+   return from_rad(math.atan2(y, x));
+}
+
+asin :: (ang: float) -> float #expand {
+   return from_rad(math.asin(ang));
+}
+
+abs :: (v: $T) -> T
+#modify { return meta.type_is_scalar(T), "Only accepting scalar types, integer and float"; } {
+   return ifx v < 0 then -v else v;
+}
+
+#scope_file
+
+// sin, cos
+meta :: #import "jc/meta";
+math :: #import "Math";
diff --git a/math/mat.jai b/math/mat.jai
index 57a8a22..66f6022 100644
--- a/math/mat.jai
+++ b/math/mat.jai
@@ -1,5 +1,24 @@
+Mat2 :: #type,distinct Vec(2 * 2, float);
 Mat4 :: #type,distinct Vec(4 * 4, float);
 
+m2 :: (_xx: float, xy: float, yx: float, yy: float) -> Mat2 #expand {
+   m: Mat2 = ---;
+   m._00 = _xx;
+   m._01 = xy;
+   m._10 = yx;
+   m._11 = yy;
+   return m;
+}
+
+m2 :: (r0: Vec2, r1: Vec2) -> Mat2 #expand {
+   return .{
+      components = .[
+         r0.x, r0.y,
+         r1.x, r1.y,
+      ]
+   };
+}
+
 m4 :: (r0: Vec4, r1: Vec4, r2: Vec4, r3: Vec4) -> Mat4 #expand {
    return .{
       components = .[
@@ -11,17 +30,66 @@ m4 :: (r0: Vec4, r1: Vec4, r2: Vec4, r3: Vec4) -> Mat4 #expand {
    };
 }
 
-m4_identity :: #bake_arguments m4d(diag = 1);
+m2_identity :: Mat2.{_00=1, _11=1};
+m4_identity :: Mat4.{_00=1, _11=1, _22=1, _33=1};
+
+m2d :: (diag: float) -> Mat2 #expand {
+   r: Mat2;
+   r._00 = diag;
+   r._11 = diag;
+   return r;
+}
 
 m4d :: (diag: float) -> Mat4 #expand {
    res: Mat4;
-   res[0][0] = diag;
-   res[1][1] = diag;
-   res[2][2] = diag;
-   res[3][3] = diag;
+   res._00 = diag;
+   res._11 = diag;
+   res._22 = diag;
+   res._33 = diag;
    return res;
 }
 
+operator [] :: inline (m: Mat2, idx: int) -> Vec2 #expand {
+   bounds_check_index(idx * 2, m.N);
+   return (m.components[idx * 2]).(*Vec4).*;
+}
+
+operator *[] :: inline (m: *Mat2, $$idx: int) -> *Vec2 #no_abc {
+   N :: Mat2.{}.N; // @note(judah): dumb workaround for not being able to access pointer type field constants
+   bounds_check_index(idx * 2, N);
+   return (*m.components[idx * 2]).(*Vec4);
+}
+
+operator []= :: inline (m: *Mat2, $$idx: int, value: Vec2) #no_abc {
+   N :: Mat2.{}.N;
+   bounds_check_index(idx * 4, N);
+
+   ptr := (*m.components[idx * 4]).(*Vec2);
+   ptr.*= value;
+}
+
+operator* :: inline (a: Mat2, b: Mat2) -> Mat2 {
+   r: Mat2 = ---;
+   r._00 = a._00*b._00 + a._01*b._10;
+   r._01 = a._00*b._01 + a._01*b._11;
+   r._10 = a._10*b._00 + a._11*b._10;
+   r._11 = a._10*b._01 + a._11*b._11;
+   return r;
+}
+
+operator* :: inline (a: Mat2, b: Vec2) -> Vec2 {
+   v: Vec2 = ---;
+   v.x = a._00*b.x + a._01*b.y;
+   v.y = a._10*b.x + a._11*b.y;
+   return v;
+}
+
+// Note(Jesse): Probably never used except for testing
+operator== :: (a: Mat2, b: Mat2) -> bool #no_abc {
+   for a.components if it != b.components[it_index] return false;
+   return true;
+}
+
 operator [] :: inline (m: Mat4, idx: int) -> Vec4 #expand {
    meta.check_bounds(idx * 4, m.N);
    return (m.components[idx * 4]).(*Vec4).*;
@@ -41,3 +109,294 @@ operator []= :: inline (m: *Mat4, $$idx: int, value: Vec4) #no_abc {
    ptr.*= value;
 }
 
+// Note(Jesse): Probably never used except for testing
+operator== :: (a: Mat4, b: Mat4) -> bool #no_abc {
+   for a.components if it != b.components[it_index] return false;
+   return true;
+}
+
+translate :: (from: Mat4, with: Vec3) -> Mat4 {
+   r := from;
+   r._03 += with.x*r._00 + with.y*r._01 + with.z*r._02;
+   r._13 += with.x*r._10 + with.y*r._11 + with.z*r._12;
+   r._23 += with.x*r._20 + with.y*r._21 + with.z*r._22;
+   return r;
+}
+
+translate :: (with: Vec3) -> Mat4 {
+   r: Mat4 = m4_identity;
+   r._03 = with.x;
+   r._13 = with.y;
+   r._23 = with.z;
+   return r;
+}
+
+make_look_at :: (camera: Vec3, at: Vec3, up_vector: Vec3) -> Mat4
+{
+   // -Z forward, Y up, X right
+   forward := normalize(at - camera);
+   right := normalize(cross(up_vector, forward)); // up_vector
+   up := normalize(cross(forward, right));
+   look_at: Mat4;
+   look_at._00 = right.x;
+   look_at._01 = right.y;
+   look_at._02 = right.z;
+   look_at._10 = up.x;
+   look_at._11 = up.y;
+   look_at._12 = up.z;
+   look_at._20 = -forward.x;
+   look_at._21 = -forward.y;
+   look_at._22 = -forward.z;
+   look_at._33 = 1;
+
+   look_at = inline translate(look_at, -camera);
+
+   return look_at;
+}
+
+operator* :: (a: Mat4, b: Mat4) -> Mat4 {
+   r: Mat4 = ---;
+   r._00 = a._00*b._00 + a._01*b._10 + a._02*b._20 + a._03*b._30;
+   r._01 = a._00*b._01 + a._01*b._11 + a._02*b._21 + a._03*b._31;
+   r._02 = a._00*b._02 + a._01*b._12 + a._02*b._22 + a._03*b._32;
+   r._03 = a._00*b._03 + a._01*b._13 + a._02*b._23 + a._03*b._33;
+
+   r._10 = a._10*b._00 + a._11*b._10 + a._12*b._20 + a._13*b._30;
+   r._11 = a._10*b._01 + a._11*b._11 + a._12*b._21 + a._13*b._31;
+   r._12 = a._10*b._02 + a._11*b._12 + a._12*b._22 + a._13*b._32;
+   r._13 = a._10*b._03 + a._11*b._13 + a._12*b._23 + a._13*b._33;
+
+   r._20 = a._20*b._00 + a._21*b._10 + a._22*b._20 + a._23*b._30;
+   r._21 = a._20*b._01 + a._21*b._11 + a._22*b._21 + a._23*b._31;
+   r._22 = a._20*b._02 + a._21*b._12 + a._22*b._22 + a._23*b._32;
+   r._23 = a._20*b._03 + a._21*b._13 + a._22*b._23 + a._23*b._33;
+
+   r._30 = a._30*b._00 + a._31*b._10 + a._32*b._20 + a._33*b._30;
+   r._31 = a._30*b._01 + a._31*b._11 + a._32*b._21 + a._33*b._31;
+   r._32 = a._30*b._02 + a._31*b._12 + a._32*b._22 + a._33*b._32;
+   r._33 = a._30*b._03 + a._31*b._13 + a._32*b._23 + a._33*b._33;
+   return r;
+}
+
+// Note(Jesse): If you want to make it symmetric go ahead, I usually don't
+operator* :: (a: Mat4, v: Vec4) -> Vec4 {
+   r: Vec4 = ---;
+   r.x = a._00*v.x + a._01*v.y + a._02*v.z + a._03*v.w;
+   r.y = a._10*v.x + a._11*v.y + a._12*v.z + a._13*v.w;
+   r.z = a._20*v.x + a._21*v.y + a._22*v.z + a._23*v.w;
+   r.w = a._30*v.x + a._31*v.y + a._32*v.z + a._33*v.w;
+   return r;
+}
+
+transpose :: (m: Mat2) -> Mat2 {
+   r: Mat2 = ---;
+   // Todo(Jesse): This can be one call to simd shuffle
+   r._00 = m._00;
+   r._01 = m._10;
+   r._10 = m._01;
+   r._11 = m._11;
+   return r;
+}
+
+transpose :: (m: Mat4) -> Mat4 {
+   r: Mat4 = ---;
+   r._00 = m._00;
+   r._01 = m._10;
+   r._02 = m._20;
+   r._03 = m._30;
+
+   r._10 = m._01;
+   r._11 = m._11;
+   r._12 = m._21;
+   r._13 = m._31;
+
+   r._20 = m._02;
+   r._21 = m._12;
+   r._22 = m._22;
+   r._23 = m._32;
+
+   r._30 = m._03;
+   r._31 = m._13;
+   r._32 = m._23;
+   r._33 = m._33;
+   return r;
+}
+
+// aspect is width/height
+// NDC Depth    depth_01
+// DirectX: 0-1 true
+// Vulkan:  0-1 true
+// Opengl: -1-1 false
+perspective_projection :: (aspect: float, fov: float, near: float, far: float, $$depth_01: bool) -> Mat4 {
+   r: Mat4 = ---;
+   // Y up, X right, Z- forward
+   fov_part := 1.0/tan(fov*0.5);
+   denom := (far - near);
+   r._00 = fov_part/aspect;
+   r._11 = fov_part;
+   r._22 = -(far + near)/denom;
+   r._23 = -2*far*near/denom;
+   r._32 = -1;
+   r._33 = 1;
+   if depth_01 {
+      // Makes depth range from 0 to 1
+      // To map -1,1 depth range to 0,1 we transform z as follows: z' = z * 0.5 + 0.5
+      r._22 = r._22 * 0.5 + r._32 * 0.5;
+      r._23 = r._23 * 0.5 + r._33 * 0.5;
+   }
+   return r;
+}
+
+perspective_projection_left_handed :: (aspect: float, fov: float, near: float, far: float, $$depth_01: bool) -> Mat4 {
+   r: Mat4 = ---;
+   // Y up, X right, Z- forward
+   fov_part := 1.0/tan(fov*0.5);
+   denom := (far - near);
+   r._00 = fov_part/aspect;
+   r._11 = fov_part;
+   r._22 = (far + near)/denom;
+   r._23 = 2*far*near/denom;
+   r._32 = -1;
+   r._33 = 1;
+   if depth_01 {
+      // Makes depth range from 0 to 1
+      // To map -1,1 depth range to 0,1 we transform z as follows: z' = z * 0.5 + 0.5
+      r._22 = r._22 * 0.5 + r._32 * 0.5;
+      r._23 = r._23 * 0.5 + r._33 * 0.5;
+   }
+   return r;
+}
+
+// Note(Jesse): Typically left/right/top/bottom are pixel coords
+// NDC Depth    depth_01
+// DirectX: 0-1 true
+// Vulkan:  0-1 true
+// Opengl: -1-1 false
+orthographic_projection :: (left: float, right: float, top: float, bottom: float, near: float, far: float, $$depth_01: bool) -> Mat4 {
+   r: Mat4 = m4_identity;
+   r._00 = 2.0/(right - left);
+   r._11 = 2.0/(top - bottom);
+   r._22 = -1.0/(far - near);
+   if depth_01 == false then r._22 *= 2;
+   r._30 = -(right + left)/(right - left);
+   r._31 = -(top + bottom)/(top - bottom);
+   r._32 = -(far + near)/(far - near);
+   r._33 = 1;
+   return r;
+}
+
+rotation_mat2 :: (angle: float) -> Mat2 {
+   m: Mat2 = ---;
+   m._00 =  cos(angle);
+   m._01 = -sin(angle);
+   m._10 =  sin(angle);
+   m._11 =  cos(angle);
+   return m;
+}
+
+rotation_mat4 :: (axis: Vec3, angle: float) -> Mat4 {
+   m := m4_identity;
+   a := normalize(axis);
+
+   sin_theta := sin(angle);
+   cos_theta := cos(angle);
+   inv_cos := 1.0 - cos_theta;
+
+   // Contributors to vector's X
+   m._00 = (a.x*a.x*inv_cos) + cos_theta;
+   m._01 = (a.x*a.y*inv_cos) + (a.z*sin_theta);
+   m._02 = (a.x*a.z*inv_cos) - (a.y*sin_theta);
+
+   // Contributors to vector's Y
+   m._10 = (a.y*a.x*inv_cos) - (a.z*sin_theta);
+   m._11 = (a.y*a.y*inv_cos) + cos_theta;
+   m._12 = (a.y*a.z*inv_cos) + (a.x*sin_theta);
+
+   // Contributors to vector's Z
+   m._20 = (a.z*a.x*inv_cos) + (a.y*sin_theta);
+   m._21 = (a.z*a.y*inv_cos) - (a.x*sin_theta);
+   m._22 = (a.z*a.z*inv_cos) + cos_theta;
+
+   return m;
+}
+
+rotation_mat4 :: (q: Quat) -> Mat4 {
+   m := m4_identity;
+
+   xs := 2*q.x;
+   ys := 2*q.y;
+   zs := 2*q.z;
+
+   wx := q.w*xs;
+   wy := q.w*ys;
+   wz := q.w*zs;
+
+   _xx := q.x*xs;
+   xy := q.x*ys;
+   xz := q.x*zs;
+
+   yy := q.y*ys;
+   yz := q.y*zs;
+   zz := q.z*zs;
+
+   m._00 = 1.0 - (yy + zz);
+   m._01 = xy - wz;
+   m._02 = xz + wy;
+
+   m._10 = xy + wz;
+   m._11 = 1.0 - (_xx + zz);
+   m._12 = yz - wx;
+
+   m._20 = xz - wy;
+   m._21 = yz + wx;
+   m._22 = 1.0 - (_xx + yy);
+
+   return m;
+}
+
+determinate :: (m: Mat2) -> float {
+   return m._00*m._11 - m._01*m._10;
+}
+
+determinate :: (m: Mat4) -> float {
+   c01 := cross(m.v4[0].xyz, m.v4[1].xyz);
+   c23 := cross(m.v4[2].xyz, m.v4[3].xyz);
+   b10 := m.v4[0].xyz*m.v4[1].w - m.v4[1].xyz*m.v4[0].w;
+   b32 := m.v4[2].xyz*m.v4[3].w - m.v4[3].xyz*m.v4[2].w;
+   return dot(c01, b32) + dot(c23, b10);
+}
+
+inverse :: (m: Mat2) -> Mat2 {
+   r: Mat2 = ---;
+   invdet := 1.0/determinate(m);
+   r._00 = invdet* m._11;
+   r._01 = invdet*-m._01;
+   r._10 = invdet*-m._10;
+   r._11 = invdet* m._00;
+   return r;
+}
+
+inverse :: (m: Mat4) -> Mat4 {
+   c01 := cross(m.v4[0].xyz, m.v4[1].xyz);
+   c23 := cross(m.v4[2].xyz, m.v4[3].xyz);
+   b10 := m.v4[0].xyz*m.v4[1].w - m.v4[1].xyz*m.v4[0].w;
+   b32 := m.v4[2].xyz*m.v4[3].w - m.v4[3].xyz*m.v4[2].w;
+   det := dot(c01, b32) + dot(c23, b10);
+
+   if abs(det) < 0.0001
+      return m4_identity;
+
+   inv_det := 1.0/det;
+   c01 *= inv_det;
+   c23 *= inv_det;
+   b10 *= inv_det;
+   b32 *= inv_det;
+
+   r: Mat4 = ---;
+   r.v4[0] = v4f(cross(m.v4[1].xyz, b32) + c23*m.v4[1].w, -dot(m.v4[1].xyz, c23));
+   r.v4[1] = v4f(cross(b32, m.v4[0].xyz) + c23*m.v4[0].w,  dot(m.v4[0].xyz, c23));
+   r.v4[2] = v4f(cross(m.v4[3].xyz, b10) + c01*m.v4[3].w, -dot(m.v4[3].xyz, c01));
+   r.v4[3] = v4f(cross(b10, m.v4[2].xyz) + c01*m.v4[2].w,  dot(m.v4[2].xyz, c01));
+
+   return transpose(r);
+}
diff --git a/math/module.jai b/math/module.jai
index 3938853..b1d50e2 100644
--- a/math/module.jai
+++ b/math/module.jai
@@ -21,14 +21,359 @@ F64_Min, F64_Max :: #run meta.lo_for(float64), #run meta.hi_for(float64);
 #load "vec.jai";
 #load "mat.jai";
 #load "ease.jai";
+#load "common.jai";
 
 #scope_module;
 
-#if RUN_TESTS {
+#if RUN_TESTS #run {
    test :: #import "jc/test";
+
+   vec2_tests();
+   vec3_tests();
+   vec4_tests();
+   vecn_tests();
+   mat2_tests();
+   mat4_tests();
+   quat_tests();
 }
 
 #scope_file;
 
 meta  :: #import "jc/meta";
 basic :: #import "Basic"; // @future
+
+math :: #import "Math";
+
+test :: #import "jc/test";
+
+vec2_tests :: () {
+   dot_print("Vec2 % tests", UNITS);
+   t: test.T;
+   {
+      a: Vec2 = v2f(0.0, 1.0);
+      b: Vec2 = v2f(1.0, 2.0);
+      test.expect(*t, a + b == v2f(1.0, 3.0));
+      test.expect(*t, b - a == v2f(1.0, 1.0));
+      test.expect(*t, a*b == v2f(0.0, 2.0));
+      test.expect(*t, a/b == v2f(0.0, 0.5));
+   }
+
+   {
+      a: Vec(2, int) = v2i(2, 1);
+      b: Vec(2, int) = v2i(1, 0);
+      test.expect(*t, a + b == v2i(3, 1));
+      test.expect(*t, b - a == v2i(-1, -1));
+      test.expect(*t, a*b == v2i(2, 0));
+      test.expect(*t, b/a == v2i(0, 0));
+   }
+   {
+      a: Vec2 = v2f(2.3, -4.1);
+      b: Vec2 = v2f(1.0, 3.6);
+      c := min(a, b);
+      test.expect(*t, c == v2f(1.0, -4.1));
+      c = max(a, b);
+      test.expect(*t, c == v2f(2.3, 3.6));
+   }
+}
+
+vec3_tests :: () {
+   dot_print("Vec3 % tests", UNITS);
+   t: test.T;
+   {
+      a: Vec3 = v3f(0.0, 1.0, 2.0);
+      b: Vec3 = v3f(1.0, 2.0, 3.0);
+      test.expect(*t, a + b == v3f(1.0, 3.0, 5.0));
+      test.expect(*t, b - a == v3f(1.0, 1.0, 1.0));
+      test.expect(*t, a*b == v3f(0.0, 2.0, 6.0));
+      test.expect(*t, a/b == v3f(0.0, 0.5, 0.66666667));
+
+      a = v3f(1.0, 1.0, 0.0);
+      b = v3f(1.0, 0.0, 0.0);
+      test.expect(*t, reflect(a, b) == v3f(1.0, -1.0, 0.0));
+      test.expect(*t, round(v3f(1.2, 1.7, 1.5)) == v3f(1.0, 2.0, 2.0));
+      test.expect(*t, round(v3f(-1.2, -1.7, -1.5)) == v3f(-1.0, -2.0, -2.0));
+
+      a = v3f(1.0, 0.0, 0.0);
+      b = v3f(0.0, 1.0, 0.0);
+      test.expect(*t, cross(a, b) == v3f(0.0, 0.0, 1.0));
+   }
+   {
+      a: Vec3 = v3f(2.3, 4.1, 9.0);
+      b: Vec3 = v3f(1.0, -3.6, 5.0);
+      c := min(a, b);
+      test.expect(*t, c == v3f(1.0, -3.6, 5.0));
+      c = max(a, b);
+      test.expect(*t, c == v3f(2.3, 4.1, 9.0));
+   }
+}
+
+vec4_tests :: () {
+   dot_print("Vec4 % tests", UNITS);
+   t: test.T;
+   {
+      a: Vec4 = v4f(2.25, 1.0, 2.0, 1.0);
+      b: Vec4 = v4f(4.0, 2.0, 3.0, 1.0);
+      test.expect(*t, a + b == v4f(6.25, 3.0, 5.0, 2.0));
+      test.expect(*t, b - a == v4f(1.75, 1.0, 1.0, 0.0));
+      test.expect(*t, a*b == v4f(9.0, 2.0, 6.0, 1.0));
+      test.expect(*t, a/b == v4f(0.5625, 0.5, 2.0/3.0, 1.0));
+   }
+}
+
+vecn_tests :: () {
+   dot_print("VecN % tests", UNITS);
+   t: test.T;
+   {
+      a: Vec(16, float);
+      b: Vec(16, float);
+      for *a {
+         it.* = xx it_index;
+      }
+      for *b {
+         it.* = xx(it_index + 1);
+      }
+      test.expect(*t, a + b == Vec(16, float).{.[1.0, 3.0, 5.0, 7.0, 9.0, 11.0, 13.0, 15.0, 17.0, 19.0, 21.0, 23.0, 25.0, 27.0, 29.0, 31.0]});
+      test.expect(*t, b - a == Vec(16, float).{.[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]});
+      test.expect(*t, a*b == Vec(16, float).{.[0.0, 2.0, 6.0, 12.0, 20.0, 30.0, 42.0, 56.0, 72.0, 90.0, 110.0, 132.0, 156.0, 182.0, 210.0, 240.0]});
+
+      test.expect(*t, min(a, b) == a);
+      test.expect(*t, max(a, b) == b);
+      test.expect(*t, max(a, 12) == Vec(16, float).{.[0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 12.0, 12.0, 12.0]});
+      test.expect(*t, min(a, 13.2) == Vec(16, float).{.[13.2, 13.2, 13.2, 13.2, 13.2, 13.2, 13.2, 13.2, 13.2, 13.2, 13.2, 13.2, 13.2, 13.2, 14.0, 15.0]});
+      test.expect(*t, clamp(a, 7.25, 12.0) == Vec(16, float).{.[7.25, 7.25, 7.25, 7.25, 7.25, 7.25, 7.25, 7.25, 8, 9, 10, 11, 12, 12, 12, 12]});
+
+      a1: Vec(16, float) = Vec(16, float).{.[1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 4.7, -1.2, -1.0, -1.5, 11.2, 14.0, 15.0, 14.0, 15.0, 65536.2]};
+      basic.print(">> a1:        %\n", a1);
+      basic.print(">> ceil(a1):  %\n", ceil(a1));
+      basic.print(">> floor(a1):  %\n", floor(a1));
+      test.expect(*t, ceil(a1) == Vec(16, float).{.[2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 5.0, -1.0, -1.0, -1.0, 12.0, 14.0, 15.0, 14.0, 15.0, 65537]});
+      test.expect(*t, floor(a1) == Vec(16, float).{.[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 4.0, -2.0, -1.0, -2.0, 11.0, 14.0, 15.0, 14.0, 15.0, 65536]});
+
+      test.expect(*t, dot(a, b) == 1360.0);
+      test.expect(*t, abs(a) == a);
+      c := a;
+      for *c { // Check making every other component negative
+         if it_index%2 == 0
+            it.* = -it.*;
+      }
+      test.expect(*t, abs(c) == a);
+      test.expect(*t, float_eq(length(normalize(a)), 1));
+      test.expect(*t, a + 2 == Vec(16, float).{.[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]});
+      test.expect(*t, a - 2 == Vec(16, float).{.[-2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]});
+      test.expect(*t, a*2 == Vec(16, float).{.[0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30]});
+   }
+}
+
+mat2_tests :: () {
+   dot_print("Mat2 % tests", UNITS);
+   t: test.T;
+   {
+      identity := m2_identity;
+      a: Mat2 = Mat2.{.[1, 2, 3, 4]};
+      test.expect(*t, a*identity == a);
+   }
+   {
+      basic.print("rotation matrix\n");
+      rotator := rotation_mat2(from_rad(PI));
+      v2 := v2f(1.0, 0.5);
+      basic.print("rotating:     %\n", v2);
+      v2 = rotator*v2;
+      basic.print("after rotate: %\n", v2);
+      test.expect(*t, v2_eq(v2, v2f(-1.0, -0.5)));
+      v2 = rotator*v2;
+      basic.print("rotate back:  %\n", v2);
+      test.expect(*t, v2_eq(v2, v2f(1.0, 0.5)));
+   }
+   {
+      basic.print("determinate\n");
+      m := m2(1.0, 3.0, 2.0, 2.0);
+      basic.print("m2: %\n", m);
+      det := determinate(m);
+      basic.print("determinate: %\n", det);
+      test.expect(*t, det == -4);
+   }
+   {
+      basic.print("inverse test\n");
+      rotator := rotation_mat2(from_rad(PI));
+      inv_rotator := inverse(rotator);
+      v2 := v2f(0.1, 1.0);
+      basic.print("   : %\n", rotator);
+      basic.print("inv: %\n", inv_rotator);
+
+      basic.print("inv_mat*mat = %\n", inv_rotator*rotator);
+      test.expect(*t, inv_rotator*rotator == m2_identity);
+
+      test.expect(*t, inv_rotator*(rotator*v2) == v2);
+      basic.print("v2  : %\n", v2);
+      v2 = rotator*v2;
+      basic.print("rot : %\n", v2);
+      v2 = inv_rotator*v2;
+      basic.print("inv : %\n", v2);
+   }
+}
+
+mat4_tests :: () {
+   dot_print("Mat4 % tests", UNITS);
+   t: test.T;
+   {
+      identity := m4_identity;
+      a: Mat4 = Mat4.{.[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]};
+      test.expect(*t, a*identity == a);
+   }
+   {
+      UP_VECTOR :: Vec3.{.[0.0, 1.0, 0.0]};
+      camera := v3f(1.0, 0.0, 1.0);
+      looking_at := v3f(0.0, 0.0, 0.0);
+      look_at := make_look_at(camera, looking_at, UP_VECTOR);
+      basic.print("lookat matrix: %\n", look_at);
+      basic.print("lookat from 0,1,0: %\n", look_at*v4f(0.0, 1.0, 0.0, 1.0));
+   }
+   {
+      translator := translate(v3f(10.0, 5.0, 2.0));
+      v3 := v3f(1.0, 2.0, 1.0);
+      v4 := v4f(v3, 1.0);
+      test.expect(*t, v4 == v4f(1.0, 2.0, 1.0, 1.0));
+      test.expect(*t, translator*v4 == v4f(11.0, 7.0, 3.0, 1.0));
+   }
+   {
+      basic.print("rotation matrix\n");
+      rotator := rotation_mat4(v3f(0.0, 1.0, 0.0), from_rad(PI));
+      v4 := v4f(1.0, 0.5, 0.1, 1.0);
+      basic.print("rotating:     %\n", v4);
+      v4 = rotator*v4;
+      basic.print("after rotate: %\n", v4);
+      test.expect(*t, v4_eq(v4, v4f(-1.0, 0.5, -0.1, 1.0)));
+      v4 = rotator*v4;
+      basic.print("rotate back:  %\n", v4);
+      test.expect(*t, v4_eq(v4, v4f(1.0, 0.5, 0.1, 1.0)));
+   }
+   {
+      rotator_x := rotation_mat4(v3f(1.0, 0.0, 0.0), from_rad(0.4*PI));
+      rotator_y := rotation_mat4(v3f(0.0, 1.0, 0.0), from_rad(PI));
+      camera_rotator := rotator_x*rotator_y;
+      basic.print("big rotator: %\n", camera_rotator);
+      det := determinate(camera_rotator);
+      basic.print("determinate: %\n", det);
+   }
+   {
+      basic.print("inverse test\n");
+      rotator := rotation_mat4(v3f(0.0, 0.0, 1.0), from_rad(PI));
+      inv_rotator := inverse(rotator);
+      v4 := v4f(0.1, 1.0, 0.0, 1.0);
+      basic.print("   : %\n", rotator);
+      basic.print("inv: %\n", inv_rotator);
+
+      basic.print("inv_mat*mat = %\n", inv_rotator*rotator);
+      test.expect(*t, inv_rotator*rotator == m4_identity);
+
+      basic.print("v4  : %\n", v4);
+      v4 = rotator*v4;
+      basic.print("rot : %\n", v4);
+      v4 = inv_rotator*v4;
+      basic.print("inv : %\n", v4);
+   }
+}
+
+quat_tests :: () {
+   dot_print("Quaternion % tests", UNITS);
+   t: test.T;
+   {
+      qi := quat_identity;
+      basic.print("qi lensq: %\n", length_squared(qi));
+      basic.print("qi len  : %\n", length(qi));
+      q: Quat = rotation_quat(from_rad(PI), v3f(0.0, 1.0, 0.0));
+      q2: Quat = rotation_quat(from_rad(PI), v3f(1.0, 0.0, 1.0));
+      basic.print("q       : %\n", q);
+      basic.print("q2      : %\n", q2);
+      basic.print("dot     : %\n", dot(q, q2));
+      basic.print("mul     : %\n", q*q2);
+      qc := conjugate(q);
+      basic.print("conjugate : %\n", qc);
+      inv_q := inverse(q);
+      basic.print("inverse   : %\n", inv_q);
+      test.expect(*t, q*inv_q == qi);
+
+      q1 := quat(2, 0, 0, 0);
+      q2  = quat(1, 1, -1, 0);
+      basic.print("q1: %\n", q1);
+      basic.print("q2: %\n", q2);
+      basic.print("q1*q2: %\n", q1*q2);
+      basic.print("dot(q2,q2): %\n", dot(q2, q2));
+      basic.print("dot(q1,q2): %\n", dot(q1, q2));
+      basic.print("q1*q2*conj(q1): %\n", q1*q2*conjugate(q1));
+      c := q1*q2*conjugate(q1);
+      test.expect(*t, float_eq(c.w, 0.0));
+
+      q = rotation_quat(from_rad(PI/4.0), v3f(0.0, 0.0, 1.0));
+      p := v3f(2.0, 0.0, 0.0);
+      basic.print("q: %\n", q);
+      basic.print("p: %\n", p);
+      c1 := q*quat(p, 0)*conjugate(q);
+      basic.print("q*quat(p, 0)*conjugate(q): %\n", c1);
+      c2 := q*p;
+      basic.print("q*p                      : %\n", c2);
+      test.expect(*t, v3_eq(c2, v3f(math.sqrt(2.0), math.sqrt(2.0), 0.0)));
+      test.expect(*t, v3_eq(c1.xyz, c2));
+
+      basic.print("quaternion to matrix rotation\n");
+      q = rotation_quat(from_rad(PI), v3f(0.0, 1.0, 0.0));
+      basic.print("q: %\n", q);
+      m := rotation_mat4(q);
+      print_mat4(m);
+      p1 := v4f(2.0, 0.0, 0.0, 1.0);
+      basic.print("p : %\n", p1);
+      basic.print("p': %\n", m*p1);
+      test.expect(*t, v4_eq(m*p1, v4f(-2.0, 0.0, -0.0, 1.0)));
+
+      q1 = rotation_quat(from_rad(PI), v3f(0.0, 1.0, 0.0));
+      q2 = rotation_quat(from_rad(2.0*PI), v3f(0.0, 1.0, 0.0));
+      perc := 0.5;
+      q = slerp(q1, q2, perc);
+      basic.print("q1: %\n", q1);
+      basic.print("q2: %\n", q2);
+      basic.print("slerp q1, q2, %: %\n", perc, q);
+
+      q = rotation_quat(from_rad(PI/4.0), v3f(0.0, 0.0, 1.0));
+      print_mat4(rotation_mat4(q));
+   }
+}
+
+dot_print :: (fmt: string, args: ..Any, width: int = 30) {
+   dots := ".................................................................";
+   str := basic.tprint(fmt, args);
+   dots.count = width - str.count;
+   basic.print("%0%", str, dots);
+   basic.print("\n");
+}
+
+check :: (cond: bool, loc := #caller_location) -> bool {
+   if !cond {
+      basic.print("\tTest failed at %:%\n", loc.fully_pathed_filename, loc.line_number);
+   }
+   return cond;
+}
+
+v2_eq :: (a: Vec2, b: Vec2) -> bool {
+   return float_eq(a.x, b.x) && float_eq(a.y, b.y);
+}
+
+v3_eq :: (a: Vec3, b: Vec3) -> bool {
+   return float_eq(a.x, b.x) && float_eq(a.y, b.y) && float_eq(a.z, b.z);
+}
+
+v4_eq :: (a: Vec4, b: Vec4) -> bool {
+   return float_eq(a.x, b.x) && float_eq(a.y, b.y) && float_eq(a.z, b.z) && float_eq(a.w, b.w);
+}
+
+// Smallest difference where a float is basically that value
+EPSILON :: 0.001;
+float_eq :: (f: float, with: float) -> bool {
+   return f > with - EPSILON && f < with + EPSILON;
+}
+
+print_mat4 :: (m: Mat4) {
+   basic.print("| % % % % |\n", m._00, m._01, m._02, m._03);
+   basic.print("| % % % % |\n", m._10, m._11, m._12, m._13);
+   basic.print("| % % % % |\n", m._20, m._21, m._22, m._23);
+   basic.print("| % % % % |\n", m._30, m._31, m._32, m._33);
+}
diff --git a/math/vec.jai b/math/vec.jai
index 192b63c..0753c63 100644
--- a/math/vec.jai
+++ b/math/vec.jai
@@ -66,6 +66,38 @@ Vec :: struct(N: int, T: Type)
          }
       }
 
+      // compound accessors
+      basic.append(*b, "#place components;\n");
+      if N >= 4 {
+         basic.append(*b, "union {\n");
+         basic.append(*b, "\txy:  Vec2 = ---;\n");
+         basic.append(*b, "\txyz: Vec3 = ---;\n");
+         basic.append(*b, "};");
+      }
+      else if N >= 3 {
+         basic.append(*b, "xy: T = ---;\n");
+      }
+
+      // Vec4 row accessors (Possibly becoming/used as SIMD lanes)
+      if N/4 > 1 {
+         basic.append(*b, "#place components;\n");
+         basic.print_to_builder(*b, "v4: [%]Vec4 = ---;\n", N/4);
+      }
+
+      // Matrix Accessors
+      if N == 4 { // Mat2
+         basic.append(*b, "#place components;\n");
+         for i: 0..3 {
+            basic.print_to_builder(*b, "_%0%: T = ---;\n", i/2, i%2);
+         }
+      }
+      if N == 16 { // Mat4
+         basic.append(*b, "#place components;\n");
+         for i: 0..15 {
+            basic.print_to_builder(*b, "_%0%: T = ---;\n", i/4, i%4);
+         }
+      }
+
       return basic.builder_to_string(*b);
    };
 }
@@ -100,118 +132,268 @@ for_expansion :: (v: *Vec, body: Code, flags: For_Flags) #expand {
 
 operator + :: inline (l: Vec, r: Vec(l.N, l.T)) -> Vec(l.N, l.T) #no_abc {
    res: Vec(l.N, l.T) = ---;
-   // @todo(judah): unroll for N <= 4
-   for l res[it_index] = it + r[it_index];
+   #if l.N <= 4 {
+      res.x = l.x + r.x;
+      res.y = l.y + r.y;
+      #if l.N >= 3 then res.z = l.z + r.z;
+      #if l.N == 4 then res.w = l.w + r.w; // @todo(jesse): SIMD
+   }
+   else {
+      for l res[it_index] = it + r[it_index];
+   }
    return res;
 }
 
 operator - :: inline (l: Vec, r: Vec(l.N, l.T)) -> Vec(l.N, l.T) #no_abc {
    res: Vec(l.N, l.T) = ---;
-   // @todo(judah): unroll for N <= 4
-   for l res[it_index] = it - r[it_index];
+   #if l.N <= 4 {
+      res.x = l.x - r.x;
+      res.y = l.y - r.y;
+      #if l.N >= 3 then res.z = l.z - r.z;
+      #if l.N == 4 then res.w = l.w - r.w; // @todo(jesse): SIMD
+   }
+   else {
+      for l res[it_index] = it - r[it_index];
+   }
    return res;
 }
 
 operator * :: inline (l: Vec, r: Vec(l.N, l.T)) -> Vec(l.N, l.T) #no_abc {
    res: Vec(l.N, l.T) = ---;
-   // @todo(judah): unroll for N <= 4
-   for l res[it_index] = it * r[it_index];
+   #if l.N <= 4 {
+      res.x = l.x*r.x;
+      res.y = l.y*r.y;
+      #if l.N >= 3 then res.z = l.z*r.z;
+      #if l.N == 4 then res.w = l.w*r.w; // @todo(jesse): SIMD
+   }
+   else {
+      for l res[it_index] = it * r[it_index];
+   }
    return res;
 }
 
 operator / :: inline (l: Vec, r: Vec(l.N, l.T)) -> Vec(l.N, l.T) #no_abc {
    res: Vec(l.N, l.T) = ---;
-   // @todo(judah): unroll for N <= 4
-   for l res[it_index] = it / r[it_index];
+   #if l.N <= 4 {
+      res.x = l.x/r.x;
+      res.y = l.y/r.y;
+      #if l.N >= 3 then res.z = l.z/r.z;
+      #if l.N == 4 then res.w = l.w/r.w; // @todo(jesse): SIMD
+   }
+   else {
+      for l res[it_index] = it / r[it_index];
+   }
    return res;
 }
 
 operator + :: inline (l: Vec, r: $R) -> Vec(l.N, l.T) #no_abc #symmetric
 #modify { return meta.type_is_scalar(R), "type is not integer or float"; } {
    res: Vec(l.N, l.T) = ---;
-   for l res[it_index] = it + r;
+   #if l.N <= 4 {
+      res.x = l.x + r;
+      res.y = l.y + r;
+      #if l.N >= 3 then res.z = l.z + r;
+      #if l.N == 4 then res.w = l.w + r; // @todo(jesse): SIMD
+   }
+   else {
+      for l res[it_index] = it + r;
+   }
    return res;
 }
 
 operator - :: inline (l: Vec, r: $R) -> Vec(l.N, l.T) #no_abc
 #modify { return meta.type_is_scalar(R), "type is not integer or float"; } {
    res: Vec(l.N, l.T) = ---;
-   for l res[it_index] = it - r;
+   #if l.N <= 4 {
+      res.x = l.x - r;
+      res.y = l.y - r;
+      #if l.N >= 3 then res.z = l.z - r;
+      #if l.N == 4 then res.w = l.w - r; // @todo(jesse): SIMD
+   }
+   else {
+      for l res[it_index] = it - r;
+   }
    return res;
 }
 operator - :: inline (l: $R, r: Vec) -> Vec(l.N, l.T) #no_abc
 #modify { return meta.type_is_scalar(R), "type is not integer or float"; } {
    res: Vec(l.N, l.T) = ---;
-   for l res[it_index] = r - it;
+   #if l.N <= 4 {
+      res.x = l - r.x;
+      res.y = l - r.y;
+      #if r.N >= 3 then res.z = l - r.z;
+      #if r.N == 4 then res.w = l - r.w; // @todo(jesse): SIMD
+   }
+   else {
+      for l res[it_index] = r - it;
+   }
+   return res;
+}
+
+operator- :: inline(v: Vec) -> Vec(v.N, v.T) #no_abc {
+   res: Vec(v.N, v.T) = ---;
+   #if v.N <= 4 {
+      res.x = -v.x;
+      res.y = -v.y;
+      #if v.N >= 3 then res.z = -v.z;
+      #if v.N == 4 then res.w = -v.w; // @todo(jesse): SIMD
+   }
+   else {
+      for v res[it_index] = -it;
+   }
    return res;
 }
 
 operator * :: inline (l: Vec, r: $R) -> Vec(l.N, l.T) #no_abc #symmetric
 #modify { return meta.type_is_scalar(R), "type is not integer or float"; } {
    res: Vec(l.N, l.T) = ---;
-   for l res[it_index] = it*r;
+   #if l.N <= 4 {
+      res.x = l.x*r;
+      res.y = l.y*r;
+      #if l.N >= 3 then res.z = l.z*r;
+      #if l.N == 4 then res.w = l.w*r; // @todo(jesse): SIMD
+   }
+   else {
+      for l res[it_index] = it*r;
+   }
    return res;
 }
 
 operator / :: inline (l: Vec, r: $R) -> Vec(l.N, l.T) #no_abc
 #modify { return meta.type_is_scalar(R), "type is not integer or float"; } {
    res: Vec(l.N, l.T) = ---;
-   for l res[it_index] = it/r;
+   #if l.N <= 4 {
+      res.x = l.x/r;
+      res.y = l.y/r;
+      #if l.N >= 3 then res.z = l.z/r;
+      #if l.N == 4 then res.w = l.w/r; // @todo(jesse): SIMD
+   }
+   else {
+      for l res[it_index] = it/r;
+   }
    return res;
 }
 
 operator == :: inline (l: Vec, r: Vec(l.N, l.T)) -> bool #no_abc {
-   for l if it != r[it_index] return false;
-   return true;
+   #if l.N <= 4 {
+      res: bool = l.x == r.x && l.y == r.y;
+      #if l.N >= 3 then res &= l.z == r.z;
+      #if l.N == 4 then res &= l.w == r.w; // @todo(jesse): SIMD
+      return res;
+   }
+   else {
+      for l if it != r[it_index] return false;
+      return true;
+   }
 }
 
 min :: (l: Vec, r: Vec(l.N, l.T)) -> Vec(l.N, l.T) #no_abc {
    res: Vec(l.N, l.T) = ---;
-   n := l.N - 1;
-   while n >= 0 {
-      if l[n] < r[n]
-         res[n] = l[n];
-      else
-         res[n] = r[n];
-      n -= 1;
+   #if l.N <= 4 {
+      res.x = ifx l.x < r.x then l.x else r.x;
+      res.y = ifx l.y < r.y then l.y else r.y;
+      #if l.N >= 3 then res.z = ifx l.z < r.z then l.z else r.z;
+      #if l.N == 4 then res.w = ifx l.w < r.w then l.w else r.w; // @todo(jesse): SIMD
+   }
+   else {
+      n := l.N - 1;
+      while n >= 0 {
+         if l[n] < r[n]
+            res[n] = l[n];
+         else
+            res[n] = r[n];
+         n -= 1;
+      }
    }
    return res;
 }
 
 max :: (l: Vec, r: Vec(l.N, l.T)) -> Vec(l.N, l.T) #no_abc {
    res: Vec(l.N, l.T) = ---;
-   n := l.N - 1;
-   while n >= 0 {
-      if l[n] > r[n]
-         res[n] = l[n];
-      else
-         res[n] = r[n];
-      n -= 1;
+   #if l.N <= 4 {
+      res.x = ifx l.x > r.x then l.x else r.x;
+      res.y = ifx l.y > r.y then l.y else r.y;
+      #if l.N >= 3 then res.z = ifx l.z > r.z then l.z else r.z;
+      #if l.N == 4 then res.w = ifx l.w > r.w then l.w else r.w; // @todo(jesse): SIMD
+   }
+   else {
+      n := l.N - 1;
+      while n >= 0 {
+         if l[n] > r[n]
+            res[n] = l[n];
+         else
+            res[n] = r[n];
+         n -= 1;
+      }
    }
    return res;
 }
 
-ceil :: (l: Vec, x: l.T) -> Vec(l.N, l.T) #no_abc {
+min :: (l: Vec, r: l.T) -> Vec(l.N, l.T) #no_abc {
+   res: Vec(l.N, l.T) = ---;
+   #if l.N <= 4 {
+      res.x = ifx l.x > r then l.x else r;
+      res.y = ifx l.y > r then l.y else r;
+      #if l.N >= 3 then res.z = ifx l.z > r.z then l.z else r;
+      #if l.N == 4 then res.w = ifx l.w > r.w then l.w else r; // @todo(jesse): SIMD
+   }
+   else {
+      n := l.N - 1;
+      while n >= 0 {
+         if l[n] > r
+            res[n] = l[n];
+         else
+            res[n] = r;
+         n -= 1;
+      }
+   }
+   return res;
+}
+
+max :: (l: Vec, r: l.T) -> Vec(l.N, l.T) #no_abc {
+   res: Vec(l.N, l.T) = ---;
+   #if l.N <= 4 {
+      res.x = ifx l.x < r then l.x else r;
+      res.y = ifx l.y < r then l.y else r;
+      #if l.N >= 3 then res.z = ifx l.z < r then l.z else r;
+      #if l.N == 4 then res.w = ifx l.w < r then l.w else r; // @todo(jesse): SIMD
+   }
+   else {
+      n := l.N - 1;
+      while n >= 0 {
+         if l[n] < r
+            res[n] = l[n];
+         else
+            res[n] = r;
+         n -= 1;
+      }
+   }
+   return res;
+}
+
+// @todo(jesse): SIMD
+ceil :: (l: Vec) -> Vec(l.N, l.T) #no_abc {
    r: Vec(l.N, l.T) = ---;
    n := l.N - 1;
+   sign: float;
    while n >= 0 {
-      if x < l[n]
-         r[n] = x;
-      else
-         r[n] = l[n];
+      int_part := l[n].(int);
+      add_part := ifx l[n] > int_part.(float) then 1.0 else 0.0;
+      r[n] = int_part + add_part;
       n -= 1;
    }
    return r;
 }
 
-floor :: (l: Vec, x: l.T) -> Vec(l.N, l.T) #no_abc {
+// @todo(jesse): SIMD
+floor :: (l: Vec) -> Vec(l.N, l.T) #no_abc {
    r: Vec(l.N, l.T) = ---;
    n := l.N - 1;
+   sign: float;
    while n >= 0 {
-      if x > l[n]
-         r[n] = x;
-      else
-         r[n] = l[n];
+      int_part := l[n].(int);
+      sub_part := ifx l[n] < int_part.(float) then 1.0 else 0.0;
+      r[n] = int_part - sub_part;
       n -= 1;
    }
    return r;
@@ -243,11 +425,11 @@ dot :: (a: Vec, b: Vec(a.N, a.T)) -> a.T #no_abc {
 }
 
 length_squared :: (v: Vec) -> float #no_abc {
-   return dot(v, v);
+   return inline dot(v, v);
 }
 
 length :: (v: Vec) -> float #no_abc {
-   return math.sqrt(dot(v, v));
+   return math.sqrt(inline dot(v, v));
 }
 
 abs :: (v: Vec) -> Vec(v.N, v.T) #no_abc {
@@ -263,9 +445,11 @@ abs :: (v: Vec) -> Vec(v.N, v.T) #no_abc {
    return r;
 }
 
+// @todo(Jesse): SIMD
 norm :: normalize;
 normalize :: (v: Vec) -> Vec(v.N, v.T) #no_abc {
-   return v/length(v);
+   inv_len := 1.0/length(v);
+   return inv_len*v;
 }
 
 lerp :: (a: Vec, b: Vec(a.N, a.T), t: float) -> Vec(a.N, a.T) #no_abc {
@@ -307,7 +491,7 @@ round :: (v: Vec($N, $T)) -> Vec(N, T) #no_abc
 Vec2 :: Vec(2, float);
 Vec3 :: Vec(3, float);
 Vec4 :: Vec(4, float);
-Quat :: Vec4;
+Quat :: #type,distinct Vec4; // Note(Jesse): I Had to make this distinct, otherwise operators stomp on eachother
 
 v2f :: (x: $T = 0, y: T = 0) -> Vec2
 #modify { return meta.type_is_float(T), "use v2i for integer arguments"; }
@@ -339,18 +523,122 @@ v4f :: (x: $T = 0, y: T = 0, z: T = 0, w: T = 0) -> Vec4
    return .{ x = x, y = y, z = z, w = w };
 }
 
+v4f :: (v: Vec3, $$w: float) -> Vec4 #expand {
+   return .{ xyz=v, w=w };
+}
+
 v4i :: (x: $T = 0, y: T = 0, z: T = 0, w: T = 0) -> Vec(4, T)
 #modify { return meta.type_is_integer(T), "use v4f for float arguments"; }
 #expand {
    return .{ x = x, y = y, z = z, w = w };
 }
 
-quat :: (x: float = 0, y: float = 0, z: float = 0, w: float = 0) -> Quat #expand {
+quat :: (x: float = 0, y: float = 0, z: float = 0, w: float = 1) -> Quat #expand {
    return .{ x = x, y = y, z = z, w = w };
 }
 
+quat :: (xyz: Vec3, w: float) -> Quat #expand {
+   return .{xyz=xyz, w=w};
+}
 
-#scope_file;
+quat_identity :: Quat.{x=0, y=0, z=0, w=1};
+
+operator* :: (a: Quat, b: Quat) -> Quat {
+   r: Quat = ---;
+   r.xyz = a.w*b.xyz + b.w*a.xyz + cross(a.xyz, b.xyz);
+   r.w = a.w*b.w - dot(a.xyz, b.xyz);
+   return r;
+}
+
+operator* :: (q: Quat, v: Vec3) -> Vec3 #symmetric {
+   r: Vec3 = 2*cross(q.xyz, v);
+   return v + q.w*r + cross(q.xyz, r);
+}
+
+operator* :: (q: Quat, r: $T) -> Quat #symmetric {
+   return .{x=q.x*r, y=q.y*r, z=q.z*r, w=q.w*r};
+}
+
+operator+ :: (a: Quat, b: Quat) -> Quat {
+   return .{x=a.x*b.x, y=a.y*b.y, z=a.z*b.z, w=a.w*b.w};
+}
+
+operator/ :: (l: Quat, r: $T) -> Quat {
+   return .{xyz=l.xyz/r, w=l.w/r};
+}
+
+operator== :: (l: Quat, r: Quat) -> bool {
+   return l.x == r.x && l.y == r.y && l.z == r.z && l.w == r.w;
+}
+
+conjugate :: (q: Quat) -> Quat #expand {
+   r: Quat = ---;
+   r.xyz = -q.xyz;
+   r.w = q.w;
+   return r;
+}
+
+// Note(Jesse): iff q is a unit vector
+inverse_unit :: (q: Quat) -> Quat {
+   return inline conjugate(q);
+}
+inverse :: (q: Quat) -> Quat {
+   return conjugate(q)/length(q);
+}
+
+rotation_quat :: (theta: float, axis: Vec3) -> Quat {
+   q: Quat = ---;
+   a := normalize(axis);
+   q.xyz = sin(theta/2.0)*a;
+   q.w = cos(theta/2.0);
+   return q;
+}
+
+operator- :: (q: Quat) -> Quat {
+   r: Quat = ---;
+   r.x = -q.x;
+   r.y = -q.y;
+   r.z = -q.z;
+   r.w = -q.w;
+   return r;
+}
+
+length_squared :: (q: Quat) -> float {
+   return dot(q, q);
+}
+
+length :: (q: Quat) -> float {
+   return math.sqrt(dot(q, q));
+}
+
+dot :: (a: Quat, b: Quat) -> float {
+   return a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w;
+}
+
+slerp :: (a: Quat, _b: Quat, t: float) -> Quat {
+   b := _b;
+   cos_a := dot(a, b);
+   if cos_a < 0.0 {
+      cos_a = -cos_a;
+      b = -b;
+   }
+   alpha := math.acos(cos_a);
+   sin_a := math.sin(alpha);
+
+   p1 := math.sin((1-t)*alpha)/sin_a;
+   p2 := math.sin(t*alpha)/sin_a;
+   return a*p1 + b*p2;
+}
+
+cross :: (a: Vec3, b: Vec3) -> Vec3 {
+   v: Vec3 = ---;
+   v.x = a.y*b.z - a.z*b.y;
+   v.y = a.z*b.x - a.x*b.z;
+   v.z = a.x*b.y - a.y*b.x;
+   return v;
+}
+
+#scope_file
 
 meta :: #import "jc/meta";
 math :: #import "Math"; // @future