use_matrix/
lib.rs

1#![forbid(unsafe_code)]
2#![doc = include_str!("../README.md")]
3
4//! Small matrix primitives and operations for `RustUse`.
5
6#[doc(inline)]
7pub use matrix2::Matrix2;
8#[doc(inline)]
9pub use matrix3::Matrix3;
10#[doc(inline)]
11pub use matrix4::Matrix4;
12
13/// Two-dimensional matrix primitives and operations.
14pub mod matrix2 {
15    use core::num::FpCategory;
16    use core::ops::{Add, Div, Mul, Neg, Sub};
17
18    use use_vector::Vector2;
19
20    #[inline]
21    fn dot2(a0: f64, b0: f64, a1: f64, b1: f64) -> f64 {
22        a0.mul_add(b0, a1 * b1)
23    }
24
25    /// A 2x2 matrix stored in row-major order.
26    ///
27    /// # Examples
28    ///
29    /// ```
30    /// use use_matrix::Matrix2;
31    ///
32    /// let matrix = Matrix2::new(
33    ///     1.0, 2.0,
34    ///     3.0, 4.0,
35    /// );
36    ///
37    /// assert_eq!(matrix.trace(), 5.0);
38    /// assert_eq!(matrix.determinant(), -2.0);
39    /// ```
40    #[derive(Debug, Clone, Copy, PartialEq)]
41    pub struct Matrix2 {
42        pub m00: f64,
43        pub m01: f64,
44        pub m10: f64,
45        pub m11: f64,
46    }
47
48    impl Matrix2 {
49        /// The zero matrix.
50        pub const ZERO: Self = Self::new(0.0, 0.0, 0.0, 0.0);
51
52        /// The identity matrix.
53        ///
54        /// # Examples
55        ///
56        /// ```
57        /// use use_matrix::Matrix2;
58        ///
59        /// let matrix = Matrix2::new(
60        ///     4.0, 7.0,
61        ///     2.0, 6.0,
62        /// );
63        ///
64        /// assert_eq!(Matrix2::IDENTITY * matrix, matrix);
65        /// assert_eq!(matrix * Matrix2::IDENTITY, matrix);
66        /// ```
67        pub const IDENTITY: Self = Self::new(1.0, 0.0, 0.0, 1.0);
68
69        /// Creates a matrix from row-major entries.
70        #[must_use]
71        #[allow(clippy::similar_names, clippy::too_many_arguments)]
72        pub const fn new(m00: f64, m01: f64, m10: f64, m11: f64) -> Self {
73            Self { m00, m01, m10, m11 }
74        }
75
76        /// Creates a matrix from row arrays.
77        #[must_use]
78        pub const fn from_rows(rows: [[f64; 2]; 2]) -> Self {
79            Self::new(rows[0][0], rows[0][1], rows[1][0], rows[1][1])
80        }
81
82        /// Creates a matrix from column arrays.
83        #[must_use]
84        pub const fn from_cols(cols: [[f64; 2]; 2]) -> Self {
85            Self::new(cols[0][0], cols[1][0], cols[0][1], cols[1][1])
86        }
87
88        /// Returns the transpose.
89        #[must_use]
90        pub const fn transpose(self) -> Self {
91            Self::new(self.m00, self.m10, self.m01, self.m11)
92        }
93
94        /// Returns the determinant.
95        #[must_use]
96        pub const fn determinant(self) -> f64 {
97            (self.m00 * self.m11) - (self.m01 * self.m10)
98        }
99
100        /// Returns the trace.
101        #[must_use]
102        pub const fn trace(self) -> f64 {
103            self.m00 + self.m11
104        }
105
106        /// Returns the inverse when the determinant is finite and non-zero.
107        ///
108        /// Returns `None` for singular matrices or matrices with non-finite
109        /// determinants.
110        ///
111        /// # Examples
112        ///
113        /// ```
114        /// use use_matrix::Matrix2;
115        ///
116        /// let matrix = Matrix2::new(
117        ///     4.0, 7.0,
118        ///     2.0, 6.0,
119        /// );
120        ///
121        /// let inverse = matrix.inverse().expect("matrix should invert");
122        ///
123        /// assert!((inverse.m00 - 0.6).abs() < 1.0e-12);
124        /// assert!((inverse.m11 - 0.4).abs() < 1.0e-12);
125        /// ```
126        #[must_use]
127        pub fn inverse(self) -> Option<Self> {
128            let determinant = self.determinant();
129
130            if !determinant.is_finite() || matches!(determinant.classify(), FpCategory::Zero) {
131                return None;
132            }
133
134            Some(Self::new(
135                self.m11 / determinant,
136                -self.m01 / determinant,
137                -self.m10 / determinant,
138                self.m00 / determinant,
139            ))
140        }
141    }
142
143    impl Add for Matrix2 {
144        type Output = Self;
145
146        fn add(self, rhs: Self) -> Self::Output {
147            Self::new(
148                self.m00 + rhs.m00,
149                self.m01 + rhs.m01,
150                self.m10 + rhs.m10,
151                self.m11 + rhs.m11,
152            )
153        }
154    }
155
156    impl Sub for Matrix2 {
157        type Output = Self;
158
159        fn sub(self, rhs: Self) -> Self::Output {
160            Self::new(
161                self.m00 - rhs.m00,
162                self.m01 - rhs.m01,
163                self.m10 - rhs.m10,
164                self.m11 - rhs.m11,
165            )
166        }
167    }
168
169    impl Mul<f64> for Matrix2 {
170        type Output = Self;
171
172        fn mul(self, rhs: f64) -> Self::Output {
173            Self::new(
174                self.m00 * rhs,
175                self.m01 * rhs,
176                self.m10 * rhs,
177                self.m11 * rhs,
178            )
179        }
180    }
181
182    impl Div<f64> for Matrix2 {
183        type Output = Self;
184
185        fn div(self, rhs: f64) -> Self::Output {
186            Self::new(
187                self.m00 / rhs,
188                self.m01 / rhs,
189                self.m10 / rhs,
190                self.m11 / rhs,
191            )
192        }
193    }
194
195    impl Neg for Matrix2 {
196        type Output = Self;
197
198        fn neg(self) -> Self::Output {
199            Self::new(-self.m00, -self.m01, -self.m10, -self.m11)
200        }
201    }
202
203    impl Mul for Matrix2 {
204        type Output = Self;
205
206        fn mul(self, rhs: Self) -> Self::Output {
207            Self::new(
208                dot2(self.m00, rhs.m00, self.m01, rhs.m10),
209                dot2(self.m00, rhs.m01, self.m01, rhs.m11),
210                dot2(self.m10, rhs.m00, self.m11, rhs.m10),
211                dot2(self.m10, rhs.m01, self.m11, rhs.m11),
212            )
213        }
214    }
215
216    impl Mul<Vector2> for Matrix2 {
217        type Output = Vector2;
218
219        fn mul(self, rhs: Vector2) -> Self::Output {
220            let x_component = rhs.x();
221            let y_component = rhs.y();
222
223            Vector2::new(
224                dot2(self.m00, x_component, self.m01, y_component),
225                dot2(self.m10, x_component, self.m11, y_component),
226            )
227        }
228    }
229}
230
231/// Three-dimensional matrix primitives and operations.
232pub mod matrix3 {
233    use core::num::FpCategory;
234    use core::ops::{Add, Div, Mul, Neg, Sub};
235
236    use use_vector::Vector3;
237
238    #[inline]
239    fn dot3(a0: f64, b0: f64, a1: f64, b1: f64, a2: f64, b2: f64) -> f64 {
240        a0.mul_add(b0, a1.mul_add(b1, a2 * b2))
241    }
242
243    #[inline]
244    fn mul_sub(a: f64, b: f64, c: f64, d: f64) -> f64 {
245        a.mul_add(b, -(c * d))
246    }
247
248    /// A 3x3 matrix stored in row-major order.
249    ///
250    /// # Examples
251    ///
252    /// ```
253    /// use use_matrix::Matrix3;
254    ///
255    /// let matrix = Matrix3::new(
256    ///     1.0, 2.0, 3.0,
257    ///     0.0, 1.0, 4.0,
258    ///     5.0, 6.0, 0.0,
259    /// );
260    ///
261    /// assert_eq!(matrix.determinant(), 1.0);
262    /// ```
263    #[derive(Debug, Clone, Copy, PartialEq)]
264    pub struct Matrix3 {
265        pub m00: f64,
266        pub m01: f64,
267        pub m02: f64,
268        pub m10: f64,
269        pub m11: f64,
270        pub m12: f64,
271        pub m20: f64,
272        pub m21: f64,
273        pub m22: f64,
274    }
275
276    impl Matrix3 {
277        /// The zero matrix.
278        pub const ZERO: Self = Self::new(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
279
280        /// The identity matrix.
281        ///
282        /// # Examples
283        ///
284        /// ```
285        /// use use_matrix::Matrix3;
286        ///
287        /// let matrix = Matrix3::new(
288        ///     2.0, 1.0, 0.0,
289        ///     0.0, 3.0, 4.0,
290        ///     0.0, 0.0, 5.0,
291        /// );
292        ///
293        /// assert_eq!(Matrix3::IDENTITY * matrix, matrix);
294        /// ```
295        pub const IDENTITY: Self = Self::new(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0);
296
297        /// Creates a matrix from row-major entries.
298        #[must_use]
299        #[allow(clippy::similar_names, clippy::too_many_arguments)]
300        pub const fn new(
301            m00: f64,
302            m01: f64,
303            m02: f64,
304            m10: f64,
305            m11: f64,
306            m12: f64,
307            m20: f64,
308            m21: f64,
309            m22: f64,
310        ) -> Self {
311            Self {
312                m00,
313                m01,
314                m02,
315                m10,
316                m11,
317                m12,
318                m20,
319                m21,
320                m22,
321            }
322        }
323
324        /// Creates a matrix from row arrays.
325        #[must_use]
326        pub const fn from_rows(rows: [[f64; 3]; 3]) -> Self {
327            Self::new(
328                rows[0][0], rows[0][1], rows[0][2], rows[1][0], rows[1][1], rows[1][2], rows[2][0],
329                rows[2][1], rows[2][2],
330            )
331        }
332
333        /// Creates a matrix from column arrays.
334        #[must_use]
335        pub const fn from_cols(cols: [[f64; 3]; 3]) -> Self {
336            Self::new(
337                cols[0][0], cols[1][0], cols[2][0], cols[0][1], cols[1][1], cols[2][1], cols[0][2],
338                cols[1][2], cols[2][2],
339            )
340        }
341
342        /// Returns the transpose.
343        #[must_use]
344        pub const fn transpose(self) -> Self {
345            Self::new(
346                self.m00, self.m10, self.m20, self.m01, self.m11, self.m21, self.m02, self.m12,
347                self.m22,
348            )
349        }
350
351        /// Returns the determinant.
352        #[must_use]
353        pub const fn determinant(self) -> f64 {
354            self.m00 * ((self.m11 * self.m22) - (self.m12 * self.m21))
355                - self.m01 * ((self.m10 * self.m22) - (self.m12 * self.m20))
356                + self.m02 * ((self.m10 * self.m21) - (self.m11 * self.m20))
357        }
358
359        /// Returns the trace.
360        #[must_use]
361        pub const fn trace(self) -> f64 {
362            self.m00 + self.m11 + self.m22
363        }
364
365        /// Returns the inverse when the determinant is finite and non-zero.
366        ///
367        /// Returns `None` for singular matrices or matrices with non-finite
368        /// determinants.
369        ///
370        /// # Examples
371        ///
372        /// ```
373        /// use use_matrix::Matrix3;
374        ///
375        /// let matrix = Matrix3::new(
376        ///     1.0, 2.0, 3.0,
377        ///     0.0, 1.0, 4.0,
378        ///     5.0, 6.0, 0.0,
379        /// );
380        ///
381        /// let inverse = matrix.inverse().expect("matrix should invert");
382        ///
383        /// assert_eq!(matrix * inverse, Matrix3::IDENTITY);
384        /// ```
385        #[must_use]
386        pub fn inverse(self) -> Option<Self> {
387            let determinant = self.determinant();
388
389            if !determinant.is_finite() || matches!(determinant.classify(), FpCategory::Zero) {
390                return None;
391            }
392
393            let c00 = mul_sub(self.m11, self.m22, self.m12, self.m21);
394            let c01 = mul_sub(self.m12, self.m20, self.m10, self.m22);
395            let c02 = mul_sub(self.m10, self.m21, self.m11, self.m20);
396            let c10 = mul_sub(self.m02, self.m21, self.m01, self.m22);
397            let c11 = mul_sub(self.m00, self.m22, self.m02, self.m20);
398            let c12 = mul_sub(self.m01, self.m20, self.m00, self.m21);
399            let c20 = mul_sub(self.m01, self.m12, self.m02, self.m11);
400            let c21 = mul_sub(self.m02, self.m10, self.m00, self.m12);
401            let c22 = mul_sub(self.m00, self.m11, self.m01, self.m10);
402
403            Some(Self::new(c00, c10, c20, c01, c11, c21, c02, c12, c22) / determinant)
404        }
405    }
406
407    impl Add for Matrix3 {
408        type Output = Self;
409
410        fn add(self, rhs: Self) -> Self::Output {
411            Self::new(
412                self.m00 + rhs.m00,
413                self.m01 + rhs.m01,
414                self.m02 + rhs.m02,
415                self.m10 + rhs.m10,
416                self.m11 + rhs.m11,
417                self.m12 + rhs.m12,
418                self.m20 + rhs.m20,
419                self.m21 + rhs.m21,
420                self.m22 + rhs.m22,
421            )
422        }
423    }
424
425    impl Sub for Matrix3 {
426        type Output = Self;
427
428        fn sub(self, rhs: Self) -> Self::Output {
429            Self::new(
430                self.m00 - rhs.m00,
431                self.m01 - rhs.m01,
432                self.m02 - rhs.m02,
433                self.m10 - rhs.m10,
434                self.m11 - rhs.m11,
435                self.m12 - rhs.m12,
436                self.m20 - rhs.m20,
437                self.m21 - rhs.m21,
438                self.m22 - rhs.m22,
439            )
440        }
441    }
442
443    impl Mul<f64> for Matrix3 {
444        type Output = Self;
445
446        fn mul(self, rhs: f64) -> Self::Output {
447            Self::new(
448                self.m00 * rhs,
449                self.m01 * rhs,
450                self.m02 * rhs,
451                self.m10 * rhs,
452                self.m11 * rhs,
453                self.m12 * rhs,
454                self.m20 * rhs,
455                self.m21 * rhs,
456                self.m22 * rhs,
457            )
458        }
459    }
460
461    impl Div<f64> for Matrix3 {
462        type Output = Self;
463
464        fn div(self, rhs: f64) -> Self::Output {
465            Self::new(
466                self.m00 / rhs,
467                self.m01 / rhs,
468                self.m02 / rhs,
469                self.m10 / rhs,
470                self.m11 / rhs,
471                self.m12 / rhs,
472                self.m20 / rhs,
473                self.m21 / rhs,
474                self.m22 / rhs,
475            )
476        }
477    }
478
479    impl Neg for Matrix3 {
480        type Output = Self;
481
482        fn neg(self) -> Self::Output {
483            Self::new(
484                -self.m00, -self.m01, -self.m02, -self.m10, -self.m11, -self.m12, -self.m20,
485                -self.m21, -self.m22,
486            )
487        }
488    }
489
490    impl Mul for Matrix3 {
491        type Output = Self;
492
493        fn mul(self, rhs: Self) -> Self::Output {
494            Self::new(
495                dot3(self.m00, rhs.m00, self.m01, rhs.m10, self.m02, rhs.m20),
496                dot3(self.m00, rhs.m01, self.m01, rhs.m11, self.m02, rhs.m21),
497                dot3(self.m00, rhs.m02, self.m01, rhs.m12, self.m02, rhs.m22),
498                dot3(self.m10, rhs.m00, self.m11, rhs.m10, self.m12, rhs.m20),
499                dot3(self.m10, rhs.m01, self.m11, rhs.m11, self.m12, rhs.m21),
500                dot3(self.m10, rhs.m02, self.m11, rhs.m12, self.m12, rhs.m22),
501                dot3(self.m20, rhs.m00, self.m21, rhs.m10, self.m22, rhs.m20),
502                dot3(self.m20, rhs.m01, self.m21, rhs.m11, self.m22, rhs.m21),
503                dot3(self.m20, rhs.m02, self.m21, rhs.m12, self.m22, rhs.m22),
504            )
505        }
506    }
507
508    impl Mul<Vector3> for Matrix3 {
509        type Output = Vector3;
510
511        fn mul(self, rhs: Vector3) -> Self::Output {
512            let x_component = rhs.x();
513            let y_component = rhs.y();
514            let z_component = rhs.z();
515
516            Vector3::new(
517                dot3(
518                    self.m00,
519                    x_component,
520                    self.m01,
521                    y_component,
522                    self.m02,
523                    z_component,
524                ),
525                dot3(
526                    self.m10,
527                    x_component,
528                    self.m11,
529                    y_component,
530                    self.m12,
531                    z_component,
532                ),
533                dot3(
534                    self.m20,
535                    x_component,
536                    self.m21,
537                    y_component,
538                    self.m22,
539                    z_component,
540                ),
541            )
542        }
543    }
544}
545
546/// Four-dimensional matrix primitives and operations.
547pub mod matrix4 {
548    use core::ops::{Add, Div, Mul, Neg, Sub};
549
550    use use_vector::Vector4;
551
552    #[inline]
553    fn dot4(a0: f64, b0: f64, a1: f64, b1: f64, a2: f64, b2: f64, a3: f64, b3: f64) -> f64 {
554        a0.mul_add(b0, a1.mul_add(b1, a2.mul_add(b2, a3 * b3)))
555    }
556
557    /// A 4x4 matrix stored in row-major order.
558    ///
559    /// # Examples
560    ///
561    /// ```
562    /// use use_matrix::Matrix4;
563    ///
564    /// let matrix = Matrix4::from_rows([
565    ///     [2.0, 1.0, 0.0, 0.0],
566    ///     [0.0, 3.0, 4.0, 0.0],
567    ///     [0.0, 0.0, 5.0, 6.0],
568    ///     [0.0, 0.0, 0.0, 7.0],
569    /// ]);
570    ///
571    /// assert_eq!(matrix.determinant(), 210.0);
572    /// ```
573    #[derive(Debug, Clone, Copy, PartialEq)]
574    pub struct Matrix4 {
575        pub m00: f64,
576        pub m01: f64,
577        pub m02: f64,
578        pub m03: f64,
579        pub m10: f64,
580        pub m11: f64,
581        pub m12: f64,
582        pub m13: f64,
583        pub m20: f64,
584        pub m21: f64,
585        pub m22: f64,
586        pub m23: f64,
587        pub m30: f64,
588        pub m31: f64,
589        pub m32: f64,
590        pub m33: f64,
591    }
592
593    impl Matrix4 {
594        /// The zero matrix.
595        pub const ZERO: Self = Self::new(
596            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
597        );
598
599        /// The identity matrix.
600        ///
601        /// # Examples
602        ///
603        /// ```
604        /// use use_matrix::Matrix4;
605        ///
606        /// let matrix = Matrix4::from_rows([
607        ///     [1.0, 2.0, 3.0, 4.0],
608        ///     [5.0, 6.0, 7.0, 8.0],
609        ///     [2.0, 0.0, 1.0, 3.0],
610        ///     [4.0, 1.0, 0.0, 2.0],
611        /// ]);
612        ///
613        /// assert_eq!(Matrix4::IDENTITY * matrix, matrix);
614        /// ```
615        pub const IDENTITY: Self = Self::new(
616            1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0,
617        );
618
619        /// Creates a matrix from row-major entries.
620        #[must_use]
621        #[allow(clippy::similar_names, clippy::too_many_arguments)]
622        pub const fn new(
623            m00: f64,
624            m01: f64,
625            m02: f64,
626            m03: f64,
627            m10: f64,
628            m11: f64,
629            m12: f64,
630            m13: f64,
631            m20: f64,
632            m21: f64,
633            m22: f64,
634            m23: f64,
635            m30: f64,
636            m31: f64,
637            m32: f64,
638            m33: f64,
639        ) -> Self {
640            Self {
641                m00,
642                m01,
643                m02,
644                m03,
645                m10,
646                m11,
647                m12,
648                m13,
649                m20,
650                m21,
651                m22,
652                m23,
653                m30,
654                m31,
655                m32,
656                m33,
657            }
658        }
659
660        /// Creates a matrix from row arrays.
661        #[must_use]
662        pub const fn from_rows(rows: [[f64; 4]; 4]) -> Self {
663            Self::new(
664                rows[0][0], rows[0][1], rows[0][2], rows[0][3], rows[1][0], rows[1][1], rows[1][2],
665                rows[1][3], rows[2][0], rows[2][1], rows[2][2], rows[2][3], rows[3][0], rows[3][1],
666                rows[3][2], rows[3][3],
667            )
668        }
669
670        /// Creates a matrix from column arrays.
671        #[must_use]
672        pub const fn from_cols(cols: [[f64; 4]; 4]) -> Self {
673            Self::new(
674                cols[0][0], cols[1][0], cols[2][0], cols[3][0], cols[0][1], cols[1][1], cols[2][1],
675                cols[3][1], cols[0][2], cols[1][2], cols[2][2], cols[3][2], cols[0][3], cols[1][3],
676                cols[2][3], cols[3][3],
677            )
678        }
679
680        /// Returns the transpose.
681        #[must_use]
682        pub const fn transpose(self) -> Self {
683            Self::new(
684                self.m00, self.m10, self.m20, self.m30, self.m01, self.m11, self.m21, self.m31,
685                self.m02, self.m12, self.m22, self.m32, self.m03, self.m13, self.m23, self.m33,
686            )
687        }
688
689        /// Returns the determinant.
690        #[must_use]
691        pub const fn determinant(self) -> f64 {
692            let minor00 = determinant3([
693                [self.m11, self.m12, self.m13],
694                [self.m21, self.m22, self.m23],
695                [self.m31, self.m32, self.m33],
696            ]);
697            let minor01 = determinant3([
698                [self.m10, self.m12, self.m13],
699                [self.m20, self.m22, self.m23],
700                [self.m30, self.m32, self.m33],
701            ]);
702            let minor02 = determinant3([
703                [self.m10, self.m11, self.m13],
704                [self.m20, self.m21, self.m23],
705                [self.m30, self.m31, self.m33],
706            ]);
707            let minor03 = determinant3([
708                [self.m10, self.m11, self.m12],
709                [self.m20, self.m21, self.m22],
710                [self.m30, self.m31, self.m32],
711            ]);
712
713            (self.m00 * minor00) - (self.m01 * minor01) + (self.m02 * minor02)
714                - (self.m03 * minor03)
715        }
716
717        /// Returns the trace.
718        #[must_use]
719        pub const fn trace(self) -> f64 {
720            self.m00 + self.m11 + self.m22 + self.m33
721        }
722    }
723
724    impl Add for Matrix4 {
725        type Output = Self;
726
727        fn add(self, rhs: Self) -> Self::Output {
728            Self::new(
729                self.m00 + rhs.m00,
730                self.m01 + rhs.m01,
731                self.m02 + rhs.m02,
732                self.m03 + rhs.m03,
733                self.m10 + rhs.m10,
734                self.m11 + rhs.m11,
735                self.m12 + rhs.m12,
736                self.m13 + rhs.m13,
737                self.m20 + rhs.m20,
738                self.m21 + rhs.m21,
739                self.m22 + rhs.m22,
740                self.m23 + rhs.m23,
741                self.m30 + rhs.m30,
742                self.m31 + rhs.m31,
743                self.m32 + rhs.m32,
744                self.m33 + rhs.m33,
745            )
746        }
747    }
748
749    impl Sub for Matrix4 {
750        type Output = Self;
751
752        fn sub(self, rhs: Self) -> Self::Output {
753            Self::new(
754                self.m00 - rhs.m00,
755                self.m01 - rhs.m01,
756                self.m02 - rhs.m02,
757                self.m03 - rhs.m03,
758                self.m10 - rhs.m10,
759                self.m11 - rhs.m11,
760                self.m12 - rhs.m12,
761                self.m13 - rhs.m13,
762                self.m20 - rhs.m20,
763                self.m21 - rhs.m21,
764                self.m22 - rhs.m22,
765                self.m23 - rhs.m23,
766                self.m30 - rhs.m30,
767                self.m31 - rhs.m31,
768                self.m32 - rhs.m32,
769                self.m33 - rhs.m33,
770            )
771        }
772    }
773
774    impl Mul<f64> for Matrix4 {
775        type Output = Self;
776
777        fn mul(self, rhs: f64) -> Self::Output {
778            Self::new(
779                self.m00 * rhs,
780                self.m01 * rhs,
781                self.m02 * rhs,
782                self.m03 * rhs,
783                self.m10 * rhs,
784                self.m11 * rhs,
785                self.m12 * rhs,
786                self.m13 * rhs,
787                self.m20 * rhs,
788                self.m21 * rhs,
789                self.m22 * rhs,
790                self.m23 * rhs,
791                self.m30 * rhs,
792                self.m31 * rhs,
793                self.m32 * rhs,
794                self.m33 * rhs,
795            )
796        }
797    }
798
799    impl Div<f64> for Matrix4 {
800        type Output = Self;
801
802        fn div(self, rhs: f64) -> Self::Output {
803            Self::new(
804                self.m00 / rhs,
805                self.m01 / rhs,
806                self.m02 / rhs,
807                self.m03 / rhs,
808                self.m10 / rhs,
809                self.m11 / rhs,
810                self.m12 / rhs,
811                self.m13 / rhs,
812                self.m20 / rhs,
813                self.m21 / rhs,
814                self.m22 / rhs,
815                self.m23 / rhs,
816                self.m30 / rhs,
817                self.m31 / rhs,
818                self.m32 / rhs,
819                self.m33 / rhs,
820            )
821        }
822    }
823
824    impl Neg for Matrix4 {
825        type Output = Self;
826
827        fn neg(self) -> Self::Output {
828            Self::new(
829                -self.m00, -self.m01, -self.m02, -self.m03, -self.m10, -self.m11, -self.m12,
830                -self.m13, -self.m20, -self.m21, -self.m22, -self.m23, -self.m30, -self.m31,
831                -self.m32, -self.m33,
832            )
833        }
834    }
835
836    impl Mul for Matrix4 {
837        type Output = Self;
838
839        fn mul(self, rhs: Self) -> Self::Output {
840            Self::new(
841                dot4(
842                    self.m00, rhs.m00, self.m01, rhs.m10, self.m02, rhs.m20, self.m03, rhs.m30,
843                ),
844                dot4(
845                    self.m00, rhs.m01, self.m01, rhs.m11, self.m02, rhs.m21, self.m03, rhs.m31,
846                ),
847                dot4(
848                    self.m00, rhs.m02, self.m01, rhs.m12, self.m02, rhs.m22, self.m03, rhs.m32,
849                ),
850                dot4(
851                    self.m00, rhs.m03, self.m01, rhs.m13, self.m02, rhs.m23, self.m03, rhs.m33,
852                ),
853                dot4(
854                    self.m10, rhs.m00, self.m11, rhs.m10, self.m12, rhs.m20, self.m13, rhs.m30,
855                ),
856                dot4(
857                    self.m10, rhs.m01, self.m11, rhs.m11, self.m12, rhs.m21, self.m13, rhs.m31,
858                ),
859                dot4(
860                    self.m10, rhs.m02, self.m11, rhs.m12, self.m12, rhs.m22, self.m13, rhs.m32,
861                ),
862                dot4(
863                    self.m10, rhs.m03, self.m11, rhs.m13, self.m12, rhs.m23, self.m13, rhs.m33,
864                ),
865                dot4(
866                    self.m20, rhs.m00, self.m21, rhs.m10, self.m22, rhs.m20, self.m23, rhs.m30,
867                ),
868                dot4(
869                    self.m20, rhs.m01, self.m21, rhs.m11, self.m22, rhs.m21, self.m23, rhs.m31,
870                ),
871                dot4(
872                    self.m20, rhs.m02, self.m21, rhs.m12, self.m22, rhs.m22, self.m23, rhs.m32,
873                ),
874                dot4(
875                    self.m20, rhs.m03, self.m21, rhs.m13, self.m22, rhs.m23, self.m23, rhs.m33,
876                ),
877                dot4(
878                    self.m30, rhs.m00, self.m31, rhs.m10, self.m32, rhs.m20, self.m33, rhs.m30,
879                ),
880                dot4(
881                    self.m30, rhs.m01, self.m31, rhs.m11, self.m32, rhs.m21, self.m33, rhs.m31,
882                ),
883                dot4(
884                    self.m30, rhs.m02, self.m31, rhs.m12, self.m32, rhs.m22, self.m33, rhs.m32,
885                ),
886                dot4(
887                    self.m30, rhs.m03, self.m31, rhs.m13, self.m32, rhs.m23, self.m33, rhs.m33,
888                ),
889            )
890        }
891    }
892
893    impl Mul<Vector4> for Matrix4 {
894        type Output = Vector4;
895
896        fn mul(self, rhs: Vector4) -> Self::Output {
897            let x_component = rhs.x();
898            let y_component = rhs.y();
899            let z_component = rhs.z();
900            let w_component = rhs.w();
901
902            Vector4::new(
903                dot4(
904                    self.m00,
905                    x_component,
906                    self.m01,
907                    y_component,
908                    self.m02,
909                    z_component,
910                    self.m03,
911                    w_component,
912                ),
913                dot4(
914                    self.m10,
915                    x_component,
916                    self.m11,
917                    y_component,
918                    self.m12,
919                    z_component,
920                    self.m13,
921                    w_component,
922                ),
923                dot4(
924                    self.m20,
925                    x_component,
926                    self.m21,
927                    y_component,
928                    self.m22,
929                    z_component,
930                    self.m23,
931                    w_component,
932                ),
933                dot4(
934                    self.m30,
935                    x_component,
936                    self.m31,
937                    y_component,
938                    self.m32,
939                    z_component,
940                    self.m33,
941                    w_component,
942                ),
943            )
944        }
945    }
946
947    const fn determinant3(rows: [[f64; 3]; 3]) -> f64 {
948        (rows[0][0] * ((rows[1][1] * rows[2][2]) - (rows[1][2] * rows[2][1])))
949            - (rows[0][1] * ((rows[1][0] * rows[2][2]) - (rows[1][2] * rows[2][0])))
950            + (rows[0][2] * ((rows[1][0] * rows[2][1]) - (rows[1][1] * rows[2][0])))
951    }
952}
953
954#[cfg(test)]
955mod tests {
956    use super::{Matrix2, Matrix3, Matrix4};
957    use use_vector::{Vector2, Vector3, Vector4};
958
959    fn assert_close(left: f64, right: f64) {
960        assert!((left - right).abs() < 1.0e-10, "left={left}, right={right}");
961    }
962
963    fn assert_matrix2_close(left: Matrix2, right: Matrix2) {
964        assert_close(left.m00, right.m00);
965        assert_close(left.m01, right.m01);
966        assert_close(left.m10, right.m10);
967        assert_close(left.m11, right.m11);
968    }
969
970    fn assert_matrix3_close(left: Matrix3, right: Matrix3) {
971        assert_close(left.m00, right.m00);
972        assert_close(left.m01, right.m01);
973        assert_close(left.m02, right.m02);
974        assert_close(left.m10, right.m10);
975        assert_close(left.m11, right.m11);
976        assert_close(left.m12, right.m12);
977        assert_close(left.m20, right.m20);
978        assert_close(left.m21, right.m21);
979        assert_close(left.m22, right.m22);
980    }
981
982    fn assert_matrix4_close(left: Matrix4, right: Matrix4) {
983        assert_close(left.m00, right.m00);
984        assert_close(left.m01, right.m01);
985        assert_close(left.m02, right.m02);
986        assert_close(left.m03, right.m03);
987        assert_close(left.m10, right.m10);
988        assert_close(left.m11, right.m11);
989        assert_close(left.m12, right.m12);
990        assert_close(left.m13, right.m13);
991        assert_close(left.m20, right.m20);
992        assert_close(left.m21, right.m21);
993        assert_close(left.m22, right.m22);
994        assert_close(left.m23, right.m23);
995        assert_close(left.m30, right.m30);
996        assert_close(left.m31, right.m31);
997        assert_close(left.m32, right.m32);
998        assert_close(left.m33, right.m33);
999    }
1000
1001    #[test]
1002    fn matrix2_construction_and_constants_work() {
1003        let matrix = Matrix2::new(1.0, 2.0, 3.0, 4.0);
1004
1005        assert_eq!(Matrix2::ZERO, Matrix2::new(0.0, 0.0, 0.0, 0.0));
1006        assert_eq!(Matrix2::IDENTITY, Matrix2::new(1.0, 0.0, 0.0, 1.0));
1007        assert_eq!(matrix, Matrix2::from_rows([[1.0, 2.0], [3.0, 4.0]]));
1008        assert_eq!(matrix, Matrix2::from_cols([[1.0, 3.0], [2.0, 4.0]]));
1009    }
1010
1011    #[test]
1012    fn matrix2_arithmetic_and_products_work() {
1013        let left = Matrix2::new(1.0, 2.0, 3.0, 4.0);
1014        let right = Matrix2::new(5.0, 6.0, 7.0, 8.0);
1015        let scale = Matrix2::from_rows([[2.0, 0.0], [0.0, 3.0]]);
1016
1017        assert_eq!(left + right, Matrix2::new(6.0, 8.0, 10.0, 12.0));
1018        assert_eq!(left - right, Matrix2::new(-4.0, -4.0, -4.0, -4.0));
1019        assert_eq!(left * 2.0, Matrix2::new(2.0, 4.0, 6.0, 8.0));
1020        assert_eq!(left / 2.0, Matrix2::new(0.5, 1.0, 1.5, 2.0));
1021        assert_eq!(-left, Matrix2::new(-1.0, -2.0, -3.0, -4.0));
1022        assert_eq!(left * scale, Matrix2::new(2.0, 6.0, 6.0, 12.0));
1023        assert_eq!(left * Vector2::new(1.0, -1.0), Vector2::new(-1.0, -1.0));
1024        assert_eq!(Matrix2::IDENTITY * left, left);
1025        assert_eq!(left * Matrix2::IDENTITY, left);
1026    }
1027
1028    #[test]
1029    fn matrix2_transpose_determinant_trace_and_inverse_work() {
1030        let matrix = Matrix2::new(4.0, 7.0, 2.0, 6.0);
1031
1032        assert_eq!(matrix.transpose(), Matrix2::new(4.0, 2.0, 7.0, 6.0));
1033        assert_eq!(matrix.transpose().transpose(), matrix);
1034        assert_close(matrix.determinant(), 10.0);
1035        assert_close(matrix.trace(), 10.0);
1036
1037        let inverse = matrix.inverse().expect("matrix should invert");
1038
1039        assert_matrix2_close(inverse, Matrix2::new(0.6, -0.7, -0.2, 0.4));
1040        assert_matrix2_close(matrix * inverse, Matrix2::IDENTITY);
1041        assert!(Matrix2::new(1.0, 2.0, 2.0, 4.0).inverse().is_none());
1042        assert!(
1043            Matrix2::new(f64::INFINITY, 0.0, 0.0, 1.0)
1044                .inverse()
1045                .is_none()
1046        );
1047    }
1048
1049    #[test]
1050    fn matrix3_construction_and_constants_work() {
1051        let matrix = Matrix3::new(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0);
1052
1053        assert_eq!(
1054            Matrix3::ZERO,
1055            Matrix3::new(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
1056        );
1057        assert_eq!(
1058            Matrix3::IDENTITY,
1059            Matrix3::new(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
1060        );
1061        assert_eq!(
1062            matrix,
1063            Matrix3::from_rows([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0],])
1064        );
1065        assert_eq!(
1066            matrix,
1067            Matrix3::from_cols([[1.0, 4.0, 7.0], [2.0, 5.0, 8.0], [3.0, 6.0, 9.0],])
1068        );
1069    }
1070
1071    #[test]
1072    fn matrix3_arithmetic_and_products_work() {
1073        let left = Matrix3::new(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 10.0);
1074        let right = Matrix3::new(9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0);
1075        let scale = Matrix3::from_rows([[2.0, 0.0, 0.0], [0.0, 3.0, 0.0], [0.0, 0.0, 4.0]]);
1076
1077        assert_eq!(
1078            left + right,
1079            Matrix3::new(10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 11.0)
1080        );
1081        assert_eq!(
1082            left - right,
1083            Matrix3::new(-8.0, -6.0, -4.0, -2.0, 0.0, 2.0, 4.0, 6.0, 9.0)
1084        );
1085        assert_eq!(
1086            left * 2.0,
1087            Matrix3::new(2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 20.0)
1088        );
1089        assert_eq!(
1090            left / 2.0,
1091            Matrix3::new(0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0)
1092        );
1093        assert_eq!(
1094            -left,
1095            Matrix3::new(-1.0, -2.0, -3.0, -4.0, -5.0, -6.0, -7.0, -8.0, -10.0)
1096        );
1097        assert_eq!(
1098            left * scale,
1099            Matrix3::new(2.0, 6.0, 12.0, 8.0, 15.0, 24.0, 14.0, 24.0, 40.0)
1100        );
1101        assert_eq!(
1102            left * Vector3::new(1.0, 0.0, -1.0),
1103            Vector3::new(-2.0, -2.0, -3.0)
1104        );
1105        assert_eq!(Matrix3::IDENTITY * left, left);
1106        assert_eq!(left * Matrix3::IDENTITY, left);
1107    }
1108
1109    #[test]
1110    fn matrix3_transpose_determinant_trace_and_inverse_work() {
1111        let matrix = Matrix3::new(1.0, 2.0, 3.0, 0.0, 1.0, 4.0, 5.0, 6.0, 0.0);
1112
1113        assert_eq!(
1114            matrix.transpose(),
1115            Matrix3::new(1.0, 0.0, 5.0, 2.0, 1.0, 6.0, 3.0, 4.0, 0.0)
1116        );
1117        assert_eq!(matrix.transpose().transpose(), matrix);
1118        assert_close(matrix.determinant(), 1.0);
1119        assert_close(matrix.trace(), 2.0);
1120
1121        let inverse = matrix.inverse().expect("matrix should invert");
1122
1123        assert_eq!(
1124            inverse,
1125            Matrix3::new(-24.0, 18.0, 5.0, 20.0, -15.0, -4.0, -5.0, 4.0, 1.0)
1126        );
1127        assert_matrix3_close(matrix * inverse, Matrix3::IDENTITY);
1128        assert!(
1129            Matrix3::new(1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 7.0, 8.0, 9.0)
1130                .inverse()
1131                .is_none()
1132        );
1133        assert!(
1134            Matrix3::new(f64::INFINITY, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
1135                .inverse()
1136                .is_none()
1137        );
1138    }
1139
1140    #[test]
1141    fn matrix4_construction_and_constants_work() {
1142        let matrix = Matrix4::from_rows([
1143            [1.0, 2.0, 3.0, 4.0],
1144            [5.0, 6.0, 7.0, 8.0],
1145            [2.0, 0.0, 1.0, 3.0],
1146            [4.0, 1.0, 0.0, 2.0],
1147        ]);
1148
1149        assert_eq!(
1150            Matrix4::ZERO,
1151            Matrix4::new(
1152                0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1153            )
1154        );
1155        assert_eq!(
1156            Matrix4::IDENTITY,
1157            Matrix4::new(
1158                1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0,
1159            )
1160        );
1161        assert_eq!(
1162            matrix,
1163            Matrix4::from_cols([
1164                [1.0, 5.0, 2.0, 4.0],
1165                [2.0, 6.0, 0.0, 1.0],
1166                [3.0, 7.0, 1.0, 0.0],
1167                [4.0, 8.0, 3.0, 2.0],
1168            ])
1169        );
1170    }
1171
1172    #[test]
1173    fn matrix4_arithmetic_and_products_work() {
1174        let left = Matrix4::from_rows([
1175            [1.0, 2.0, 3.0, 4.0],
1176            [5.0, 6.0, 7.0, 8.0],
1177            [2.0, 0.0, 1.0, 3.0],
1178            [4.0, 1.0, 0.0, 2.0],
1179        ]);
1180        let right = Matrix4::from_rows([
1181            [4.0, 3.0, 2.0, 1.0],
1182            [0.0, 1.0, 2.0, 3.0],
1183            [1.0, 0.0, 1.0, 0.0],
1184            [2.0, 1.0, 0.0, 2.0],
1185        ]);
1186        let scale = Matrix4::from_rows([
1187            [2.0, 0.0, 0.0, 0.0],
1188            [0.0, 3.0, 0.0, 0.0],
1189            [0.0, 0.0, 4.0, 0.0],
1190            [0.0, 0.0, 0.0, 5.0],
1191        ]);
1192
1193        assert_eq!(
1194            left + right,
1195            Matrix4::from_rows([
1196                [5.0, 5.0, 5.0, 5.0],
1197                [5.0, 7.0, 9.0, 11.0],
1198                [3.0, 0.0, 2.0, 3.0],
1199                [6.0, 2.0, 0.0, 4.0],
1200            ])
1201        );
1202        assert_eq!(
1203            left - right,
1204            Matrix4::from_rows([
1205                [-3.0, -1.0, 1.0, 3.0],
1206                [5.0, 5.0, 5.0, 5.0],
1207                [1.0, 0.0, 0.0, 3.0],
1208                [2.0, 0.0, 0.0, 0.0],
1209            ])
1210        );
1211        assert_eq!(
1212            left * 2.0,
1213            Matrix4::from_rows([
1214                [2.0, 4.0, 6.0, 8.0],
1215                [10.0, 12.0, 14.0, 16.0],
1216                [4.0, 0.0, 2.0, 6.0],
1217                [8.0, 2.0, 0.0, 4.0],
1218            ])
1219        );
1220        assert_eq!(
1221            left / 2.0,
1222            Matrix4::from_rows([
1223                [0.5, 1.0, 1.5, 2.0],
1224                [2.5, 3.0, 3.5, 4.0],
1225                [1.0, 0.0, 0.5, 1.5],
1226                [2.0, 0.5, 0.0, 1.0],
1227            ])
1228        );
1229        assert_eq!(
1230            -left,
1231            Matrix4::from_rows([
1232                [-1.0, -2.0, -3.0, -4.0],
1233                [-5.0, -6.0, -7.0, -8.0],
1234                [-2.0, 0.0, -1.0, -3.0],
1235                [-4.0, -1.0, 0.0, -2.0],
1236            ])
1237        );
1238        assert_eq!(
1239            left * scale,
1240            Matrix4::from_rows([
1241                [2.0, 6.0, 12.0, 20.0],
1242                [10.0, 18.0, 28.0, 40.0],
1243                [4.0, 0.0, 4.0, 15.0],
1244                [8.0, 3.0, 0.0, 10.0],
1245            ])
1246        );
1247        assert_eq!(
1248            left * Vector4::new(1.0, 0.0, -1.0, 2.0),
1249            Vector4::new(6.0, 14.0, 7.0, 8.0)
1250        );
1251        assert_eq!(Matrix4::IDENTITY * left, left);
1252        assert_eq!(left * Matrix4::IDENTITY, left);
1253    }
1254
1255    #[test]
1256    fn matrix4_transpose_determinant_and_trace_work() {
1257        let matrix = Matrix4::from_rows([
1258            [1.0, 2.0, 3.0, 4.0],
1259            [5.0, 6.0, 7.0, 8.0],
1260            [2.0, 0.0, 1.0, 3.0],
1261            [4.0, 1.0, 0.0, 2.0],
1262        ]);
1263        let triangular = Matrix4::from_rows([
1264            [2.0, 1.0, 0.0, 0.0],
1265            [0.0, 3.0, 4.0, 0.0],
1266            [0.0, 0.0, 5.0, 6.0],
1267            [0.0, 0.0, 0.0, 7.0],
1268        ]);
1269
1270        assert_eq!(
1271            matrix.transpose(),
1272            Matrix4::from_rows([
1273                [1.0, 5.0, 2.0, 4.0],
1274                [2.0, 6.0, 0.0, 1.0],
1275                [3.0, 7.0, 1.0, 0.0],
1276                [4.0, 8.0, 3.0, 2.0],
1277            ])
1278        );
1279        assert_eq!(matrix.transpose().transpose(), matrix);
1280        assert_close(matrix.trace(), 10.0);
1281        assert_close(triangular.determinant(), 210.0);
1282        assert_matrix4_close(Matrix4::IDENTITY * triangular, triangular);
1283    }
1284}
use_matrix/lib.rs

use_matrix/
lib.rs
