commit a9ee1cc8c2b8b3898c53a930cf59f91259c5a17e
parent 74ccbf67849b247a4c46bb76ba0cce1f09fdfe94
Author: Christian Ermann <christianermann@gmail.com>
Date: Sun, 8 Jun 2025 12:56:21 -0700
Add vector/matrix math
Diffstat:
| M | build.zig | | | 9 | +++++++++ |
| A | src/math.zig | | | 495 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
2 files changed, 504 insertions(+), 0 deletions(-)
diff --git a/build.zig b/build.zig
@@ -60,4 +60,13 @@ pub fn build(b: *std.Build) void {
// This will evaluate the `run` step rather than the default, which is "install".
const run_step = b.step("run", "Run the app");
run_step.dependOn(&run_cmd.step);
+
+ const test_step = b.step("test", "Run unit tests");
+ const tests = b.addTest(.{
+ .root_source_file = b.path("src/math.zig"),
+ .target = target,
+ .optimize = optimize,
+ });
+ const run_tests = b.addRunArtifact(tests);
+ test_step.dependOn(&run_tests.step);
}
diff --git a/src/math.zig b/src/math.zig
@@ -0,0 +1,495 @@
+const std = @import("std");
+
+const Vec4 = @Vector(4, f32);
+const Mat4 = [4]Vec4;
+
+pub fn dot(u: Vec4, v: Vec4) f32 {
+ return @reduce(.Add, u * v);
+}
+
+// Adapted from https://geometrian.com/programming/tutorials/cross-product/index.php
+pub fn vecMul(u: Vec4, v: Vec4) Vec4 {
+ const mask = @Vector(4, i32){ 1, 2, 0, 3 };
+ const t1 = u * @shuffle(f32, v, undefined, mask);
+ const t2 = v * @shuffle(f32, u, undefined, mask);
+ const t3 = t1 - t2;
+ return @shuffle(f32, t3, undefined, mask);
+}
+
+pub fn scale(u: Vec4, s: f32) Vec4 {
+ return u * @as(Vec4, @splat(s));
+}
+
+pub fn norm(u: Vec4) f32 {
+ return @sqrt(dot(u, u));
+}
+
+pub fn normalize(u: Vec4) Vec4 {
+ return u / @as(Vec4, @splat(norm(u)));
+}
+
+pub fn matMul(m: Mat4, n: Mat4) Mat4 {
+ var result: Mat4 = undefined;
+ inline for (0..4) |row| {
+ const x = @as(Vec4, @splat(m[row][0]));
+ const y = @as(Vec4, @splat(m[row][1]));
+ const z = @as(Vec4, @splat(m[row][2]));
+ const w = @as(Vec4, @splat(m[row][3]));
+ const xy = @mulAdd(Vec4, x, n[0], y * n[1]);
+ const zw = @mulAdd(Vec4, z, n[2], w * n[3]);
+ result[row] = xy + zw;
+ }
+ return result;
+}
+
+pub fn matMulVec(m: Mat4, u: Vec4) Vec4 {
+ return Vec4{
+ dot(m[0], u),
+ dot(m[1], u),
+ dot(m[2], u),
+ dot(m[3], u),
+ };
+}
+
+test "matMulVec" {
+ const u = matMulVec(
+ .{
+ .{ 2, 0, 0, 1 },
+ .{ 0, 3, 0, 1 },
+ .{ 0, 0, 4, 1 },
+ .{ 0, 0, 1, 5 },
+ },
+ .{ 1, 1, 1, 1 },
+ );
+ try std.testing.expectApproxEqAbs(3.0, u[0], 0.000001);
+ try std.testing.expectApproxEqAbs(4.0, u[1], 0.000001);
+ try std.testing.expectApproxEqAbs(5.0, u[2], 0.000001);
+ try std.testing.expectApproxEqAbs(6.0, u[3], 0.000001);
+}
+
+pub fn project(m: Mat4, u: Vec4) Vec4 {
+ const v = matMulVec(m, u);
+ return scale(v, 1.0 / v[3]);
+}
+
+pub fn lookToLH(eye: Vec4, dir: Vec4, up: Vec4) Mat4 {
+ const f = normalize(dir);
+ const s = normalize(vecMul(up, f));
+ const u = normalize(vecMul(f, s));
+ var result = Mat4{ s, u, f, .{ 0, 0, 0, 1 } };
+ result[0][3] = -dot(s, eye);
+ result[1][3] = -dot(u, eye);
+ result[2][3] = -dot(f, eye);
+ return result;
+}
+
+test "lookToLH" {
+ const eye = Vec4{ 5.0, 6.0, 7.0, 1.0 };
+ const dir = Vec4{ 0.5, 1.0, 0.0, 0.0 };
+ const up_ref = Vec4{ 0.0, 1.0, 0.0, 0.0 };
+ const view = lookToLH(eye, dir, up_ref);
+
+ const forward = normalize(dir);
+ const right = normalize(vecMul(up_ref, forward));
+ const up = normalize(vecMul(forward, right));
+
+ const eye_cs = matMulVec(view, eye);
+ try std.testing.expectApproxEqAbs(0.0, eye_cs[0], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, eye_cs[1], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, eye_cs[2], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, eye_cs[3], 0.000001);
+
+ const forward_cs = matMulVec(view, eye + forward);
+ try std.testing.expectApproxEqAbs(0.0, forward_cs[0], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, forward_cs[1], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, forward_cs[2], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, forward_cs[3], 0.000001);
+ const backward_cs = matMulVec(view, eye - forward);
+ try std.testing.expectApproxEqAbs(0.0, backward_cs[0], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, backward_cs[1], 0.000001);
+ try std.testing.expectApproxEqAbs(-1.0, backward_cs[2], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, backward_cs[3], 0.000001);
+
+ const right_cs = matMulVec(view, eye + right);
+ try std.testing.expectApproxEqAbs(1.0, right_cs[0], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, right_cs[1], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, right_cs[2], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, right_cs[3], 0.000001);
+ const left_cs = matMulVec(view, eye - right);
+ try std.testing.expectApproxEqAbs(-1.0, left_cs[0], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, left_cs[1], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, left_cs[2], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, left_cs[3], 0.000001);
+
+ const up_cs = matMulVec(view, eye + up);
+ try std.testing.expectApproxEqAbs(0.0, up_cs[0], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, up_cs[1], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, up_cs[2], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, up_cs[3], 0.000001);
+ const down_cs = matMulVec(view, eye - up);
+ try std.testing.expectApproxEqAbs(0.0, down_cs[0], 0.000001);
+ try std.testing.expectApproxEqAbs(-1.0, down_cs[1], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, down_cs[2], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, down_cs[3], 0.000001);
+}
+
+pub fn invLookToLH(eye: Vec4, dir: Vec4, up: Vec4) Mat4 {
+ const f = normalize(dir);
+ const s = normalize(vecMul(up, f));
+ const u = normalize(vecMul(f, s));
+ return .{
+ .{ s[0], u[0], f[0], eye[0] },
+ .{ s[1], u[1], f[1], eye[1] },
+ .{ s[2], u[2], f[2], eye[2] },
+ .{ 0, 0, 0, 1 },
+ };
+}
+
+test "invLookToLH" {
+ const eye = Vec4{ 5.0, 6.0, 7.0, 1.0 };
+ const dir = Vec4{ 0.5, 1.0, 0.0, 0.0 };
+ const up = Vec4{ 0.0, 1.0, 0.0, 0.0 };
+ const view = lookToLH(eye, dir, up);
+ const inv_view = invLookToLH(eye, dir, up);
+ const result = matMul(view, inv_view);
+
+ try std.testing.expectApproxEqAbs(1.0, result[0][0], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[0][1], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[0][2], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[0][3], 0.000001);
+
+ try std.testing.expectApproxEqAbs(0.0, result[1][0], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, result[1][1], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[1][2], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[1][3], 0.000001);
+
+ try std.testing.expectApproxEqAbs(0.0, result[2][0], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[2][1], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, result[2][2], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[2][3], 0.000001);
+
+ try std.testing.expectApproxEqAbs(0.0, result[3][0], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[3][1], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[3][2], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, result[3][3], 0.000001);
+}
+
+pub fn lookAtLH(eye: Vec4, target: Vec4, up: Vec4) Mat4 {
+ return lookToLH(eye, target - eye, up);
+}
+
+fn perspectiveLH(fovy: f32, aspect: f32, near: f32, far: f32) Mat4 {
+ // we'll split the formation of our perspective matrix into 5 steps:
+ // 1. translate the frustum apex to (0, 0, 0)
+ // 2. map depth to (NDC_Z_MIN, NDC_Z_MAX) or (NDC_Z_MAX, NDC_Z_MIN)
+ // 3. project onto the near plane
+ // 4. map x to (NDC_X_MIN, NDC_X_MAX) and y to (NDC_Y_MIN, NDC_Y_MAX)
+ // 5. combine transformations
+ //
+ // we can skip step 1 since our parameterization can only describe a
+ // frustum with its apex at the origin.
+ //
+ // for step 2, we'll use NDC_Z_MIN=0 and NDC_Z_MAX=1. this works for
+ // metal, directx, and wgpu. for opengl and vulkan, you'll either need to
+ // change your configuration or re-calculate the coefficients used in the
+ // non-linear depth mapping function with NDC_Z_MIN=-1 and NDC_Z_MAX=1. we
+ // will also choose to map the near plane to NDC_Z_MAX and the far plane to
+ // NDC_Z_MIN. this leads to a nicer distribution of values within the
+ // z-buffer if you're dealing with wide ranges of depths. you'll need to
+ // update your configuration or re-calculate the coefficients using a
+ // different mapping.
+ //
+ // for step 3, you could project onto a different plane. i've heard of
+ // people doing this but never done it myself. i'm not aware of what the
+ // advantages or disadvantages would be.
+ //
+ // for step 4, NDC_X/Y_MIN=-1 and NDC_X/Y_MAX=1. i'm not aware of any APIs
+ // that use different values, but if you use one you'll need to recalculate
+ // the depth mapping coefficients.
+ //
+ // in all of the steps, i've assumed a left-handed camera coordinate system
+ // with +z pointing into the screen. to adapt this approach to a
+ // right-handed camera coordinate system, you'll need to be careful about
+ // minus signs anywhere z is used.
+ //
+ // 2. non-linear depth mapping
+ //
+ // d(z) = (a/z) + b
+ // d(z) = (a/z) + b*(z/z)
+ // d(z) = (a + b*z) / z
+ //
+ // as a matrix:
+ // | 1 0 0 0 |
+ // | 0 1 0 0 |
+ // | 0 0 b a |
+ // | 0 0 1 0 |
+ //
+ // this matrix takes advantage of the homogenous divide performed
+ // automatically by most graphics pipelines by setting w=z.
+ //
+ // we want to use a reversed z-buffer so:
+ // d(near) = 1 d(far) = 0
+ // 1 = (a/near)+b 0 = (a/far)+b
+ // 1 = (a/near)-(a/far)
+ // 1 = a(1/near-1/far)
+ // a = 1/(1/near-1/far)
+ // a = near*far/(far-near)
+ // b = -(near*far/(far-near))/far
+ // b = near/(near-far)
+ //
+ // we'll split the formation of our "perspective" matrix into 5 steps:
+ // 1. translate the frustum apex to (0, 0, 0)
+ // 2. map depth to (NDC_Z_MIN, NDC_Z_MAX) or (NDC_Z_MAX, NDC_Z_MIN)
+ // 3. project onto the near plane
+ // 4. map x to (NDC_X_MIN, NDC_X_MAX) and y to (NDC_Y_MIN, NDC_Y_MAX)
+ // 5. combine transformations
+ //
+ // we can skip step 1 since our parameterization can only describe a
+ // frustum with its apex at the origin.
+ //
+ // for step 2, we'll use NDC_Z_MIN=0 and NDC_Z_MAX=1. this works for
+ // metal, directx, and wgpu. for opengl and vulkan, you'll either need to
+ // change your configuration or re-calculate the coefficients used in the
+ // non-linear depth mapping function with NDC_Z_MIN=-1 and NDC_Z_MAX=1. we
+ // will also choose to map the near plane to NDC_Z_MAX and the far plane to
+ // NDC_Z_MIN. this leads to a nicer distribution of values within the
+ // z-buffer if you're dealing with wide ranges of depths. you'll need to
+ // update your configuration or re-calculate the coefficients using a
+ // different mapping.
+ //
+ // for step 3, you could project onto a different plane. i've heard of
+ // people doing this but never done it myself. i'm not aware of what the
+ // advantages or disadvantages would be.
+ //
+ // for step 4, NDC_X/Y_MIN=-1 and NDC_X/Y_MAX=1. i'm not aware of any APIs
+ // that use different values, but if you use one you'll need to recalculate
+ // the depth mapping coefficients.
+ //
+ // in all of the steps, i've assumed a left-handed camera coordinate system
+ // with +z pointing into the screen. to adapt this approach to a
+ // right-handed camera coordinate system, you'll need to be careful about
+ // minus signs anywhere z is used.
+ //
+ // 2. non-linear depth mapping
+ //
+ // d(z) = (a/z) + b
+ // d(z) = (a/z) + b*(z/z)
+ // d(z) = (a + b*z) / z
+ //
+ // as a matrix:
+ // | 1 0 0 0 |
+ // | 0 1 0 0 |
+ // | 0 0 b a |
+ // | 0 0 1 0 |
+ //
+ // this matrix takes advantage of the homogenous divide performed
+ // automatically by most graphics pipelines by setting w=z.
+ //
+ // we want to use a reversed z-buffer so:
+ // d(near) = 1 d(far) = 0
+ // 1 = (a/near)+b 0 = (a/far)+b
+ // 1 = (a/near)-(a/far)
+ // 1 = a(1/near-1/far)
+ // a = 1/(1/near-1/far)
+ // a = near*far/(far-near)
+ // b = -(near*far/(far-near))/far
+ // b = near/(near-far)
+ //
+ // 3. near plane projection:
+ // ^ y
+ // ^ y' |
+ // | |
+ // ------------------------------>|
+ // camera near z
+ //
+ // this projection maps (x, y, z) to (x', y', near). the diagram is
+ // intended to show that the values of y and y' are related by similar
+ // right triangles. this means that
+ // y/z = y'/near
+ // y' = (y*near)/z
+ //
+ // as a matrix:
+ // | near 0 0 0 |
+ // | 0 near 0 0 |
+ // | 0 0 1 0 |
+ // | 0 0 1 0 |
+ //
+ // we utilize the homogenous divide here too, just like we did for our
+ // depth mapping.
+ //
+ // 4. scale x and y
+ //
+ // h = near * tan(fovy/2)
+ // w = near * tan(fovy/2) * aspect
+ //
+ // x' = x / w
+ // y' = y / h
+ //
+ // as a matrix:
+ // | 1/w 0 0 0 |
+ // | 0 1/h 0 0 |
+ // | 0 0 1 0 |
+ // | 0 0 0 1 |
+ //
+ // 5. combine transformations
+ //
+ // combining our previous steps gives us the matrix:
+ // | near/w 0 0 0 |
+ // | 0 near/h 0 0 |
+ // | 0 0 b a |
+ // | 0 0 1 0 |
+ //
+ // and now we can simplify the x and y terms:
+ // near/w = 1 / (tan(fovy/2) * aspect)
+ // near/h = 1 / tan(fovy/2)
+ //
+ // 3. near plane projection:
+ // ^ y
+ // ^ y' |
+ // | |
+ // ------------------------------>|
+ // camera near z
+ //
+ // this projection maps (x, y, z) to (x', y', near). the diagram is
+ // intended to show that the values of y and y' are related by similar
+ // right triangles. this means that
+ // y/z = y'/near
+ // y' = (y*near)/z
+ //
+ // as a matrix:
+ // | near 0 0 0 |
+ // | 0 near 0 0 |
+ // | 0 0 1 0 |
+ // | 0 0 1 0 |
+ //
+ // we utilize the homogenous divide here too, just like we did for our
+ // depth mapping.
+ //
+ // 4. scale x and y
+ //
+ // h = near * tan(fovy/2)
+ // w = near * tan(fovy/2) * aspect
+ //
+ // x' = x / w
+ // y' = y / h
+ //
+ // as a matrix:
+ // | 1/w 0 0 0 |
+ // | 0 1/h 0 0 |
+ // | 0 0 1 0 |
+ // | 0 0 0 1 |
+ //
+ // 5. combine transformations
+ //
+ // combining our previous steps gives us the matrix:
+ // | near/w 0 0 0 |
+ // | 0 near/h 0 0 |
+ // | 0 0 b a |
+ // | 0 0 1 0 |
+ //
+ // and now we can simplify the x and y terms:
+ // near/w = 1 / (tan(fovy/2) * aspect)
+ // near/h = 1 / tan(fovy/2)
+ //
+ // the implementation is pretty simple:
+ const a = near * far / (far - near);
+ const b = near / (near - far);
+ const tan_half_fov = @tan(std.math.degreesToRadians(fovy * 0.5));
+ const x = 1.0 / (aspect * tan_half_fov);
+ const y = 1.0 / tan_half_fov;
+ return Mat4{
+ .{ x, 0, 0, 0 },
+ .{ 0, y, 0, 0 },
+ .{ 0, 0, b, a },
+ .{ 0, 0, 1, 0 },
+ };
+}
+
+test "perspectiveLH" {
+ const fovy = 60.0;
+ const aspect = 16.0 / 9.0;
+ const near = 0.1;
+ const far = 10.0;
+ const proj = perspectiveLH(fovy, aspect, near, far);
+ const tan_half_fov = @tan(std.math.degreesToRadians(fovy * 0.5));
+
+ const near_h = tan_half_fov * near;
+ const near_w = near_h * aspect;
+ const near_c = project(proj, Vec4{ 0.0, 0.0, near, 1.0 });
+ try std.testing.expectApproxEqAbs(0.0, near_c[0], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, near_c[1], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, near_c[2], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, near_c[3], 0.000001);
+ const near_tr = project(proj, Vec4{ near_w, near_h, near, 1.0 });
+ try std.testing.expectApproxEqAbs(1.0, near_tr[0], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, near_tr[1], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, near_tr[2], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, near_tr[3], 0.000001);
+ const near_bl = project(proj, Vec4{ -near_w, -near_h, near, 1.0 });
+ try std.testing.expectApproxEqAbs(-1.0, near_bl[0], 0.000001);
+ try std.testing.expectApproxEqAbs(-1.0, near_bl[1], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, near_bl[2], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, near_bl[3], 0.000001);
+
+ const far_h = tan_half_fov * far;
+ const far_w = far_h * aspect;
+ const far_c = project(proj, Vec4{ 0.0, 0.0, far, 1.0 });
+ try std.testing.expectApproxEqAbs(0.0, far_c[0], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, far_c[1], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, far_c[2], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, far_c[3], 0.000001);
+ const far_tr = project(proj, Vec4{ far_w, far_h, far, 1.0 });
+ try std.testing.expectApproxEqAbs(1.0, far_tr[0], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, far_tr[1], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, far_tr[2], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, far_tr[3], 0.000001);
+ const far_bl = project(proj, Vec4{ -far_w, -far_h, far, 1.0 });
+ try std.testing.expectApproxEqAbs(-1.0, far_bl[0], 0.000001);
+ try std.testing.expectApproxEqAbs(-1.0, far_bl[1], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, far_tr[2], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, far_tr[3], 0.000001);
+}
+
+fn invPerspectiveLH(fovy: f32, aspect: f32, near: f32, far: f32) Mat4 {
+ const inv_a = (far - near) / (near * far);
+ const b = near / (near - far);
+ const tan_half_fov = @tan(std.math.degreesToRadians(fovy * 0.5));
+ const inv_x = aspect * tan_half_fov;
+ const inv_y = tan_half_fov;
+ return Mat4{
+ .{ inv_x, 0, 0, 0 },
+ .{ 0, inv_y, 0, 0 },
+ .{ 0, 0, 0, 1 },
+ .{ 0, 0, inv_a, -inv_a * b },
+ };
+}
+
+test "invPerspectiveLH" {
+ const fovy = 60.0;
+ const aspect = 16.0 / 9.0;
+ const near = 0.1;
+ const far = 10.0;
+ const proj = perspectiveLH(fovy, aspect, near, far);
+ const inv_proj = invPerspectiveLH(fovy, aspect, near, far);
+ const result = matMul(proj, inv_proj);
+
+ try std.testing.expectApproxEqAbs(1.0, result[0][0], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[0][1], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[0][2], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[0][3], 0.000001);
+
+ try std.testing.expectApproxEqAbs(0.0, result[1][0], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, result[1][1], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[1][2], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[1][3], 0.000001);
+
+ try std.testing.expectApproxEqAbs(0.0, result[2][0], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[2][1], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, result[2][2], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[2][3], 0.000001);
+
+ try std.testing.expectApproxEqAbs(0.0, result[3][0], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[3][1], 0.000001);
+ try std.testing.expectApproxEqAbs(0.0, result[3][2], 0.000001);
+ try std.testing.expectApproxEqAbs(1.0, result[3][3], 0.000001);
+}
