Add matrix inverse calculation

src/vec3.zig

@@ -32,6 +32,10 @@ pub fn vec3Dot(a: Vec3F64, b: Vec3F64) f64 {
     return @reduce(.Add, a.data * b.data);
 }
 
+pub fn vec3MulScalar(a: Vec3F64, b: f64) Vec3F64 {
+    return .init(a.x() * b, a.y() * b, a.z() * b);
+}
+
 pub fn vec3Cross(a: Vec3F64, b: Vec3F64) Vec3F64 {
     const x_component = a.y() * b.z() - a.z() * b.y();
     const y_component = a.z() * b.x() - a.x() * b.z();
@@ -83,6 +87,45 @@ pub fn mat3Mul(a: Mat3F64, b: Mat3F64) Mat3F64 {
     return .{ .row1 = c_row1, .row2 = c_row2, .row3 = c_row3 };
 }
 
+pub fn mat3Det(a: Mat3F64) f64 {
+    return a.row1.x() * (a.row2.y() * a.row3.z() - a.row2.z() * a.row3.y()) -
+        a.row1.y() * (a.row2.x() * a.row3.z() - a.row2.z() * a.row3.x()) +
+        a.row1.y() * (a.row2.x() * a.row3.y() - a.row2.y() * a.row3.x());
+}
+
+pub fn mat3Inv(a: Mat3F64) Mat3F64 {
+    const inv_det = 1.0 / mat3Det(a);
+    const row1 = vec3MulScalar(
+        .init(
+            a.row2.y() * a.row3.z() - a.row2.z() * a.row3.y(),
+            a.row1.z() * a.row3.y() - a.row1.y() * a.row3.z(),
+            a.row1.y() * a.row2.z() - a.row1.z() * a.row2.y(),
+        ),
+        inv_det,
+    );
+    const row2 = vec3MulScalar(
+        .init(
+            a.row2.z() * a.row3.x() + a.row2.x() * a.row3.z(),
+            a.row1.x() * a.row3.z() - a.row1.z() * a.row3.x(),
+            a.row1.z() * a.row2.x() - a.row1.x() * a.row2.z(),
+        ),
+        inv_det,
+    );
+    const row3 = vec3MulScalar(
+        .init(
+            a.row2.x() * a.row3.y() + a.row2.y() * a.row3.x(),
+            a.row1.y() * a.row3.x() - a.row1.x() * a.row3.y(),
+            a.row1.x() * a.row2.y() - a.row1.y() * a.row2.x(),
+        ),
+        inv_det,
+    );
+    return .{
+        .row1 = row1,
+        .row2 = row2,
+        .row3 = row3,
+    };
+}
+
 pub fn vec3MulMat3(a: Mat3F64, b: Vec3F64) Vec3F64 {
     return .init(
         vec3Dot(a.row1, b),
@@ -267,6 +310,36 @@ test "Mat3 mul works" {
     try std.testing.expect(c.row3.z() == 9);
 }
 
+test "mat determinant works" {
+    const a: Mat3F64 = .{
+        .row1 = Vec3F64.init(1, 2, 3),
+        .row2 = Vec3F64.init(0, 1, 5),
+        .row3 = Vec3F64.init(0, 0, 9),
+    };
+    try std.testing.expectApproxEqAbs(9, mat3Det(a), 1e-12);
+}
+
+test "mat inverse works" {
+    const a: Mat3F64 = .{
+        .row1 = Vec3F64.init(1, 2, 3),
+        .row2 = Vec3F64.init(0, 1, 5),
+        .row3 = Vec3F64.init(0, 0, 9),
+    };
+
+    const ainv = mat3Inv(a);
+    std.debug.print("{}", .{ainv});
+    std.debug.print("{}", .{ainv.row2.z()});
+
+    try std.testing.expectApproxEqAbs(1, ainv.row1.x(), 1e-12);
+    try std.testing.expectApproxEqAbs(-2, ainv.row1.y(), 1e-12);
+    try std.testing.expectApproxEqAbs(7.0 / 9.0, ainv.row1.z(), 1e-12);
+    try std.testing.expectApproxEqAbs(0, ainv.row2.x(), 1e-12);
+    try std.testing.expectApproxEqAbs(1, ainv.row2.y(), 1e-12);
+    try std.testing.expectApproxEqAbs(-5.0 / 9.0, ainv.row2.z(), 1e-12);
+    try std.testing.expectApproxEqAbs(0, ainv.row3.x(), 1e-12);
+    try std.testing.expectApproxEqAbs(0, ainv.row3.y(), 1e-12);
+    try std.testing.expectApproxEqAbs(1.0 / 9.0, ainv.row3.z(), 1e-12);
+}
 test "vec3 mul mat3 works" {
     const b: Mat3F64 = .{
         .row1 = Vec3F64.init(0, 1, 0),