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 = .[
         r0.x, r0.y, r0.z, r0.w,
         r1.x, r1.y, r1.z, r1.w,
         r2.x, r2.y, r2.z, r2.w,
         r3.x, r3.y, r3.z, r3.w,
      ],
   };
}

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._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).*;
}

operator *[] :: inline (m: *Mat4, $$idx: int) -> *Vec4 #no_abc {
   N :: Mat4.{}.N; // @note(judah): dumb workaround for not being able to access pointer type field constants
   meta.check_bounds(idx * 4, N);
   return (*m.components[idx * 4]).(*Vec4);
}

operator []= :: inline (m: *Mat4, $$idx: int, value: Vec4) #no_abc {
   N :: Mat4.{}.N;
   meta.check_bounds(idx * 4, N);

   ptr := (*m.components[idx * 4]).(*Vec4);
   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);
}


#scope_file;

#if #exists(RUN_TESTS) #run {
   test :: #import "jc/meta/test";

   test.run(basic.tprint("%: Mat2", UNITS), t => {
      {
         identity := m2_identity;
         a: Mat2 = Mat2.{.[1, 2, 3, 4]};
         test.expect(t, a*identity == a);
      }
      {
         rotator := rotation_mat2(from_rad(PI));
         v2 := v2f(1.0, 0.5);
         v2 = rotator*v2;
         test.expect(t, v2_eq(v2, v2f(-1.0, -0.5)));
         v2 = rotator*v2;
         test.expect(t, v2_eq(v2, v2f(1.0, 0.5)));
      }
      {
         m := m2(1.0, 3.0, 2.0, 2.0);
         det := determinate(m);
         test.expect(t, det == -4);
      }
      {
         rotator := rotation_mat2(from_rad(PI));
         inv_rotator := inverse(rotator);
         v2 := v2f(0.1, 1.0);

         test.expect(t, inv_rotator*rotator == m2_identity);
         test.expect(t, inv_rotator*(rotator*v2) == v2);
         v2 = rotator*v2;
         v2 = inv_rotator*v2;
      }
   });

   test.run(basic.tprint("%: Mat4", UNITS), 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);
      }
      {
         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));
      }
      {
         rotator := rotation_mat4(v3f(0.0, 1.0, 0.0), from_rad(PI));
         v4 := v4f(1.0, 0.5, 0.1, 1.0);
         v4 = rotator*v4;
         test.expect(t, v4_eq(v4, v4f(-1.0, 0.5, -0.1, 1.0)));
         v4 = rotator*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;
         det := determinate(camera_rotator);
      }
      {
         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);

         test.expect(t, inv_rotator*rotator == m4_identity);

         v4 = rotator*v4;
         v4 = inv_rotator*v4;
      }
   });
}