diff --git a/rtchallenge/src/lib.rs b/rtchallenge/src/lib.rs
index afabcff..58b9c78 100644
--- a/rtchallenge/src/lib.rs
+++ b/rtchallenge/src/lib.rs
@@ -1,2 +1,3 @@
pub mod canvas;
+pub mod matrices;
pub mod tuples;
diff --git a/rtchallenge/src/matrices.rs b/rtchallenge/src/matrices.rs
new file mode 100644
index 0000000..3d3c35f
--- /dev/null
+++ b/rtchallenge/src/matrices.rs
@@ -0,0 +1,261 @@
+use std::fmt;
+use std::ops::{Add, Div, Index, Mul, Neg, Sub};
+
+#[derive(Default, Clone, Copy)]
+/// Matrix4x4 represents a 4x4 matrix in row-major form. So, element `m[i][j]` corresponds to mi,j
+/// where `i` is the row number and `j` is the column number.
+pub struct Matrix4x4 {
+ m: [[f32; 4]; 4],
+}
+
+impl From<[f32; 16]> for Matrix4x4 {
+ fn from(t: [f32; 16]) -> Self {
+ Matrix4x4 {
+ m: [
+ [t[0], t[1], t[2], t[3]],
+ [t[4], t[5], t[6], t[7]],
+ [t[8], t[9], t[10], t[11]],
+ [t[12], t[13], t[14], t[15]],
+ ],
+ }
+ }
+}
+
+impl Matrix4x4 {
+ /// Create a `Matrix4x4` containing the identity, all zeros with ones along the diagonal.
+ pub fn identity() -> Matrix4x4 {
+ Matrix4x4::new(
+ [1., 0., 0., 0.],
+ [0., 1., 0., 0.],
+ [0., 0., 1., 0.],
+ [0., 0., 0., 1.],
+ )
+ }
+
+ /// Create a `Matrix4x4` with each of the given rows.
+ pub fn new(r0: [f32; 4], r1: [f32; 4], r2: [f32; 4], r3: [f32; 4]) -> Matrix4x4 {
+ Matrix4x4 {
+ m: [r0, r1, r2, r3],
+ }
+ }
+
+ /// Transpose self, returning a new matrix that has been reflected across the diagonal.
+ /// # Examples
+ ///
+ /// ```
+ /// use rtchallenge::matrices::Matrix4x4;
+ ///
+ /// let m = Matrix4x4::new(
+ /// [2., 0., 0., 0.],
+ /// [3., 1., 0., 0.],
+ /// [4., 0., 1., 0.],
+ /// [5., 6., 7., 1.],
+ /// );
+ /// let m_t = Matrix4x4::new(
+ /// [2., 3., 4., 5.],
+ /// [0., 1., 0., 6.],
+ /// [0., 0., 1., 7.],
+ /// [0., 0., 0., 1.],
+ /// );
+ /// assert_eq!(m.transpose(), m_t);
+ pub fn transpose(&self) -> Matrix4x4 {
+ let m = self.m;
+ Matrix4x4 {
+ m: [
+ [m[0][0], m[1][0], m[2][0], m[3][0]],
+ [m[0][1], m[1][1], m[2][1], m[3][1]],
+ [m[0][2], m[1][2], m[2][2], m[3][2]],
+ [m[0][3], m[1][3], m[2][3], m[3][3]],
+ ],
+ }
+ }
+
+ /// Returns a new matrix that is the inverse of self. If self is A, inverse returns A-1, where
+ /// AA-1 = I.
+ /// This implementation uses a numerically stable Gauss–Jordan elimination routine to compute the inverse.
+ ///
+ /// # Examples
+ ///
+ /// ```
+ /// use rtchallenge::matrices::Matrix4x4;
+ ///
+ /// let i = Matrix4x4::identity();
+ /// assert_eq!(i.inverse() * i, i);
+ ///
+ /// let m = Matrix4x4::new(
+ /// [2., 0., 0., 0.],
+ /// [0., 3., 0., 0.],
+ /// [0., 0., 4., 0.],
+ /// [0., 0., 0., 1.],
+ /// );
+ /// assert_eq!(m.inverse() * m, i);
+ /// assert_eq!(m * m.inverse(), i);
+ /// ```
+ pub fn inverse(&self) -> Matrix4x4 {
+ // TODO(wathiede): how come the C++ version doesn't need to deal with non-invertable
+ // matrix.
+ let mut indxc: [usize; 4] = Default::default();
+ let mut indxr: [usize; 4] = Default::default();
+ let mut ipiv: [usize; 4] = Default::default();
+ let mut minv = self.m;
+
+ for i in 0..4 {
+ let mut irow: usize = 0;
+ let mut icol: usize = 0;
+ let mut big: f32 = 0.;
+ // Choose pivot
+ for j in 0..4 {
+ if ipiv[j] != 1 {
+ for (k, ipivk) in ipiv.iter().enumerate() {
+ if *ipivk == 0 {
+ if minv[j][k].abs() >= big {
+ big = minv[j][k].abs();
+ irow = j;
+ icol = k;
+ }
+ } else if *ipivk > 1 {
+ eprintln!("Singular matrix in MatrixInvert");
+ }
+ }
+ }
+ }
+ ipiv[icol] += 1;
+ // Swap rows _irow_ and _icol_ for pivot
+ if irow != icol {
+ // Can't figure out how to make swap work here.
+ #[allow(clippy::manual_swap)]
+ for k in 0..4 {
+ let tmp = minv[irow][k];
+ minv[irow][k] = minv[icol][k];
+ minv[icol][k] = tmp;
+ }
+ }
+ indxr[i] = irow;
+ indxc[i] = icol;
+ if minv[icol][icol] == 0. {
+ eprintln!("Singular matrix in MatrixInvert");
+ }
+
+ // Set $m[icol][icol]$ to one by scaling row _icol_ appropriately
+ let pivinv: f32 = minv[icol][icol].recip();
+ minv[icol][icol] = 1.;
+ for j in 0..4 {
+ minv[icol][j] *= pivinv;
+ }
+
+ // Subtract this row from others to zero out their columns
+ for j in 0..4 {
+ if j != icol {
+ let save = minv[j][icol];
+ minv[j][icol] = 0.;
+ for k in 0..4 {
+ minv[j][k] -= minv[icol][k] * save;
+ }
+ }
+ }
+ }
+ // Swap columns to reflect permutation
+ for j in (0..4).rev() {
+ if indxr[j] != indxc[j] {
+ for mi in &mut minv {
+ mi.swap(indxr[j], indxc[j])
+ }
+ }
+ }
+ Matrix4x4 { m: minv }
+ }
+}
+
+impl fmt::Debug for Matrix4x4 {
+ fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+ if f.alternate() {
+ write!(
+ f,
+ "{:?}\n {:?}\n {:?}\n {:?}",
+ self.m[0], self.m[1], self.m[2], self.m[3]
+ )
+ } else {
+ write!(
+ f,
+ "[{:?} {:?} {:?} {:?}]",
+ self.m[0], self.m[1], self.m[2], self.m[3]
+ )
+ }
+ }
+}
+
+impl Mul for Matrix4x4 {
+ type Output = Matrix4x4;
+
+ /// Implement matrix multiplication for `Matrix4x4`.
+ ///
+ /// # Examples
+ /// ```
+ /// use rtchallenge::matrices::Matrix4x4;
+ ///
+ /// let i = Matrix4x4::identity();
+ /// let m1 = Matrix4x4::identity();
+ /// let m2 = Matrix4x4::identity();
+ ///
+ /// assert_eq!(m1 * m2, i);
+ /// ```
+ fn mul(self, m2: Matrix4x4) -> Matrix4x4 {
+ let m1 = self;
+ let mut r: Matrix4x4 = Default::default();
+ for i in 0..4 {
+ for j in 0..4 {
+ r.m[i][j] = m1.m[i][0] * m2.m[0][j]
+ + m1.m[i][1] * m2.m[1][j]
+ + m1.m[i][2] * m2.m[2][j]
+ + m1.m[i][3] * m2.m[3][j];
+ }
+ }
+ r
+ }
+}
+
+impl PartialEq for Matrix4x4 {
+ fn eq(&self, rhs: &Matrix4x4) -> bool {
+ let l = self.m;
+ let r = rhs.m;
+ for i in 0..4 {
+ for j in 0..4 {
+ let d = (l[i][j] - r[i][j]).abs();
+ if d > f32::EPSILON {
+ return false;
+ }
+ }
+ }
+ true
+ }
+}
+
+impl Index<(usize, usize)> for Matrix4x4 {
+ type Output = f32;
+ fn index(&self, (row, col): (usize, usize)) -> &Self::Output {
+ &self.m[row][col]
+ }
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+
+ #[test]
+ fn construct4x4() {
+ let m = Matrix4x4::new(
+ [1., 2., 3., 4.],
+ [5.5, 6.5, 7.5, 8.5],
+ [9., 10., 11., 12.],
+ [13.5, 14.5, 15.5, 16.5],
+ );
+
+ assert_eq!(m[(0, 0)], 1.);
+ assert_eq!(m[(0, 3)], 4.);
+ assert_eq!(m[(1, 0)], 5.5);
+ assert_eq!(m[(1, 2)], 7.5);
+ assert_eq!(m[(2, 2)], 11.);
+ assert_eq!(m[(3, 0)], 13.5);
+ assert_eq!(m[(3, 2)], 15.5);
+ }
+}