Added collision detection for several shapes and SIMD matrix multiplies

math/common.jai

@@ -92,6 +92,22 @@ abs :: (v: $T) -> T
    return ifx v < 0 then -v else v;
 }
 
+sign :: (x: int) -> int {
+   if x < 0
+      return -1;
+   else if x > 0
+      return 1;
+   return 0;
+}
+
+sign :: (x: float) -> float {
+   if x < 0.0
+      return -1.0;
+   else if x > 0.0
+      return 1.0;
+   return 0.0;
+}
+
 #scope_file
 
 // sin, cos

math/mat.jai

@@ -1,3 +1,4 @@
+
 Mat2 :: #type,distinct Vec(2 * 2, float);
 Mat4 :: #type,distinct Vec(4 * 4, float);
 
@@ -155,27 +156,7 @@ make_look_at :: (camera: Vec3, at: Vec3, up_vector: Vec3) -> Mat4
 }
 
 operator* :: (a: Mat4, b: Mat4) -> Mat4 {
-   r: Mat4 = ---;
+   return mul_sse(a, b);
-   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
@@ -402,8 +383,196 @@ inverse :: (m: Mat4) -> Mat4 {
 }
 
 
+// Note(Jesse): While this is surprisingly faster in debug, in release the naive
+// implementation and the sse version are disappointingly very similar in cycle count.
+// This is even more confusing as the naive when optimized still uses single single arithmatic,
+// while we're using packed single arithmatic
+mul_sse :: inline (a: Mat4, b: Mat4) -> Mat4 {
+   r: Mat4 = ---;
+
+   #asm {
+      r0: vec;
+      r1: vec;
+      r2: vec;
+      r3: vec;
+
+      // Loading
+      movups a_0:, [*a     ];
+      movups a_1:, [*a + 16];
+      movups a_2:, [*a + 32];
+      movups a_3:, [*a + 48];
+      movups b_0:, [*b     ];
+      movups b_1:, [*b + 16];
+      movups b_2:, [*b + 32];
+      movups b_3:, [*b + 48];
+
+      // row0
+      // We're using pshufd because shufps references 'a' in the lower, and 'b' in the upper.
+      // pshufd only reference a_0. So shufps gets from both the dest and the src registers
+      pshufd r0, a_0, 0x00;
+      pshufd r1, a_0, 0x55;
+      pshufd r2, a_0, 0xaa;
+      pshufd r3, a_0, 0xff;
+      mulps r0, b_0;
+      mulps r1, b_1;
+      mulps r2, b_2;
+      mulps r3, b_3;
+      addps r0, r1;
+      addps r2, r3;
+      addps r0, r2;
+      movups [*r], r0;
+
+      // row1
+      pshufd r0, a_1, 0x00;
+      pshufd r1, a_1, 0x55;
+      pshufd r2, a_1, 0xaa;
+      pshufd r3, a_1, 0xff;
+      mulps r0, b_0;
+      mulps r1, b_1;
+      mulps r2, b_2;
+      mulps r3, b_3;
+      addps r0, r1;
+      addps r2, r3;
+      addps r0, r2;
+      movups [*r + 16], r0;
+
+      // row2
+      pshufd r0, a_2, 0x00;
+      pshufd r1, a_2, 0x55;
+      pshufd r2, a_2, 0xaa;
+      pshufd r3, a_2, 0xff;
+      mulps r0, b_0;
+      mulps r1, b_1;
+      mulps r2, b_2;
+      mulps r3, b_3;
+      addps r0, r1;
+      addps r2, r3;
+      addps r0, r2;
+      movups [*r + 32], r0;
+
+      // row3
+      pshufd r0, a_3, 0x00;
+      pshufd r1, a_3, 0x55;
+      pshufd r2, a_3, 0xaa;
+      pshufd r3, a_3, 0xff;
+      mulps r0, b_0;
+      mulps r1, b_1;
+      mulps r2, b_2;
+      mulps r3, b_3;
+      addps r0, r1;
+      addps r2, r3;
+      addps r0, r2;
+      movups [*r + 48], r0;
+   }
+
+   // 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): This procedure will crash if it tries to store or access an element across a cache line
+mul_sse_aligned :: inline (dst: *Mat4, a: *Mat4, b: *Mat4) {
+   // Todo(Jesse): Find out what we want to do for debug checks
+   // assert(is_aligned(dst, 64) && is_aligned(a, 64) && is_aligned(b, 64), "If you use this procedure, dst, a, and b's base pointer must be aligned to a 64 byte boundary");
+
+   #asm {
+      r0: vec;
+      r1: vec;
+      r2: vec;
+      r3: vec;
+
+      // Loading
+      movaps a_0:, [a     ];
+      movaps a_1:, [a + 16];
+      movaps a_2:, [a + 32];
+      movaps a_3:, [a + 48];
+      movaps b_0:, [b     ];
+      movaps b_1:, [b + 16];
+      movaps b_2:, [b + 32];
+      movaps b_3:, [b + 48];
+
+      // row0
+      // We're using pshufd because shufps references 'a' in the lower, and 'b' in the upper.
+      // pshufd only reference a_0. So shufps gets from both the dest and the src registers
+      pshufd r0, a_0, 0x00;
+      pshufd r1, a_0, 0x55;
+      pshufd r2, a_0, 0xaa;
+      pshufd r3, a_0, 0xff;
+      mulps r0, b_0;
+      mulps r1, b_1;
+      mulps r2, b_2;
+      mulps r3, b_3;
+      addps r0, r1;
+      addps r2, r3;
+      addps r0, r2;
+      movaps [dst], r0;
+
+      // row1
+      pshufd r0, a_1, 0x00;
+      pshufd r1, a_1, 0x55;
+      pshufd r2, a_1, 0xaa;
+      pshufd r3, a_1, 0xff;
+      mulps r0, b_0;
+      mulps r1, b_1;
+      mulps r2, b_2;
+      mulps r3, b_3;
+      addps r0, r1;
+      addps r0, r2;
+      addps r0, r3;
+      movaps [dst + 16], r0;
+
+      // row2
+      pshufd r0, a_2, 0x00;
+      pshufd r1, a_2, 0x55;
+      pshufd r2, a_2, 0xaa;
+      pshufd r3, a_2, 0xff;
+      mulps r0, b_0;
+      mulps r1, b_1;
+      mulps r2, b_2;
+      mulps r3, b_3;
+      addps r0, r1;
+      addps r0, r2;
+      addps r0, r3;
+      movaps [dst + 32], r0;
+
+      // row3
+      pshufd r0, a_3, 0x00;
+      pshufd r1, a_3, 0x55;
+      pshufd r2, a_3, 0xaa;
+      pshufd r3, a_3, 0xff;
+      mulps r0, b_0;
+      mulps r1, b_1;
+      mulps r2, b_2;
+      mulps r3, b_3;
+      addps r0, r1;
+      addps r0, r2;
+      addps r0, r3;
+      movaps [dst + 48], r0;
+   }
+}
+
 #scope_file;
 
