lars/matrix/
mat3.rs

1//! 3×3 Matrix utilities.
2//!
3//! Provides a small, self-contained 3×3 matrix type [`Mat3`] with
4//! support for basic linear algebra operations, including addition,
5//! subtraction, scalar and matrix multiplication, and inversion.
6//!
7//! This type is designed to pair naturally with the [`Vec3] struct
8//! for 3D linear transformations.
9
10use std::ops::Mul;
11use derive_more::{Constructor, Add, Sub, Div};
12use crate::Vec3;
13
14/// a 3×3 matrix of `f64` values.
15///
16/// The matrix is stored in **row-major order**:
17///
18/// ```text
19/// | a  b  c |
20/// | c  d  e |
21/// | g  h  i |
22/// ```
23///
24
25/// ```
26#[derive(Constructor, Copy, Clone, Debug, Add, Sub, PartialOrd, Div)]
27pub struct Mat3 {
28    /// First row, first column element.
29    pub a: f64,
30    /// First row, second column element.
31    pub b: f64,
32    /// First row, third column element.
33    pub c: f64,
34    /// Second row, first column element.
35    pub d: f64,
36    /// Second row, second column element.
37    pub e: f64,
38    /// Second row, third column element.
39    pub f: f64,
40    /// Third row, first column element.
41    pub g: f64,
42    /// Third row, second column element.
43    pub h: f64,
44    /// Third row, third column element.
45    pub i: f64,
46}
47
48impl Mat3 {
49    /// The **identity matrix**:
50    ///
51    /// ```text
52    /// | 1  0  0 |
53    /// | 0  1  0 |
54    /// | 0  0  1 |
55    /// ```
56    pub const IDENTITY: Mat3 = Mat3 {
57        a: 1.0,
58        b: 0.0,
59        c: 0.0,
60        d: 0.0,
61        e: 1.0,
62        f: 0.0,
63        g: 0.0,
64        h: 0.0,
65        i: 1.0,
66    };
67    /// The **zero matrix**:
68    ///
69    /// ```text
70    /// | 0  0  0 |
71    /// | 0  0  0 |
72    /// | 0  0  0 |
73    /// ```
74    pub const ZERO: Mat3 = Mat3 {
75        a: 0.0,
76        b: 0.0,
77        c: 0.0,
78        d: 0.0,
79        e: 0.0,
80        f: 0.0,
81        g: 0.0,
82        h: 0.0,
83        i: 0.0,
84    };
85    /// Returns the **determinant** of the matrix.
86    ///
87    /// Computed as:
88    /// \[
89    /// \det(M) = a(ei - fh) - b(di - fg) + c(dh - eg)
90    /// \]
91    ///
92    /// # Examples
93    /// ```
94    /// use lars::Mat3;
95    /// let m = Mat3::new(1.0, 2.0, 3.0, 3.0, 2.0, 1.0, 2.0, 1.0, 3.0);
96    /// assert_eq!(m.determinant(), -12.0);
97    /// ```
98    pub fn determinant(&self) -> f64 {
99        self.a * (self.e * self.i - self.f * self.h) - self.b * (self.d * self.i - self.f * self.g) + self.c * (self.d * self.h - self.e * self.g)
100    }
101    /// Returns the **inverse** of the matrix, if it exists.
102    ///
103    /// Computed as:
104    /// M⁻¹ = (1/det(M)) * adj(M)
105    /// /
106    ///          1        | ei - fh   ch - bi   bf - ce |
107    /// M⁻¹ = -------  x  | fg - di   ai - cg   cd - af |
108    ///        det(M)     | dh - eg   bg - ah   ae - bd |
109    ///
110    /// # Panics
111    /// Panics if the matrix is singular (determinant = 0).
112    ///
113    /// # Examples
114    /// ```
115    /// use lars::Mat3;
116    /// let m = Mat3::new(1.0, 2.0, 3.0, 3.0, 2.0, 1.0, 2.0, 1.0, 3.0);
117    /// assert_eq!(m.inverse(), Mat3::new(-5.0, 3.0, 4.0, 7.0, 3.0, -8.0, 1.0, -3.0, 4.0)/12.0);
118    pub fn inverse(&self) -> Mat3 {
119        let det = self.determinant();
120        if det == 0.0 {
121            panic!("Matrix is singular and cannot be inverted.");
122        }
123        let inv_det = 1.0 / det;
124
125        Mat3 {
126            a: (self.e * self.i - self.f * self.h) * inv_det,
127            b: (self.c * self.h - self.b * self.i) * inv_det,
128            c: (self.b * self.f - self.c * self.e) * inv_det,
129            d: (self.f * self.g - self.d * self.i) * inv_det,
130            e: (self.a * self.i - self.c * self.g) * inv_det,
131            f: (self.c * self.d - self.a * self.f) * inv_det,
132            g: (self.d * self.h - self.e * self.g) * inv_det,
133            h: (self.b * self.g - self.a * self.h) * inv_det,
134            i: (self.a * self.e - self.b * self.d) * inv_det,
135        }
136    }
137
138}
139
140const EPSILON: f64 = 1e-9;
141
142impl PartialEq for Mat3 {
143    fn eq(&self, other: &Self) -> bool {
144        (self.a - other.a).abs() < EPSILON &&
145            (self.b - other.b).abs() < EPSILON &&
146            (self.c - other.c).abs() < EPSILON &&
147            (self.d - other.d).abs() < EPSILON &&
148            (self.e - other.e).abs() < EPSILON &&
149            (self.f - other.f).abs() < EPSILON &&
150            (self.g - other.g).abs() < EPSILON &&
151            (self.h - other.h).abs() < EPSILON &&
152            (self.i - other.i).abs() < EPSILON
153    }
154}
155
156impl Mul<Vec3> for Mat3 {
157    type Output = Vec3;
158
159    fn mul(self, rhs: Vec3) -> Vec3 {
160        Vec3 {
161            x: self.a * rhs.x + self.b * rhs.y + self.c * rhs.z,
162            y: self.d * rhs.x + self.e * rhs.y + self.f * rhs.z,
163            z: self.g * rhs.x + self.h * rhs.y + self.i * rhs.z,
164        }
165    }
166}
167
168impl Mul<Mat3> for Mat3 {
169    type Output = Mat3;
170
171    fn mul(self, rhs: Mat3) -> Mat3 {
172        Mat3 {
173            a: self.a * rhs.a + self.b * rhs.d + self.c * rhs.g,
174            b: self.a * rhs.b + self.b * rhs.e + self.c * rhs.h,
175            c: self.a * rhs.c + self.b * rhs.f + self.c * rhs.i,
176            d: self.d * rhs.a + self.e * rhs.d + self.f * rhs.g,
177            e: self.d * rhs.b + self.e * rhs.e + self.f * rhs.h,
178            f: self.d * rhs.c + self.e * rhs.f + self.f * rhs.i,
179            g: self.g * rhs.a + self.h * rhs.d + self.i * rhs.g,
180            h: self.g * rhs.b + self.h * rhs.e + self.i * rhs.h,
181            i: self.g * rhs.c + self.h * rhs.f + self.i * rhs.i,
182        }
183    }
184}
185
186#[cfg(test)]
187mod tests {
188    use super::*;
189
190    #[test]
191    fn test_add() {
192        let m = Mat3::new(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0);
193        assert_eq!(m + m, Mat3::new(2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0));
194    }
195
196    #[test]
197    fn sub() {
198        let m = Mat3::new(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0);
199        assert_eq!(m - m, Mat3::ZERO);
200    }
201
202    #[test]
203    fn test_mat_mul() {
204        let a = Mat3::IDENTITY;
205        let b = Mat3::new(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0);
206        assert_eq!(a * b, b);
207    }
208
209    #[test]
210    fn test_determinant() {
211        let m = Mat3::new(1.0, 2.0, 3.0, 3.0, 2.0, 1.0, 2.0, 1.0, 3.0);
212        assert_eq!(m.determinant(), -12.0);
213    }
214
215    #[test]
216    fn test_inverse() {
217        let m = Mat3::new(1.0, 2.0, 3.0, 3.0, 2.0, 1.0, 2.0, 1.0, 3.0);
218        assert_eq!(m.inverse(), Mat3::new(-5.0, 3.0, 4.0, 7.0, 3.0, -8.0, 1.0, -3.0, 4.0)/12.0);
219    }
220}