/*
   Vec is a generic set of N named values of type T (aka. a mathematical vector)

   The values can be accessed via array index or their common component names:

      - u, v, d
      - x, y, z, w
      - r, g, b, a

      - x,     y,     width, height
      - min_x, min_y, max_x, max_y
      - c0,    c1,    c2,    c3, ... cN

   For most use cases, opt for the named variants of this type:

      Vec2 : Vec(2, float)
      Vec3 : Vec(3, float)
      Vec4 : Vec(4, float)
      Quat : Vec(4, float)

   Sadly, Vec does *NOT* solve the problem of interfacing with
   external vector types (including Jai's 'Math' module).

   To fix this, create two helper macros:

      to_jai :: (v: Vec4) -> jmath.Vector4 #expand {
         return v.(jmath.Vector4,force);
      }
      to_jai :: (v: *Vec4) -> *jmath.Vector4 #expand {
         return v.(*jmath.Vector4);
      }
*/
Vec :: struct(N: int, T: Type)
#modify {
   info := T.(*Type_Info);
   return info.type == .INTEGER || info.type == .FLOAT, "Vec T must be a numeric type (int or float)";
} {
   #assert (N > 0) "Vec N cannot be <= 0";

   // Vecs are backed by an array internally. The #insert block below
   // generates #place'd unions of named fields or cN fields when N is > 4.
   #as components: [N]T;

   #insert -> string {
      b: basic.String_Builder;
      basic.append(*b, "#place components;\n");

      fields :: []string.[
         string.[ "x", "r", "u", "min_x", ""        ],
         string.[ "y", "g", "v", "min_y", ""        ],
         string.[ "z", "b", "d", "max_x", "width"   ],
         string.[ "w", "a", "" , "max_y", "height"  ],
      ];

      for i: 0..N - 1 {
         if i < fields.count {
            basic.append(*b, "union {\n");
            for field: fields[i] if field.count != 0 {
               basic.print_to_builder(*b, "\t%: T = ---;\n", field);
            }
            basic.print_to_builder(*b, "\tc%: T = ---;\n", i);
            basic.append(*b, "};\n");
         }
         else {
            basic.print_to_builder(*b, "c%: T = ---;\n", i);
         }
      }

      // 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);
   };
}

operator [] :: (v: Vec, $$idx: int) -> v.T #no_abc {
   CheckBounds(idx, v.N);
   return v.components[idx];
}

operator *[] :: (v: *Vec, $$idx: int) -> *v.T #no_abc {
   CheckBounds(idx, v.N);
   return *v.components[idx];
}

operator []= :: (v: *Vec, $$idx: int, value: v.T) #no_abc {
   CheckBounds(idx, v.N);
   v.components[idx] = value;
}