+is_aligned :: (p: *void, bound: int) -> bool {
+   return cast(int)p & (bound - 1) == 0;
+}
+
 #if #exists(RUN_TESTS) #run {
    test :: #import "jc/meta/test";
 

math/shape.jai

@@ -1,8 +1,7 @@
 Rect :: #type,distinct Vec(4, RECT_TYPE);
 
 Circle :: struct {
-   x,y: float;
+   x,y,r: float;
-   r: float;
    #place x;
    pos: Vec2;
 }
@@ -13,8 +12,14 @@ Line :: struct {
    a: Vec2;
    #place x1;
    b: Vec2;
+   #place x0;
+   e: [2]Vec2;
 }
 
+Triangle :: [3]Vec2;
+
+ORIGIN :: Vec2.{.[0, 0]};
+
 make_rect :: (x: RECT_TYPE, y: RECT_TYPE, w: RECT_TYPE, h: RECT_TYPE) -> Rect {
    r: Rect = ---;
    r.x = x;
@@ -78,6 +83,20 @@ inside :: (c: Circle, p: Vec2) -> bool #symmetric {
    return dp.x*dp.x + dp.y*dp.y <= c.r*c.r;
 }
 
+inside :: (t: Triangle, p: Vec2) -> bool #symmetric {
+   cay := t[2].y - t[0].y;
+   pay := p.y - t[0].y;
+   cax := t[2].x - t[0].x;
+   bay := t[1].y - t[0].y;
+   bax := t[1].x - t[0].x;
+
+   denom := bay*cax - bax*cay;
+   w1 := t[0].x*cay + pay*cax - p.x*cay;
+   w1 /= denom;
+   w2 := (p.y - t[0].y - w1*bay)/cay;
+   return w1 >= 0 && w2 >= 0 && w1 + w2 <= 1;
+}
+
 collides :: (c1: Circle, c2: Circle) -> bool {
    dp := c2.pos - c1.pos;
    return c1.r + c2.r >= length(dp);
@@ -85,18 +104,18 @@ collides :: (c1: Circle, c2: Circle) -> bool {
 
 // Note(Jesse): This is using sdfs. Very elegant
 collides :: (c: Circle, r: Rect) -> bool #symmetric {
-   // We need to 'transpose' all the math to the origin (centered at 0,0)
+   // We need to 'transpose' all the math to the origin the center of the circle acting as the origin (0, 0)
    r_center: Vec2 = get_center(r);
    p := r_center - c.pos;
-   d: Vec2 = abs(p) - v2f(r.width/2, r.height/2);
+   d: Vec2 = abs(p) - Vec2.{.[r.width/2, r.height/2]};
    dist := length(max(d, 0.0)) + min(max(d.x, d.y), 0.0);
    dist -= c.r; // 'adding' the distance of the circle
    return dist <= 0.0;
 }
 
-// Note(Jesse): Minkowski difference
+// Note(Jesse): Minkowski difference but instead of the origin we check if the other rectangle center
 collides :: (r1: Rect, r2: Rect) -> bool {
-   r: Rect = make_rect(r1.x - r2.width/2, r1.y - r2.width/2, r1.width + r2.width, r1.height + r2.height);
+   r: Rect = Rect.{.[r1.x - r2.width/2, r1.y - r2.width/2, r1.width + r2.width, r1.height + r2.height]};
    return inside(r, get_center(r2));
 }
 
@@ -124,22 +143,224 @@ collides :: (s1: Line, s2: Line) -> bool {
    return (a >= 0 && a <= 1) && (b >= 0 && b <= 1);
 }
 
-// Todo(Jesse): I want to find a better line segment rectangle algorithm.
+// Cohen-sutherland's line clipping algorithm
-//    one option is a line clipping algorithm. Cohen-Sutherland algorithm
+collides :: (r: Rect, l: Line) -> bool #symmetric {
-//    would be the best option. It does accepts/rejects fastest
+   INSIDE :: 0b0000;
-//    but still has a while in it, and a few ifs.
+   LEFT   :: 0b0001;
-collides :: (r: Rect, seg: Line) -> bool #symmetric {
+   RIGHT  :: 0b0010;
-   p1 := Vec2.{r.x,           r.y};
+   BOTTOM :: 0b0100;
-   p2 := Vec2.{r.x + r.width, r.y};
+   TOP    :: 0b1000;
-   p3 := Vec2.{r.x,           r.y + r.height};
-   p4 := Vec2.{r.x + r.width, r.y + r.height};
 
-   l1 := Line.{a=p1, b=p2};
+   compute_out_code :: (x: float, y: float) -> int #expand {
-   l2 := Line.{a=p1, b=p3};
+      code := INSIDE;
-   l3 := Line.{a=p3, b=p4};
+      if x < xmin
-   l4 := Line.{a=p2, b=p4};
+         code |= LEFT;
+      else if x > xmax
+         code |= RIGHT;
+      if y < ymin
+         code |= BOTTOM;
+      else if y > ymax
+         code |= TOP;
+      return code;
+   }
 
