feo_math/linear_algebra/
matrix4.rs

1use crate::{F32Fmt, SignOps, /*Sum,*/ Two};
2
3use {
4    super::{
5        SqMatrix,
6        vector4::Vector4,
7    },
8    crate::{
9        Construct,
10        One, Zero
11    },
12    std::{
13        fmt,
14        ops::{Add, Mul, Sub, Div, Rem, Neg}
15    },
16};
17
18/// A 4x4 matrix with row major order.
19#[derive(PartialEq, Clone, Copy)]
20pub struct Matrix4<T>{pub m:[[T; 4]; 4]}
21
22impl<T> Matrix4<T> {
23    /// A cleaner way of making a new Matrix4.
24    /// Returns a new Matrix4<T>
25    /// # Arguments
26    /// * `row0` - A 4 element array representing the first row of the matrix
27    /// * `row2` - A 4 element array representing the second row of the matrix
28    /// * `row3` - A 4 element array representing the third row of the matrix
29    /// * `row4` - A 4 element array representing the fourth row of the matrix
30    pub fn new(
31            row0: [T; 4], 
32            row1: [T; 4], 
33            row2: [T; 4], 
34            row3: [T; 4]) -> Self {
35        Matrix4 {
36            m: [
37                row0,
38                row1,
39                row2,
40                row3
41            ]
42        }
43    }
44    
45    #[allow(clippy::too_many_arguments)]
46    pub fn new_val(
47        r0c0: T, r0c1: T, r0c2: T, r0c3: T,
48        r1c0: T, r1c1: T, r1c2: T, r1c3: T,
49        r2c0: T, r2c1: T, r2c2: T, r2c3: T,
50        r3c0: T, r3c1: T, r3c2: T, r3c3: T) -> Self {
51        Matrix4 {
52            m: [
53                [r0c0, r0c1, r0c2, r0c3],
54                [r1c0, r1c1, r1c2, r1c3],
55                [r2c0, r2c1, r2c2, r2c3],
56                [r3c0, r3c1, r3c2, r3c3]
57            ]
58        }
59    }
60
61    /// one: T equivalent of one zero: T equivalent of zero
62    pub fn identity() -> Self
63    where T: Zero + One {
64        Self::ONE
65    }
66
67    #[allow(dead_code)]
68    pub fn determinant(&self) -> T 
69    where T: Add<T, Output = T> + Mul<T, Output = T> + Sub<T, Output = T> + Copy {
70        self.m[0][0] * (
71                (self.m[1][1] * self.m[2][2] * self.m[3][3]) +
72                (self.m[1][2] * self.m[2][3] * self.m[3][1]) +
73                (self.m[1][3] * self.m[2][1] * self.m[3][2])
74                - (self.m[1][3] * self.m[2][2] * self.m[3][1])
75                - (self.m[1][2] * self.m[2][1] * self.m[3][3])
76                - (self.m[1][1] * self.m[2][3] * self.m[3][2])
77            ) - self.m[1][0] * (
78                (self.m[0][1] * self.m[2][2] * self.m[3][3]) +
79                (self.m[0][2] * self.m[2][3] * self.m[3][1]) +
80                (self.m[0][3] * self.m[2][1] * self.m[3][2])
81                - (self.m[0][3] * self.m[2][2] * self.m[3][1])
82                - (self.m[0][2] * self.m[2][1] * self.m[3][3])
83                - (self.m[0][1] * self.m[2][3] * self.m[3][2])
84            ) + self.m[2][0] * (
85                (self.m[0][1] * self.m[1][2] * self.m[3][3]) +
86                (self.m[0][2] * self.m[1][3] * self.m[3][1]) +
87                (self.m[0][3] * self.m[1][1] * self.m[3][2])
88                - (self.m[0][3] * self.m[1][2] * self.m[3][1])
89                - (self.m[0][2] * self.m[1][1] * self.m[3][3])
90                - (self.m[0][1] * self.m[1][3] * self.m[3][2])
91            ) - self.m[3][0] * (
92                (self.m[0][1] * self.m[1][2] * self.m[2][3]) +
93                (self.m[0][2] * self.m[1][3] * self.m[2][1]) +
94                (self.m[0][3] * self.m[1][1] * self.m[2][2])
95                - (self.m[0][3] * self.m[1][2] * self.m[2][1])
96                - (self.m[0][2] * self.m[1][1] * self.m[2][3])
97                - (self.m[0][1] * self.m[1][3] * self.m[2][2])
98            )
99    }
100
101    // transpose
102    // matrix (a, b, c, d)
103    //        (e, f, g, h)
104    //        (i, j, k, l)
105    //        (m, n, o, p)
106    //
107    // matrix (a, d, g)
108    //        (b, e, h)
109    //        (c, f, i) // nd
110    pub fn transpose(&self) -> Matrix4<T>
111    where T: Copy {
112        Matrix4 {
113            m: [
114                [self.m[0][0], self.m[1][0], self.m[2][0], self.m[3][0]],
115                [self.m[0][1], self.m[1][1], self.m[2][1], self.m[3][1]],
116                [self.m[0][2], self.m[1][2], self.m[2][2], self.m[3][2]],
117                [self.m[0][3], self.m[1][3], self.m[2][3], self.m[3][3]]
118            ]
119        }
120    }
121
122    #[allow(clippy::possible_missing_comma)]
123    pub fn cofactor(self) -> Self
124    where T: Mul<T, Output = T> + Neg<Output = T> + Add<T, Output = T> + Sub<T, Output = T> + Copy {
125        Matrix4 {
126            m: [
127                [
128                    (self.m[1][1] * self.m[2][2] * self.m[3][3]) +
129                    (self.m[2][1] * self.m[3][2] * self.m[1][3]) +
130                    (self.m[3][1] * self.m[1][2] * self.m[2][3])
131                    - (self.m[3][1] * self.m[2][2] * self.m[1][3])
132                    - (self.m[2][1] * self.m[1][2] * self.m[3][3])
133                    - (self.m[1][1] * self.m[3][2] * self.m[2][3]),
134
135                    - (self.m[1][0] * self.m[2][2] * self.m[3][3])
136                    - (self.m[2][0] * self.m[3][2] * self.m[1][3])
137                    - (self.m[3][0] * self.m[1][2] * self.m[2][3]) + 
138                    (self.m[3][0] * self.m[2][2] * self.m[1][3]) +
139                    (self.m[2][0] * self.m[1][2] * self.m[3][3]) +
140                    (self.m[1][0] * self.m[3][2] * self.m[2][3]),
141
142                    (self.m[1][0] * self.m[2][1] * self.m[3][3]) + 
143                    (self.m[2][0] * self.m[3][1] * self.m[1][3]) + 
144                    (self.m[3][0] * self.m[1][1] * self.m[2][3])
145                    - (self.m[3][0] * self.m[2][1] * self.m[1][3]) 
146                    - (self.m[2][0] * self.m[1][1] * self.m[3][3])
147                    - (self.m[1][0] * self.m[3][1] * self.m[2][3]),
148
149                    - (self.m[1][0] * self.m[2][1] * self.m[3][2])
150                    - (self.m[2][0] * self.m[3][1] * self.m[1][2])
151                    - (self.m[3][0] * self.m[1][1] * self.m[2][2]) + 
152                    (self.m[3][0] * self.m[2][1] * self.m[1][2]) + 
153                    (self.m[2][0] * self.m[1][1] * self.m[3][2]) + 
154                    (self.m[1][0] * self.m[3][1] * self.m[2][2])
155                ], 
156                [
157                    - (self.m[0][1] * self.m[2][2] * self.m[3][3])
158                    - (self.m[2][1] * self.m[3][2] * self.m[0][3])
159                    - (self.m[3][1] * self.m[0][2] * self.m[2][3]) +
160                    (self.m[3][1] * self.m[2][2] * self.m[0][3]) +
161                    (self.m[2][1] * self.m[0][2] * self.m[3][3]) +
162                    (self.m[0][1] * self.m[3][2] * self.m[2][3]),
163
164                    (self.m[0][0] * self.m[2][2] * self.m[3][3]) +
165                    (self.m[2][0] * self.m[3][2] * self.m[0][3]) +
166                    (self.m[3][0] * self.m[0][2] * self.m[2][3])
167                    - (self.m[3][0] * self.m[2][2] * self.m[0][3])
168                    - (self.m[2][0] * self.m[0][2] * self.m[3][3])
169                    - (self.m[0][0] * self.m[3][2] * self.m[2][3]),
170
171                    - (self.m[0][0] * self.m[2][1] * self.m[3][3]) // ERROR
172                    - (self.m[2][0] * self.m[3][1] * self.m[0][3])
173                    - (self.m[3][0] * self.m[0][1] * self.m[2][3]) +
174                    (self.m[3][0] * self.m[2][1] * self.m[0][3]) +
175                    (self.m[2][0] * self.m[0][1] * self.m[3][3]) +
176                    (self.m[0][0] * self.m[3][1] * self.m[2][3]), // FIX: self.m[1][3] to self.m[2][3] worked
177
178                    (self.m[0][0] * self.m[2][1] * self.m[3][2]) +
179                    (self.m[2][0] * self.m[3][1] * self.m[0][2]) +
180                    (self.m[3][0] * self.m[0][1] * self.m[2][2])
181                    - (self.m[3][0] * self.m[2][1] * self.m[0][2])
182                    - (self.m[2][0] * self.m[0][1] * self.m[3][2])
183                    - (self.m[0][0] * self.m[3][1] * self.m[2][2])
184                ], 
185                [
186                    (self.m[0][1] * self.m[1][2] * self.m[3][3]) +
187                    (self.m[1][1] * self.m[3][2] * self.m[0][3]) +
188                    (self.m[3][1] * self.m[0][2] * self.m[1][3])
189                    - (self.m[3][1] * self.m[1][2] * self.m[0][3])
190                    - (self.m[1][1] * self.m[0][2] * self.m[3][3])
191                    - (self.m[0][1] * self.m[3][2] * self.m[1][3]),
192
193                    - (self.m[0][0] * self.m[1][2] * self.m[3][3])
194                    - (self.m[1][0] * self.m[3][2] * self.m[0][3])
195                    - (self.m[3][0] * self.m[0][2] * self.m[1][3]) +
196                    (self.m[3][0] * self.m[1][2] * self.m[0][3]) +
197                    (self.m[1][0] * self.m[0][2] * self.m[3][3]) +
198                    (self.m[0][0] * self.m[3][2] * self.m[1][3]),
199
200                    (self.m[0][0] * self.m[1][1] * self.m[3][3]) + 
201                    (self.m[1][0] * self.m[3][1] * self.m[0][3]) + 
202                    (self.m[3][0] * self.m[0][1] * self.m[1][3])
203                    - (self.m[3][0] * self.m[1][1] * self.m[0][3]) 
204                    - (self.m[1][0] * self.m[0][1] * self.m[3][3])
205                    - (self.m[0][0] * self.m[3][1] * self.m[1][3]),
206
207                    - (self.m[0][0] * self.m[1][1] * self.m[3][2])
208                    - (self.m[1][0] * self.m[3][1] * self.m[0][2])
209                    - (self.m[3][0] * self.m[0][1] * self.m[1][2]) + 
210                    (self.m[3][0] * self.m[1][1] * self.m[0][2]) + 
211                    (self.m[1][0] * self.m[0][1] * self.m[3][2]) + 
212                    (self.m[0][0] * self.m[3][1] * self.m[1][2])
213                ], 
214                [
215                    - (self.m[0][1] * self.m[1][2] * self.m[2][3])
216                    - (self.m[1][1] * self.m[2][2] * self.m[0][3])
217                    - (self.m[2][1] * self.m[0][2] * self.m[1][3]) +
218                    (self.m[2][1] * self.m[1][2] * self.m[0][3]) +
219                    (self.m[1][1] * self.m[0][2] * self.m[2][3]) +
220                    (self.m[0][1] * self.m[2][2] * self.m[1][3]),
221
222                    (self.m[0][0] * self.m[1][2] * self.m[2][3]) +
223                    (self.m[1][0] * self.m[2][2] * self.m[0][3]) +
224                    (self.m[2][0] * self.m[0][2] * self.m[1][3])
225                    - (self.m[2][0] * self.m[1][2] * self.m[0][3])
226                    - (self.m[1][0] * self.m[0][2] * self.m[2][3])
227                    - (self.m[0][0] * self.m[2][2] * self.m[1][3]),
228
229                    - (self.m[0][0] * self.m[1][1] * self.m[2][3])
230                    - (self.m[1][0] * self.m[2][1] * self.m[0][3])
231                    - (self.m[2][0] * self.m[0][1] * self.m[1][3]) +
232                    (self.m[2][0] * self.m[1][1] * self.m[0][3]) +
233                    (self.m[1][0] * self.m[0][1] * self.m[2][3]) +
234                    (self.m[0][0] * self.m[2][1] * self.m[1][3]),
235
236                    (self.m[0][0] * self.m[1][1] * self.m[2][2]) +
237                    (self.m[1][0] * self.m[2][1] * self.m[0][2]) +
238                    (self.m[2][0] * self.m[0][1] * self.m[1][2]) 
239                    - (self.m[2][0] * self.m[1][1] * self.m[0][2])
240                    - (self.m[1][0] * self.m[0][1] * self.m[2][2])
241                    - (self.m[0][0] * self.m[2][1] * self.m[1][2])
242                ]
243            ],
244        }
245    }
246    
247    /// Get the adjugate matrix of a 3x3 matrix.
248    pub fn adjugate(self) -> Self
249    where T: Mul<T, Output = T> + Neg<Output = T> + Add<T, Output = T> + Sub<T, Output = T> + Copy {
250        // The adjugate matrix is equal to cofactor after it is transposed.
251        self.cofactor().transpose()
252    }
253
254    /// Get the inverse of a 4x4 matrix.
255    ///
256    /// # Arguments
257    /// * `self` - The matrix the function is being called for
258    ///
259    /// Inserts a value computed from `f` into the option if it is [`None`], then
260    /// returns a mutable reference to the contained value.
261    ///
262    /// # Examples
263    /// 
264    /// ```rust
265    /// use feo_math::linear_algebra::matrix4::Matrix4;
266    /// let x = Matrix4::new(
267    ///   [1.0,12.0,13.0,14.0],
268    ///   [21.0,2.0,23.0,24.0],
269    ///   [31.0,32.0,3.0,34.0],
270    ///   [41.0,42.0,43.0,4.0]);
271    /// let x_inverse = x.inverse();
272    /// assert_eq!(x_inverse, Matrix4::new(
273    ///   [-57.0/932.0,17.0/932.0,53.0/4660.0,37.0/4660.0],
274    ///   [139.0/4660.0,-8.0/233.0,51.0/4660.0,2.0/233.0],
275    ///   [25.0/932.0,69.0/4660.0,-21.0/932.0,41.0/4660.0],
276    ///   [59.0/2330.0,67.0/4660.0,5.0/466.0,-15.0/932.0])); // TODO is this transposed should not be !!!!!!!!!!!!!!!
277    /// ```
278    ///
279    pub fn inverse(&self) -> Matrix4<T>
280    where T: Add<T, Output = T> + Mul<T, Output = T> + Sub<T, Output = T> + Div<T, Output = T> + Neg<Output = T> + Copy { // TODO result
281        let det: T = self.determinant();
282
283        // Adjugate matrix
284        let adjugate_mat = self.adjugate();
285
286        adjugate_mat / det
287    }
288}
289
290impl<T> Construct<T> for Matrix4<T> where T: Construct<T> {}
291impl<T> SqMatrix<T, Vector4<T>> for Matrix4<T> where T: Construct<T> {}
292
293impl<T> fmt::Debug for Matrix4<T>
294where T: fmt::Debug + Copy {
295    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
296        write!(f,"\n[ {:?}\t {:?}\t {:?}\t {:?}\t ]
297                  \n[ {:?}\t {:?}\t {:?}\t {:?}\t ]
298                  \n[ {:?}\t {:?}\t {:?}\t {:?}\t ]
299                  \n[ {:?}\t {:?}\t {:?}\t {:?}\t ]",
300                  self.m[0][0], self.m[0][1], self.m[0][2], self.m[0][3],
301                  self.m[1][0], self.m[1][1], self.m[1][2], self.m[1][3],
302                  self.m[2][0], self.m[2][1], self.m[2][2], self.m[2][3],
303                  self.m[3][0], self.m[3][1], self.m[3][2], self.m[3][3])
304    }
305}
306
307impl<T> Add for Matrix4<T>
308where T: Add<T, Output = T> + Copy {
309    type Output = Self;
310
311    /// Matrix4<T> + T
312    /// = Matrix4<T>
313    /// ```rust
314    /// use feo_math::linear_algebra::matrix4::Matrix4;
315    /// let mat0 = Matrix4::new(
316    ///     [1, 2, 3, 4],
317    ///     [3, 4, 1, 2],
318    ///     [4, 1, 2, 3],
319    ///     [2, 3, 4, 1]
320    /// );
321    /// let mat1 = Matrix4::new(
322    ///     [1, 2, 3, 4],
323    ///     [3, 4, 2, 2],
324    ///     [4, 1, 2, 3],
325    ///     [2, 3, 4, 1]
326    /// );
327    /// let expected = Matrix4::new(
328    ///     [2, 4, 6, 8],
329    ///     [6, 8, 3, 4],
330    ///     [8, 2, 4, 6],
331    ///     [4, 6, 8, 2]
332    /// );
333    /// assert_eq!(mat0 + mat1, expected);
334    /// ```
335    fn add(self, rhs: Self) -> Self::Output {
336        Matrix4{
337            m: [
338                [self.m[0][0] + rhs.m[0][0], self.m[0][1] + rhs.m[0][1], self.m[0][2] + rhs.m[0][2], self.m[0][3] + rhs.m[0][3]],
339                [self.m[1][0] + rhs.m[1][0], self.m[1][1] + rhs.m[1][1], self.m[1][2] + rhs.m[1][2], self.m[1][3] + rhs.m[1][3]],
340                [self.m[2][0] + rhs.m[2][0], self.m[2][1] + rhs.m[2][1], self.m[2][2] + rhs.m[2][2], self.m[2][3] + rhs.m[2][3]],
341                [self.m[3][0] + rhs.m[3][0], self.m[3][1] + rhs.m[3][1], self.m[3][2] + rhs.m[3][2], self.m[3][3] + rhs.m[3][3]],
342            ]
343        }
344    }
345}
346
347impl<T> Add<T> for Matrix4<T>
348where T: Add<T, Output = T> + Copy {
349    type Output = Self;
350
351    /// Matrix4<T> + T
352    /// = Matrix4<T>
353    /// ```rust
354    /// use feo_math::linear_algebra::matrix4::Matrix4;
355    /// let mat = Matrix4::new(
356    ///     [1, 2, 3, 4],
357    ///     [3, 4, 1, 2],
358    ///     [4, 1, 2, 3],
359    ///     [2, 3, 4, 1]
360    /// );
361    /// let val = 2;
362    /// let expected = Matrix4::new(
363    ///     [3, 4, 5, 6],
364    ///     [5, 6, 3, 4],
365    ///     [6, 3, 4, 5],
366    ///     [4, 5, 6, 3]
367    /// );
368    /// assert_eq!(mat + val, expected);
369    /// ```
370    fn add(self, rhs: T) -> Self::Output {
371        Matrix4 {
372            m: [
373                [self.m[0][0] + rhs, self.m[0][1] + rhs, self.m[0][2] + rhs, self.m[0][3] + rhs],
374                [self.m[1][0] + rhs, self.m[1][1] + rhs, self.m[1][2] + rhs, self.m[1][3] + rhs],
375                [self.m[2][0] + rhs, self.m[2][1] + rhs, self.m[2][2] + rhs, self.m[2][3] + rhs],
376                [self.m[3][0] + rhs, self.m[3][1] + rhs, self.m[3][2] + rhs, self.m[3][3] + rhs],
377            ]
378        }
379    }
380}
381
382impl<T> Sub for Matrix4<T>
383where T: Sub<T, Output = T> + Copy {
384    type Output = Self;
385
386    /// Matrix4<T> - T
387    /// = Matrix4<T>
388    /// ```rust
389    /// use feo_math::linear_algebra::matrix4::Matrix4;
390    /// let mat0 = Matrix4::new(
391    ///     [1, 2, 3, 4],
392    ///     [3, 4, 1, 2],
393    ///     [4, 1, 2, 3],
394    ///     [2, 3, 4, 1]
395    /// );
396    /// let mat1 = Matrix4::new(
397    ///     [1, 2, 3, 4],
398    ///     [3, 4, 2, 2],
399    ///     [4, 1, 2, 3],
400    ///     [2, 3, 4, 1]
401    /// );
402    /// let expected = Matrix4::new(
403    ///     [0, 0,  0, 0],
404    ///     [0, 0, -1, 0],
405    ///     [0, 0,  0, 0],
406    ///     [0, 0,  0, 0]
407    /// );
408    /// assert_eq!(mat0 - mat1, expected);
409    /// ```
410    fn sub(self, rhs: Self) -> Self::Output {
411        Matrix4 {
412            m: [
413                [self.m[0][0] - rhs.m[0][0], self.m[0][1] - rhs.m[0][1], self.m[0][2] - rhs.m[0][2], self.m[0][3] - rhs.m[0][3]],
414                [self.m[1][0] - rhs.m[1][0], self.m[1][1] - rhs.m[1][1], self.m[1][2] - rhs.m[1][2], self.m[1][3] - rhs.m[1][3]],
415                [self.m[2][0] - rhs.m[2][0], self.m[2][1] - rhs.m[2][1], self.m[2][2] - rhs.m[2][2], self.m[2][3] - rhs.m[2][3]],
416                [self.m[3][0] - rhs.m[3][0], self.m[3][1] - rhs.m[3][1], self.m[3][2] - rhs.m[3][2], self.m[3][3] - rhs.m[3][3]],
417            ]
418        }
419    }
420}
421
422impl<T> Sub<T> for Matrix4<T>
423where T: Sub<T, Output = T> + Copy {
424    type Output = Self;
425
426    /// Matrix4<T> - T
427    /// = Matrix4<T>
428    /// ```rust
429    /// use feo_math::linear_algebra::matrix4::Matrix4;
430    /// let mat = Matrix4::new(
431    ///     [1, 2, 3, 4],
432    ///     [3, 4, 1, 2],
433    ///     [4, 1, 2, 3],
434    ///     [2, 3, 4, 1]
435    /// );
436    /// let val = 2;
437    /// let expected = Matrix4::new(
438    ///     [-1,  0,  1,  2],
439    ///     [ 1,  2, -1,  0],
440    ///     [ 2, -1,  0,  1],
441    ///     [ 0,  1,  2, -1]
442    /// );
443    /// assert_eq!(mat - val, expected);
444    /// ```
445    fn sub(self, rhs: T) -> Self::Output {
446        Matrix4 {
447            m: [
448                [self.m[0][0] - rhs, self.m[0][1] - rhs, self.m[0][2] - rhs, self.m[0][3] - rhs],
449                [self.m[1][0] - rhs, self.m[1][1] - rhs, self.m[1][2] - rhs, self.m[1][3] - rhs],
450                [self.m[2][0] - rhs, self.m[2][1] - rhs, self.m[2][2] - rhs, self.m[2][3] - rhs],
451                [self.m[3][0] - rhs, self.m[3][1] - rhs, self.m[3][2] - rhs, self.m[3][3] - rhs],
452            ]
453        }
454    }
455}
456
457impl<T> Mul for Matrix4<T> where T: Add<T, Output = T> + Mul<T, Output = T> + Copy{
458    type Output = Self;
459
460    /// Matrix4 * Matrix4 = Matrix4
461    /// ```rust
462    /// use feo_math::linear_algebra::matrix4::Matrix4;
463    /// let mat0 = Matrix4::new(
464    ///     [4, 6, 6, 0],
465    ///     [6, 4, 5, 0],
466    ///     [5, 6, 4, 0],
467    ///     [0, 0, 0, 1] 
468    /// );
469    /// let mat1 = Matrix4::new(
470    ///     [1, 2, 3, 0],
471    ///     [2, 3, 1, 0],
472    ///     [2, 1, 2, 0],
473    ///     [0, 0, 0, 1] 
474    /// );
475    /// let expected = Matrix4::new(
476    ///     [28, 32, 30, 0],
477    ///     [24, 29, 32, 0],
478    ///     [25, 32, 29, 0],
479    ///     [0, 0, 0, 1] 
480    /// );
481    /// assert_eq!(mat0 * mat1, expected);
482    /// ```
483    fn mul(self, rhs: Self) -> Self::Output {
484        Matrix4 {
485            m: [
486                [
487                    self.m[0][0] * rhs.m[0][0]
488                        + self.m[0][1] * rhs.m[1][0]
489                        + self.m[0][2] * rhs.m[2][0]
490                        + self.m[0][3] * rhs.m[3][0],
491                    self.m[0][0] * rhs.m[0][1]
492                        + self.m[0][1] * rhs.m[1][1]
493                        + self.m[0][2] * rhs.m[2][1]
494                        + self.m[0][3] * rhs.m[3][1],
495                    self.m[0][0] * rhs.m[0][2]
496                        + self.m[0][1] * rhs.m[1][2]
497                        + self.m[0][2] * rhs.m[2][2]
498                        + self.m[0][3] * rhs.m[3][2],
499                    self.m[0][0] * rhs.m[0][3]
500                        + self.m[0][1] * rhs.m[1][3]
501                        + self.m[0][2] * rhs.m[2][3]
502                        + self.m[0][3] * rhs.m[3][3],
503                ],
504                [
505                    self.m[1][0] * rhs.m[0][0]
506                        + self.m[1][1] * rhs.m[1][0]
507                        + self.m[1][2] * rhs.m[2][0]
508                        + self.m[1][3] * rhs.m[3][0],
509                    self.m[1][0] * rhs.m[0][1]
510                        + self.m[1][1] * rhs.m[1][1]
511                        + self.m[1][2] * rhs.m[2][1]
512                        + self.m[1][3] * rhs.m[3][1],
513                    self.m[1][0] * rhs.m[0][2]
514                        + self.m[1][1] * rhs.m[1][2]
515                        + self.m[1][2] * rhs.m[2][2]
516                        + self.m[1][3] * rhs.m[2][3],
517                    self.m[1][0] * rhs.m[0][3]
518                        + self.m[1][1] * rhs.m[1][3]
519                        + self.m[1][2] * rhs.m[2][3]
520                        + self.m[1][3] * rhs.m[3][3],
521                ],
522                [
523                    self.m[2][0] * rhs.m[0][0]
524                        + self.m[2][1] * rhs.m[1][0]
525                        + self.m[2][2] * rhs.m[2][0]
526                        + self.m[2][3] * rhs.m[3][0],
527                    self.m[2][0] * rhs.m[0][1]
528                        + self.m[2][1] * rhs.m[1][1]
529                        + self.m[2][2] * rhs.m[2][1]
530                        + self.m[2][3] * rhs.m[3][1],
531                    self.m[2][0] * rhs.m[0][2]
532                        + self.m[2][1] * rhs.m[1][2]
533                        + self.m[2][2] * rhs.m[2][2]
534                        + self.m[2][3] * rhs.m[3][2],
535                    self.m[2][0] * rhs.m[0][3]
536                        + self.m[2][1] * rhs.m[1][3]
537                        + self.m[2][2] * rhs.m[2][3]
538                        + self.m[2][3] * rhs.m[3][3],
539                ],
540                [
541                    self.m[3][0] * rhs.m[0][0]
542                        + self.m[3][1] * rhs.m[1][0]
543                        + self.m[3][2] * rhs.m[2][0]
544                        + self.m[3][3] * rhs.m[3][0],
545                    self.m[3][0] * rhs.m[0][1]
546                        + self.m[3][1] * rhs.m[1][1]
547                        + self.m[3][2] * rhs.m[2][1]
548                        + self.m[3][3] * rhs.m[3][1],
549                    self.m[3][0] * rhs.m[0][2]
550                        + self.m[3][1] * rhs.m[1][2]
551                        + self.m[3][2] * rhs.m[2][2]
552                        + self.m[3][3] * rhs.m[3][2],
553                    self.m[3][0] * rhs.m[0][3]
554                        + self.m[3][1] * rhs.m[1][3]
555                        + self.m[3][2] * rhs.m[2][3]
556                        + self.m[3][3] * rhs.m[3][3],
557                ],
558            ],
559        }
560    }
561}
562
563impl<T> Mul<Vector4<T>> for Matrix4<T> where T: Add<T, Output = T> + Mul<T, Output = T> + Copy{
564    type Output = Vector4<T>;
565
566    /// Matrix4 * Vec3d = type Vec3d
567    /// # Examples
568    /// ```rust
569    /// use feo_math::linear_algebra::vector4::Vector4;
570    /// use feo_math::linear_algebra::matrix4::Matrix4;
571    /// let mat = Matrix4::new(
572    ///     [1, 2, 3, 4],
573    ///     [3, 4, 1, 2],
574    ///     [4, 1, 2, 3],
575    ///     [2, 3, 4, 1]
576    /// );
577    /// let vec = Vector4::new(1, 4, 2, 3);
578    /// assert_eq!(mat * vec, Vector4::new(27, 27, 21, 25));
579    /// ```
580    fn mul(self, rhs: Vector4<T>) -> Self::Output {
581        Vector4::new(
582            rhs.0 * self.m[0][0] + rhs.1 * self.m[0][1] + rhs.2 * self.m[0][2] + rhs.3 * self.m[0][3],
583            rhs.0 * self.m[1][0] + rhs.1 * self.m[1][1] + rhs.2 * self.m[1][2] + rhs.3 * self.m[1][3],
584            rhs.0 * self.m[2][0] + rhs.1 * self.m[2][1] + rhs.2 * self.m[2][2] + rhs.3 * self.m[2][3],
585            rhs.0 * self.m[3][0] + rhs.1 * self.m[3][1] + rhs.2 * self.m[3][2] + rhs.3 * self.m[3][3],
586        )
587    }
588}
589
590impl<T> Mul<T> for Matrix4<T> where T: Mul<T, Output = T> + Copy{
591    type Output = Matrix4<T>;
592
593    /// Matrix4<T> * T
594    /// = Matrix4<T>
595    /// ```rust
596    /// use feo_math::linear_algebra::matrix4::Matrix4;
597    /// let mat = Matrix4::new(
598    ///     [1, 2, 3, 4],
599    ///     [3, 4, 1, 2],
600    ///     [4, 1, 2, 3],
601    ///     [2, 3, 4, 1]
602    /// );
603    /// let val = 2;
604    /// let expected = Matrix4::new(
605    ///     [2, 4, 6, 8],
606    ///     [6, 8, 2, 4],
607    ///     [8, 2, 4, 6],
608    ///     [4, 6, 8, 2]
609    /// );
610    /// assert_eq!(mat * val, expected);
611    /// ```
612    fn mul(self, rhs: T) -> Self::Output {
613        Matrix4 {
614            m: [
615                [self.m[0][0] * rhs, self.m[0][1] * rhs, self.m[0][2] * rhs, self.m[0][3] * rhs],
616                [self.m[1][0] * rhs, self.m[1][1] * rhs, self.m[1][2] * rhs, self.m[1][3] * rhs],
617                [self.m[2][0] * rhs, self.m[2][1] * rhs, self.m[2][2] * rhs, self.m[2][3] * rhs],
618                [self.m[3][0] * rhs, self.m[3][1] * rhs, self.m[3][2] * rhs, self.m[3][3] * rhs],
619            ],
620        }
621    }
622}
623
624impl<T> Div for Matrix4<T>
625where T: Add<T, Output = T> + Sub<T, Output = T> + Mul<T, Output = T> + Div<T, Output = T> + Neg<Output = T> + Copy {
626    type Output = Self;
627
628    /// # Examples
629    /// ```rust
630    /// use feo_math::linear_algebra::vector4::Vector4;
631    /// use feo_math::linear_algebra::matrix4::Matrix4;
632    /// let mat0 = Matrix4::new(
633    ///     [4, 5, 6, 0],
634    ///     [6, 4, 5, 0],
635    ///     [5, 6, 4, 0],
636    ///     [0, 0, 0, 1],
637    /// );
638    /// let mat1 = Matrix4::new(
639    ///     [1, 2, 3, 0],
640    ///     [2, 3, 1, 0],
641    ///     [3, 1, 2, 0],
642    ///     [0, 0, 0, 1]
643    /// );
644    /// let expected = Matrix4::new(
645    ///     [1, 0, 0, 0],
646    ///     [0, 0, 1, 0],
647    ///     [0, 1, 0, 0],
648    ///     [0, 0, 0, 1],
649    /// );
650    /// assert_eq!(mat0 / mat1, expected);
651    /// ```
652    fn div(self, rhs: Self) -> Self::Output {
653        // self * rhs.inverse()
654
655        let determinant = rhs.determinant();
656
657        let adjugate_mat = rhs.adjugate();
658
659        // first multiply to mitigate error with non Floatingpoint types
660        self * adjugate_mat / determinant
661    }
662}
663
664impl<T> Div<T> for Matrix4<T> where T: Div<T, Output = T> + Copy {
665    type Output = Matrix4<T>;
666
667    /// Matrix4<T> / T
668    /// = Matrix4<T>
669    /// ```rust
670    /// use feo_math::linear_algebra::matrix4::Matrix4;
671    /// let mat = Matrix4::new(
672    ///     [1, 2, 3, 4],
673    ///     [3, 4, 1, 2],
674    ///     [4, 1, 2, 3],
675    ///     [2, 3, 4, 1]
676    /// );
677    /// let val = 2;
678    /// let expected = Matrix4::new(
679    ///     [0, 1, 1, 2],
680    ///     [1, 2, 0, 1],
681    ///     [2, 0, 1, 1],
682    ///     [1, 1, 2, 0]
683    /// );
684    /// assert_eq!(mat / val, expected);
685    /// ```
686    fn div(self, rhs: T) -> Self::Output {
687        Matrix4 {
688            m: [
689                [self.m[0][0] / rhs, self.m[0][1] / rhs, self.m[0][2] / rhs, self.m[0][3] / rhs],
690                [self.m[1][0] / rhs, self.m[1][1] / rhs, self.m[1][2] / rhs, self.m[1][3] / rhs],
691                [self.m[2][0] / rhs, self.m[2][1] / rhs, self.m[2][2] / rhs, self.m[2][3] / rhs],
692                [self.m[3][0] / rhs, self.m[3][1] / rhs, self.m[3][2] / rhs, self.m[3][3] / rhs],
693            ],
694        }
695    }
696}
697
698
699impl<T> Rem for Matrix4<T>
700where T: Rem<T, Output = T> + Copy {
701    type Output = Self;
702
703    /// Matrix4<T> % Matrix4<T>
704    /// = Matrix4<T>
705    /// Finds the remainder after an ELEMENT WISE division
706    /// # Examples
707    /// ```rust
708    /// use feo_math::linear_algebra::matrix4::Matrix4;
709    /// let mat0 = Matrix4::new(
710    ///     [4, 5, 6, 3],
711    ///     [6, 4, 5, 2],
712    ///     [5, 6, 4, 1],
713    ///     [8, 1, 3, 2],
714    /// );
715    /// let mat1 = Matrix4::new(
716    ///     [1, 2, 3, 4],
717    ///     [2, 3, 1, 5],
718    ///     [3, 1, 2, 6],
719    ///     [3, 9, 8, 7],
720    /// );
721    /// let expected = Matrix4::new(
722    ///     [0, 1, 0, 3],
723    ///     [0, 1, 0, 2],
724    ///     [2, 0, 0, 1],
725    ///     [2, 1, 3, 2],
726    /// );
727    /// assert_eq!(mat0 % mat1, expected);
728    /// ```
729    fn rem(self, rhs: Self) -> Self::Output {
730        Matrix4::new(
731            [self.m[0][0] % rhs.m[0][0], self.m[0][1] % rhs.m[0][1], self.m[0][2] % rhs.m[0][2], self.m[0][3] % rhs.m[0][3]],
732            [self.m[1][0] % rhs.m[1][0], self.m[1][1] % rhs.m[1][1], self.m[1][2] % rhs.m[1][2], self.m[1][3] % rhs.m[1][3]],
733            [self.m[2][0] % rhs.m[2][0], self.m[2][1] % rhs.m[2][1], self.m[2][2] % rhs.m[2][2], self.m[2][3] % rhs.m[2][3]],
734            [self.m[3][0] % rhs.m[3][0], self.m[3][1] % rhs.m[3][1], self.m[3][2] % rhs.m[3][2], self.m[3][3] % rhs.m[3][3]],
735        )
736    }
737}
738
739impl<T> Rem<T> for Matrix4<T> where T: Rem<T, Output = T> + Copy {
740    type Output = Matrix4<T>;
741
742    /// Matrix4<T> % T
743    /// = Matrix4<T>
744    /// ```rust
745    /// use feo_math::linear_algebra::matrix4::Matrix4;
746    /// let mat = Matrix4::new(
747    ///     [1, 2, 3, 4],
748    ///     [3, 4, 1, 2],
749    ///     [4, 1, 2, 3],
750    ///     [2, 3, 4, 1]
751    /// );
752    /// let val = 2;
753    /// let expected = Matrix4::new(
754    ///     [1, 0, 1, 0],
755    ///     [1, 0, 1, 0],
756    ///     [0, 1, 0, 1],
757    ///     [0, 1, 0, 1]
758    /// );
759    /// assert_eq!(mat % val, expected);
760    /// ```
761    /// ```should_panic
762    /// use feo_math::linear_algebra::matrix4::Matrix4;
763    /// let mat = Matrix4::new(
764    ///     [1, 2, 3, 4],
765    ///     [3, 4, 1, 2],
766    ///     [4, 1, 2, 3],
767    ///     [2, 3, 4, 1]
768    /// );
769    /// mat % 0;
770    /// ```
771    fn rem(self, rhs: T) -> Self::Output {
772        Matrix4 {
773            m: [
774                [self.m[0][0] % rhs, self.m[0][1] % rhs, self.m[0][2] % rhs, self.m[0][3] % rhs],
775                [self.m[1][0] % rhs, self.m[1][1] % rhs, self.m[1][2] % rhs, self.m[1][3] % rhs],
776                [self.m[2][0] % rhs, self.m[2][1] % rhs, self.m[2][2] % rhs, self.m[2][3] % rhs],
777                [self.m[3][0] % rhs, self.m[3][1] % rhs, self.m[3][2] % rhs, self.m[3][3] % rhs],
778            ],
779        }
780    }
781}
782
783
784impl<T> Neg for Matrix4<T>
785where T: Neg<Output = T> + Copy {
786    type Output = Self;
787    
788    /// -Matrix4<T>
789    /// = Matrix4<T>
790    /// # Examples
791    /// ```rust
792    /// use feo_math::linear_algebra::matrix4::Matrix4;
793    /// let mat = Matrix4::new(
794    ///     [1, 2, 3, 4],
795    ///     [3, 4, 1, 2],
796    ///     [4, 1, 2, 3],
797    ///     [2, 3, 4, 1]
798    /// );
799    /// let expected = Matrix4::new(
800    ///     [-1, -2, -3, -4],
801    ///     [-3, -4, -1, -2],
802    ///     [-4, -1, -2, -3],
803    ///     [-2, -3, -4, -1]
804    /// );
805    /// assert_eq!(-mat, expected);
806    /// ```
807    fn neg(self) -> Self::Output {
808        Matrix4{
809            m: [
810                [-self.m[0][0], -self.m[0][1], -self.m[0][2], -self.m[0][3]],
811                [-self.m[1][0], -self.m[1][1], -self.m[1][2], -self.m[1][3]],
812                [-self.m[2][0], -self.m[2][1], -self.m[2][2], -self.m[2][3]],
813                [-self.m[3][0], -self.m[3][1], -self.m[3][2], -self.m[3][3]],
814            ],
815        }
816    }
817}
818
819impl<T> Zero for Matrix4<T> 
820where T: Zero + Copy {
821    const ZERO: Self = Matrix4 {
822        m: [
823            [T::ZERO, T::ZERO, T::ZERO, T::ZERO],
824            [T::ZERO, T::ZERO, T::ZERO, T::ZERO],
825            [T::ZERO, T::ZERO, T::ZERO, T::ZERO],
826            [T::ZERO, T::ZERO, T::ZERO, T::ZERO]
827        ]
828    };
829}
830
831impl<T> One for Matrix4<T> 
832where T: Zero + One {
833    /// The identity matrix
834    const ONE: Self = Matrix4 {
835        m: [
836            [T::ONE , T::ZERO, T::ZERO, T::ZERO],
837            [T::ZERO, T::ONE , T::ZERO, T::ZERO],
838            [T::ZERO, T::ZERO, T::ONE , T::ZERO],
839            [T::ZERO, T::ZERO, T::ZERO, T::ONE ]
840        ]
841    };
842}
843
844impl<T> Two for Matrix4<T> 
845where T: Zero + Two {
846    /// The identity matrix
847    const TWO: Self = Matrix4 {
848        m: [
849            [T::TWO , T::ZERO, T::ZERO, T::ZERO],
850            [T::ZERO, T::TWO , T::ZERO, T::ZERO],
851            [T::ZERO, T::ZERO, T::TWO , T::ZERO],
852            [T::ZERO, T::ZERO, T::ZERO, T::TWO ]
853        ]
854    };
855}
856
857
858impl<T> From<Matrix4<T>> for [[T; 4]; 4] {
859    fn from(other: Matrix4<T>) -> [[T; 4]; 4] {
860        other.m
861    }
862}
863
864impl<T> F32Fmt for Matrix4<T> 
865where T: F32Fmt + Copy { 
866    type F32Fmt = Matrix4<T::F32Fmt>;
867    #[inline] 
868    fn intoF32Fmt(self) -> Self::F32Fmt {
869        let m = &self.m;
870        Matrix4{
871            m: [
872                [m[0][0].intoF32Fmt(), m[0][1].intoF32Fmt(), m[0][2].intoF32Fmt(), m[0][3].intoF32Fmt()],
873                [m[1][0].intoF32Fmt(), m[1][1].intoF32Fmt(), m[1][2].intoF32Fmt(), m[1][3].intoF32Fmt()],
874                [m[2][0].intoF32Fmt(), m[2][1].intoF32Fmt(), m[2][2].intoF32Fmt(), m[2][3].intoF32Fmt()],
875                [m[3][0].intoF32Fmt(), m[3][1].intoF32Fmt(), m[3][2].intoF32Fmt(), m[3][3].intoF32Fmt()],
876            ]
877        }
878    } 
879    #[inline]
880    fn fromF32Fmt(f32_fmt: Self::F32Fmt) -> Self { 
881        let m = &f32_fmt.m;
882        Matrix4{
883            m: [
884                [T::fromF32Fmt(m[0][0]), T::fromF32Fmt(m[0][1]), T::fromF32Fmt(m[0][2]), T::fromF32Fmt(m[0][2])],
885                [T::fromF32Fmt(m[1][0]), T::fromF32Fmt(m[1][1]), T::fromF32Fmt(m[1][2]), T::fromF32Fmt(m[1][2])],
886                [T::fromF32Fmt(m[2][0]), T::fromF32Fmt(m[2][1]), T::fromF32Fmt(m[2][2]), T::fromF32Fmt(m[2][2])],
887                [T::fromF32Fmt(m[3][0]), T::fromF32Fmt(m[3][1]), T::fromF32Fmt(m[3][2]), T::fromF32Fmt(m[3][2])],
888            ]
889        }
890    }
891
892    fn sqrt(self) -> Self {
893        todo!()
894    }
895
896    fn cbrt(self) -> Self { 
897        todo!() 
898    }
899
900    fn f32_const_mul(self, constant: f32) -> Self {
901        Matrix4{
902            m: [
903                [self.m[0][0].f32_const_mul(constant), self.m[0][1].f32_const_mul(constant), self.m[0][2].f32_const_mul(constant), self.m[0][3].f32_const_mul(constant)],
904                [self.m[1][0].f32_const_mul(constant), self.m[1][1].f32_const_mul(constant), self.m[1][2].f32_const_mul(constant), self.m[1][3].f32_const_mul(constant)],
905                [self.m[2][0].f32_const_mul(constant), self.m[2][1].f32_const_mul(constant), self.m[2][2].f32_const_mul(constant), self.m[2][3].f32_const_mul(constant)],
906                [self.m[3][0].f32_const_mul(constant), self.m[3][1].f32_const_mul(constant), self.m[3][2].f32_const_mul(constant), self.m[3][3].f32_const_mul(constant)],
907            ]
908        }
909    }
910
911    fn sin_mul(self, _mul_by: Self) -> Self where Self: Mul<Self, Output = Self> + Sized {
912        todo!()
913    }
914
915    fn cos_mul(self, _mul_by: Self) -> Self where Self: Mul<Self, Output = Self> + Sized {
916        todo!()
917    }
918
919    fn tan_mul(self, _mul_by: Self) -> Self where Self: Mul<Self, Output = Self> + Sized {
920        todo!()
921    }
922
923    fn asin_mul(self, _mul_by: Self) -> Self where Self: Mul<Self, Output = Self> + Sized {
924        todo!()
925    }
926
927    fn acos_mul(self, _mul_by: Self) -> Self where Self: Mul<Self, Output = Self> + Sized {
928        todo!()
929    }
930
931    fn atan_mul(self, _mul_by: Self) -> Self where Self: Mul<Self, Output = Self> + Sized {
932        todo!()
933    }
934
935    fn atan2_mul(self, _other: Self, _mul_by: Self) -> Self where Self: Mul<Self, Output = Self> + Sized {
936        todo!()
937    }
938
939    fn sinh_mul(self, _mul_by: Self) -> Self where Self: Mul<Self, Output = Self> + Sized {
940        todo!()
941    }
942
943    fn cosh_mul(self, _mul_by: Self) -> Self where Self: Mul<Self, Output = Self> + Sized {
944        todo!()
945    }
946
947    fn tanh_mul(self, _mul_by: Self) -> Self where Self: Mul<Self, Output = Self> + Sized {
948        todo!()
949    }
950}
951
952impl<T> SignOps for Matrix4<T> {
953    fn ptcopysign(self, _sign: Self) -> Self {
954        todo!()
955    }
956
957    fn ptsignum(self) -> i8 {
958        todo!()
959    }
960
961    fn abs(self) -> Self {
962        todo!()
963    }
964}