feo_math/linear_algebra/
matrix3.rs

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