-   return collides(seg, l1) || collides(seg, l2) || collides(seg, l3) || collides(seg, l4);
+   xmin := r.x;
+   xmax := r.x + r.width;
+   ymin := r.y;
+   ymax := r.y + r.height;
+
+   x0 := l.a.x;
+   y0 := l.a.y;
+   x1 := l.b.x;
+   y1 := l.b.y;
+
+   oc0 := compute_out_code(l.a.x, l.a.y);
+   oc1 := compute_out_code(l.b.x, l.b.y);
+   accept: bool;
+   while true {
+      if !(oc0 | oc1) {
+         accept = true;
+         break;
+      }
+      else if oc0 & oc1 {
+         break;
+      }
+      else {
+         x, y: float = ---;
+         oc_out := ifx oc1 > oc0 then oc1 else oc0;
+
+         // find intersection point;
+         // slope = (y1 - y0)/(x1 - x0)
+         // x = x0 + (1/slope)*(ym - y0), ym is ymin or ymax
+         // y = y0 + slope*(xm - x0), xm is xmin or xmax
+         // outcode bit guarantees the denom is non-zero
+         if oc_out & TOP {
+            x = x0 + (x1 - x0)*(ymax - y0)/(y1 - y0);
+            y = ymax;
+         }
+         else if oc_out & BOTTOM {
+            x = x0 + (x1 - x0)*(ymin - y0)/(y1 - y0);
+            y = ymin;
+         }
+         else if oc_out & RIGHT {
+            y = y0 + (y1 - y0)*(xmax - x0)/(x1 - x0);
+            x = xmax;
+         }
+         else {
+            y = y0 + (y1 - y0)*(xmin - x0)/(x1 - x0);
+            x = xmin;
+         }
+
+         // Move outside points to intersection point
+         if oc_out == oc0 {
+            x0 = x;
+            y0 = y;
+            oc0 = compute_out_code(x0, y0);
+         }
+         else {
+            x1 = x;
+            y1 = y;
+            oc1 = compute_out_code(x1, y1);
+         }
+      }
+   }
+   return accept;
+}
+
+// Triangle SDF with the circle position inplace of the origin and a circle expansion of the sdf
+collides :: (t: Triangle, c: Circle) -> bool #symmetric {
+   e0 := t[1] - t[0];
+   e1 := t[2] - t[1];
+   e2 := t[0] - t[2];
+   v0 := c.pos - t[0];
+   v1 := c.pos - t[1];
+   v2 := c.pos - t[2];
+   pq0 := v0 - e0*min(max(dot(v0,e0)/dot(e0,e0), 0.0), 1.0);
+   pq1 := v1 - e1*min(max(dot(v1,e1)/dot(e1,e1), 0.0), 1.0);
+   pq2 := v2 - e2*min(max(dot(v2,e2)/dot(e2,e2), 0.0), 1.0);
+   s := sign(e0.x*e2.y - e0.y*e2.x);
+   d := min(min(Vec2.{.[dot(pq0,pq0), s*(v0.x*e0.y - v0.y*e0.x)]},
+           Vec2.{.[dot(pq1,pq1), s*(v1.x*e1.y - v1.y*e1.x)]}),
+           Vec2.{.[dot(pq2,pq2), s*(v2.x*e2.y - v2.y*e2.x)]});
+   dist := -math.sqrt(d.x)*sign(d.y);
+   return dist < c.r;
+}
+
+collides :: (t: Triangle, r: Rect) -> bool #symmetric {
+   p1 := Vec2.{.[r.x,           r.y]};
+   p2 := Vec2.{.[r.x + r.width, r.y]};
+   p3 := Vec2.{.[r.x,           r.y + r.height]};
+   p4 := Vec2.{.[r.x + r.width, r.y + r.height]};
+
+   rt1: Triangle = .[p1, p4, p2];
+   rt2: Triangle = .[p1, p3, p4];
+
+   if collides(t, rt1) || collides(t, rt2)
+      return true;
+   return false;
+}
+
+collides :: (t1: Triangle, t2: Triangle) -> bool {
+   return inline gjk(t1, t2);
+}
+
+collides :: (t: Triangle, l: Line) -> bool #symmetric {
+   return inline gjk(t, l.e);
+}
+
+// Note(Jesse): 2D
+gjk :: (s1: []Vec2, s2: []Vec2) -> bool {
+   S: Simplex2D;
+   S.a = gjk_support(s1, s2, s1[0] - s2[0]);
+   d := -S.a; // This makes a vector from S.a towards the origin (ORIGIN - S.a)
+   S.points = 1;
+   while true {
+      A := gjk_support(s1, s2, d);
+      if !same_dir(A, d) // If A is even towards the origin. Otherwise there's no way they intersect
+         return false;
+      S.points += 1;
+      S.c = S.b; // we want S.a to be the newest point. S.c is the oldest
+      S.b = S.a;
+      S.a = A;
+      if gjk_do_simplex(*S, *d)
+         return true;
+   }
+   return false;
+}
+
+#scope_file
+
+basic :: #import "Basic";
+
+// Used with gjk, not for general use
+Simplex2D :: struct {
+   a,b,c: Vec2;
+   points: int;
+}
+
+// Used with gjk. Simplified for 2D
+triple_cross :: inline (a: Vec2, b: Vec2, c: Vec2) -> Vec2 {
+   // return cross(cross(v3(a, 0), v3(b, 0)), v3(c, 0)).xy;
+   return Vec2.{.[-a.x*b.y*c.y + a.y*b.x*c.y, a.x*b.y*c.x - a.y*b.x*c.x]};
+}
+
+// Sets the direction for the next support function based on the simplex we have and how it's orientated
+// to the origin: 0, 0. Because we know which point was added last (S.a) we can deduce directions we don't
+// have to check
+gjk_do_simplex :: (S: *Simplex2D, dir: *Vec2) -> bool {
+   if S.points == {
+      case 2;
+         AB := S.b - S.a;
+         AO := -S.a;
+         if same_dir(AB, AO) {
+            dir.* = triple_cross(AB, AO, AB);
+         }
+         else {
+            dir.* = AO;
+            S.points = 1;
+         }
+      case 3;
+         AO := -S.a;
+         // Because this is 2D, we don't need a large 3 point simplex case, or a 4 point case
+         AB := S.b - S.a;
+         AC := S.c - S.a;
+         abperp := triple_cross(AC, AB, AB);
+         acperp := triple_cross(AB, AC, AC);
+         if same_dir(abperp, AO) {
+            S.points = 2;
+            dir.* = abperp;
+         }
+         else if same_dir(acperp, AO) {
+            S.b = S.c;
+            S.points = 2;
+            dir.* = acperp;
+         }
+         else {
+            return true; // origin must be inside the triangle
+         }
+  }
+  return false;
+}
+
+// Returns the vertex furthest along the dir vector for two polygons. We could split it up
+// to take something that can't be represented by a list of points, like a circle.
+gjk_support :: inline (a: []Vec2, b: []Vec2, dir: Vec2) -> Vec2 {
+   helper :: (t: []Vec2, dir: Vec2) -> Vec2 {
+      p := t[0];
+      best := dot(t[0], dir);
+      for 1..t.count - 1 {
+         v := t[it];
+         if dot(t[it], dir) > best {
+            best = dot(v, dir);
+            p = v;
+         }
+      }
+      return p;
+   }
+
+   v := helper(a,  dir);
+   w := helper(b, -dir);
+   return v - w;
 }
 
 ///////