for_expansion :: (v: *Vec, body: Code, flags: For_Flags) #expand {
   for i: 0..v.N - 1 {
      `it_index := i;
      #if flags & .POINTER {
         `it := *v.components[i];
      } else {
         `it := v.components[i];
      }

      #insert,scope(body) body;
   }
}

operator + :: inline (l: Vec, r: Vec(l.N, l.T)) -> Vec(l.N, l.T) #no_abc {
   res: Vec(l.N, l.T) = ---;
   #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) = ---;
   #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) = ---;
   #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) = ---;
   #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 type_is_scalar(R), "type is not integer or float"; } {
   res: Vec(l.N, l.T) = ---;
   #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 type_is_scalar(R), "type is not integer or float"; } {
   res: Vec(l.N, l.T) = ---;
   #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 type_is_scalar(R), "type is not integer or float"; } {
   res: Vec(l.N, l.T) = ---;
   #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 type_is_scalar(R), "type is not integer or float"; } {
   res: Vec(l.N, l.T) = ---;
   #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 type_is_scalar(R), "type is not integer or float"; } {
   res: Vec(l.N, l.T) = ---;
   #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 {
   #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) = ---;
   #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) = ---;
   #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;
}

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 {
      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;
}

// @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 {
      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;
}

clamp :: (v: Vec, low: v.T, high: v.T) -> Vec(v.N, v.T) #no_abc {
   r: Vec(v.N, v.T) = ---;
   n := v.N - 1;
   while n >= 0 {
      if v[n] < low
         r[n] = low;
      else if v[n] > high
         r[n] = high;
      else
         r[n] = v[n];
      n -= 1;
   }
   return r;
}

dot :: (a: Vec, b: Vec(a.N, a.T)) -> a.T #no_abc {
   sum: a.T;
   n := a.N - 1;
   while n >= 0 {
      sum += a[n]*b[n];
      n -= 1;
   }
   return sum;
}

length_squared :: (v: Vec) -> float #no_abc {
   return inline dot(v, v);
}

length :: (v: Vec) -> float #no_abc {
   return math.sqrt(inline dot(v, v));
}

abs :: (v: Vec) -> Vec(v.N, v.T) #no_abc {
   r: Vec(v.N, v.T) = ---;
   #if v.N <= 4 {
      r = v;
      if r.x < 0 r.x = -r.x;
      if r.y < 0 r.y = -r.y;
      #if v.N >= 3 {
         if r.z < 0 r.z = -r.z;
      }
      #if v.N == 4 {
         if r.w < 0 r.w = -r.w;
      }
   }
   else {
      n := v.N - 1;
      while n >= 0 {
         if v[n] < 0
            r[n] = -v[n];
         else
            r[n] = v[n];
         n -= 1;
      }
   }
   return r;
}

// @todo(Jesse): SIMD
norm :: normalize;
normalize :: (v: Vec) -> Vec(v.N, v.T) #no_abc {
   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 {
   r: Vec(a.N, a.T) = ---;
   n := a.N - 1;
   while n >= 0 {
      r[n] = a[n] + t*(b[n] - a[n]);
      n -= 1;
   }
   return r;
}

// Note(Jesse): I don't think this is needed for bigger vectors
reflect :: (v: Vec3, p: Vec3) -> Vec3 #no_abc {
   projection := p*dot(v, p)/length_squared(p);
   return 2*projection - v;
}
reflect :: (v: Vec2, p: Vec2) -> Vec2 #no_abc {
   projection := p*dot(v, p)/length_squared(p);
   return 2*projection - v;
}

round :: (v: Vec($N, $T)) -> Vec(N, T) #no_abc
#modify { return type_is_float(T), "Used non-float vector on round"; } {
   r: Vec(N, T) = ---;
   n := N - 1;
   while n >= 0 {
      if v[n] < 0
         r[n] = (v[n] - 0.5).(int).(float);
      else
         r[n] = (v[n] + 0.5).(int).(float);
      n -= 1;
   }
   return r;
}

// Concrete vector types (the usual cases)

Vec2 :: Vec(2, float);
Vec3 :: Vec(3, float);
Vec4 :: Vec(4, float);
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 type_is_float(T), "use v2i for integer arguments"; }
#expand {
   return .{ x = x, y = y };
}

v2i :: (x: $T = 0, y: T = 0) -> Vec(2, T)
#modify { return type_is_integer(T), "use v2f for float arguments"; }
#expand {
   return .{ x = x, y = y };
}

v3f :: (x: $T = 0, y: T = 0, z: T = 0) -> Vec3
#modify { return type_is_float(T), "use v3i for integer arguments"; }
#expand {
   return .{ x = x, y = y, z = z };
}

v3i :: (x: $T = 0, y: T = 0, z: T = 0) -> Vec(3, T)
#modify { return type_is_integer(T), "use v3f for float arguments"; }
#expand {
   return .{ x = x, y = y, z = z };
}

v4f :: (x: $T = 0, y: T = 0, z: T = 0, w: T = 0) -> Vec4
#modify { return type_is_float(T), "use v4i for integer arguments"; }
#expand {
   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 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 = 1) -> Quat #expand {
   return .{ x = x, y = y, z = z, w = w };
}

quat :: (xyz: Vec3, w: float) -> Quat #expand {
   return .{xyz=xyz, w=w};
}

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;

#if RUN_TESTS #run,stallable {
   Test(basic.tprint("%: Vec2", UNITS), t => {
      {
         a: Vec2 = v2f(0.0, 1.0);
         b: Vec2 = v2f(1.0, 2.0);

         Expect(a + b == v2f(1.0, 3.0));
         Expect(b - a == v2f(1.0, 1.0));
         Expect(a*b == v2f(0.0, 2.0));
         Expect(a/b == v2f(0.0, 0.5));
      }

      {
         a: Vec(2, int) = v2i(2, 1);
         b: Vec(2, int) = v2i(1, 0);
         Expect(a + b == v2i(3, 1));
         Expect(b - a == v2i(-1, -1));
         Expect(a*b == v2i(2, 0));
         Expect(b/a == v2i(0, 0));
      }
      {
         a: Vec2 = v2f(2.3, -4.1);
         b: Vec2 = v2f(1.0, 3.6);
         c := min(a, b);
         Expect(c == v2f(1.0, -4.1));
         c = max(a, b);
         Expect(c == v2f(2.3, 3.6));
      }
   });

   Test(basic.tprint("%: Vec3", UNITS), t => {
      {
         a: Vec3 = v3f(0.0, 1.0, 2.0);
         b: Vec3 = v3f(1.0, 2.0, 3.0);
         Expect(a + b == v3f(1.0, 3.0, 5.0));
         Expect(b - a == v3f(1.0, 1.0, 1.0));
         Expect(a*b == v3f(0.0, 2.0, 6.0));
         Expect(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);
         Expect(reflect(a, b) == v3f(1.0, -1.0, 0.0));
         Expect(round(v3f(1.2, 1.7, 1.5)) == v3f(1.0, 2.0, 2.0));
         Expect(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);
         Expect(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);
         Expect(c == v3f(1.0, -3.6, 5.0));
         c = max(a, b);
         Expect(c == v3f(2.3, 4.1, 9.0));
      }
   });

   Test(basic.tprint("%: Vec4", UNITS), t => {
      a: Vec4 = v4f(2.25, 1.0, 2.0, 1.0);
      b: Vec4 = v4f(4.0, 2.0, 3.0, 1.0);
      Expect(a + b == v4f(6.25, 3.0, 5.0, 2.0));
      Expect(b - a == v4f(1.75, 1.0, 1.0, 0.0));
      Expect(a*b == v4f(9.0, 2.0, 6.0, 1.0));
      Expect(a/b == v4f(0.5625, 0.5, 2.0/3.0, 1.0));
   });

   Test(basic.tprint("%: VecN", UNITS), t => {
      a: Vec(16, float);
      b: Vec(16, float);
      for *a {
         it.* = xx it_index;
      }
      for *b {
         it.* = xx(it_index + 1);
      }
      Expect(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]});
      Expect(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]});
      Expect(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]});

      Expect(min(a, b) == a);
      Expect(max(a, b) == b);
      Expect(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]});
      Expect(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]});
      Expect(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]};
      Expect(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]});
      Expect(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]});

      Expect(dot(a, b) == 1360.0);
      Expect(abs(a) == a);

      c := a;
      for *c { // Check making every other component negative
         if it_index%2 == 0 then it.* = -it.*;
      }

      Expect(abs(c) == a);
      Expect(float_eq(length(normalize(a)), 1));
      Expect(a + 2 == Vec(16, float).{.[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]});
      Expect(a - 2 == Vec(16, float).{.[-2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]});
      Expect(a*2 == Vec(16, float).{.[0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30]});
   });

   Test(basic.tprint("%: Quat", UNITS), t => {
      qi := quat_identity;
      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));
      qc := conjugate(q);
      inv_q := inverse(q);
      Expect(q*inv_q == qi);

      q1 := quat(2, 0, 0, 0);
      q2  = quat(1, 1, -1, 0);
      c := q1*q2*conjugate(q1);
      Expect(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);
      c1 := q*quat(p, 0)*conjugate(q);
      c2 := q*p;
      Expect(v3_eq(c2, v3f(math.sqrt(2.0), math.sqrt(2.0), 0.0)));
      Expect(v3_eq(c1.xyz, c2));

      q = rotation_quat(from_rad(PI), v3f(0.0, 1.0, 0.0));
      m := rotation_mat4(q);
      p1 := v4f(2.0, 0.0, 0.0, 1.0);
      Expect(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);
      q = rotation_quat(from_rad(PI/4.0), v3f(0.0, 0.0, 1.0));
   });
}