rust_linear_algebra/matrix/
determinant.rs

1use super::Matrix;
2use crate::num::MinusOne;
3use crate::num::One;
4use std::fmt::Debug;
5use std::ops::{Add, Div, Mul, Sub, SubAssign};
6
7impl<K> Matrix<K>
8where
9    K: Debug
10        + Copy
11        + One
12        + MinusOne
13        + Default
14        + PartialEq
15        + Add<Output = K>
16        + Sub<Output = K>
17        + SubAssign<K>
18        + Mul<Output = K>
19        + Div<Output = K>,
20{
21    // ad−bc
22    fn det_matrix_2x2(&self) -> K {
23        self.elements[0][0] * self.elements[1][1] - self.elements[0][1] * self.elements[1][0]
24    }
25
26    // a(ei−fh)−b(di−fg)+c(dh−eg)
27    fn det_matrix_3x3(&self) -> K {
28        self.elements[0][0]
29            * (self.elements[1][1] * self.elements[2][2]
30                - self.elements[1][2] * self.elements[2][1])
31            - self.elements[0][1]
32                * (self.elements[1][0] * self.elements[2][2]
33                    - self.elements[1][2] * self.elements[2][0])
34            + self.elements[0][2]
35                * (self.elements[1][0] * self.elements[2][1]
36                    - self.elements[1][1] * self.elements[2][0])
37    }
38
39    fn det_matrix_4x4(&self) -> K {
40        let a = &self.elements;
41
42        let det_m11 = a[1][1] * (a[2][2] * a[3][3] - a[2][3] * a[3][2])
43            - a[1][2] * (a[2][1] * a[3][3] - a[2][3] * a[3][1])
44            + a[1][3] * (a[2][1] * a[3][2] - a[2][2] * a[3][1]);
45
46        let det_m12 = a[1][0] * (a[2][2] * a[3][3] - a[2][3] * a[3][2])
47            - a[1][2] * (a[2][0] * a[3][3] - a[2][3] * a[3][0])
48            + a[1][3] * (a[2][0] * a[3][2] - a[2][2] * a[3][0]);
49
50        let det_m13 = a[1][0] * (a[2][1] * a[3][3] - a[2][3] * a[3][1])
51            - a[1][1] * (a[2][0] * a[3][3] - a[2][3] * a[3][0])
52            + a[1][3] * (a[2][0] * a[3][1] - a[2][1] * a[3][0]);
53
54        let det_m14 = a[1][0] * (a[2][1] * a[3][2] - a[2][2] * a[3][1])
55            - a[1][1] * (a[2][0] * a[3][2] - a[2][2] * a[3][0])
56            + a[1][2] * (a[2][0] * a[3][1] - a[2][1] * a[3][0]);
57
58        a[0][0] * det_m11 - a[0][1] * det_m12 + a[0][2] * det_m13 - a[0][3] * det_m14
59    }
60
61    fn det_matrix_nxn(&self) -> K {
62        if !self.is_row_echelon_form() {
63            let mut det = K::one();
64            let _ = self.gaussian_elimination(Some(&mut det), None);
65            return det;
66        }
67        let mut det = K::one();
68        for i in 0..self.rows() {
69            det = det * self.elements[i][i];
70        }
71        det
72    }
73
74    pub fn determinant(&self) -> K {
75        if self.cols() != self.rows() {
76            panic!("Matrix must be square");
77        }
78
79        let dim = self.cols();
80
81        match dim {
82            0 => panic!("Empty matrix"),
83            1 => self.elements[0][0],
84            2 => self.det_matrix_2x2(),
85            3 => self.det_matrix_3x3(),
86            4 => self.det_matrix_4x4(),
87            _ => self.det_matrix_nxn(),
88        }
89    }
90}