@@ -149,9 +370,48 @@ get_center :: (r: Rect) -> Vec2 {
    return v2f(r.x + r.width/2, r.y + r.height/2);
 }
 
-#scope_file
-
-basic :: #import "Basic";
-
 #if #exists(RUN_TESTS) #run {
+   // Just a formality to get proper compile errors
+
+   b: bool;
+
+   r1 := Rect.{.[100, 150, 250, 250]};
+   c1 := Circle.{150, 500, 50};
+   l1 := Line.{400, 100, 650, 150};
+   t1: Triangle;
+   t1[0] = v2f(600.0, 500.0);
+   t1[1] = v2f(650.0, 250.0);
+   t1[2] = v2f(625.0, 200.0);
+
+   p: Vec2 = v2f(100.0, 100.0);
+   r2 := Rect.{.[100, 100, 10, 10]};
+   c2 := Circle.{100, 100, 10};
+   l2 := Line.{100, 100, 200, 200};
+   t2: Triangle;
+   t2[0] = v2f(50.0,      50.0 + 20.0);
+   t2[1] = v2f(50.0 + 20.0, 50.0 - 20.0);
+   t2[2] = v2f(50.0 - 20.0, 50.0 - 20.0);
+
+   b |= inside(r1, p);
+   b |= collides(r1, r2);
+   b |= collides(r1, l2);
+   b |= collides(r1, c2);
+   b |= collides(r1, t2);
+
+   b |= inside(c1, p);
+   b |= collides(c1, r2);
+   b |= collides(c1, l2);
+   b |= collides(c1, c2);
+   b |= collides(c1, t2);
+
+   b |= collides(l1, r2);
+   b |= collides(l1, l2);
+   b |= collides(l1, c2);
+   b |= collides(l1, t2);
+
+   b |= inside(t1, p);
+   b |= collides(t1, r2);
+   b |= collides(t1, l2);
+   b |= collides(t1, c2);
+   b |= collides(t1, t2);
 }

math/vec.jai

@@ -371,7 +371,7 @@ max :: (l: Vec, r: l.T) -> Vec(l.N, l.T) #no_abc {
    return res;
 }
 
-// @todo(jesse): SIMD
+// @todo(jesse): SIMD?
 ceil :: (l: Vec) -> Vec(l.N, l.T) #no_abc {
    r: Vec(l.N, l.T) = ---;
    n := l.N - 1;
@@ -385,7 +385,7 @@ ceil :: (l: Vec) -> Vec(l.N, l.T) #no_abc {
    return r;
 }
 
-// @todo(jesse): SIMD
+// @todo(jesse): SIMD?
 floor :: (l: Vec) -> Vec(l.N, l.T) #no_abc {
    r: Vec(l.N, l.T) = ---;
    n := l.N - 1;
@@ -651,6 +651,13 @@ cross :: (a: Vec3, b: Vec3) -> Vec3 {
    return v;
 }
 
+same_dir :: (a: Vec2, b: Vec2) -> bool {
+   return dot(a, b) > 0;
+}
+same_dir :: (a: Vec3, b: Vec3) -> bool {
+   return dot(a, b) > 0;
+}
+
 #scope_file;
 
 #if RUN_TESTS #run,stallable {