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

use_matrix/
lib.rs
