use std::ops; fn det(m: &[f64], x: usize) -> f64 { if x == 1 { m[0] } else if x == 2 { m[0] * m[3] - m[1] * m[2] } else { let mut d = 0.0; for n in 0..x { let ms = m .iter() .enumerate() .skip(x) .filter(|(i, _)| (i % x) != n) .map(|(_, v)| *v) .collect::>(); d += (-1.0f64).powi(n as i32) * m[n] * det(&ms, x - 1); } d } } #[derive(Clone, Debug)] pub struct Matrix { data: Box<[f64]>, height: usize, width: usize, } impl Matrix { pub fn new(height: usize, width: usize) -> Matrix { Matrix { data: vec![0.0; height * width].into_boxed_slice(), height, width, } } pub fn transpose(&self) -> Matrix { let mut matrix = Matrix::new(self.width, self.height); for c in 0..self.width { for r in 0..self.height { matrix[(c, r)] = self[(r, c)]; } } matrix } pub fn minor(&self, row_n: usize, col_n: usize) -> Matrix { let mut matrix = Matrix::new(self.height - 1, self.width - 1); let mut nr = 0; for r in 0..self.height { if r == row_n { continue; } let mut nc = 0; for c in 0..self.width { if c == col_n { continue; } matrix[(nr, nc)] = self[(r, c)]; nc += 1; } nr += 1; } matrix } pub fn determinant(&self) -> f64 { debug_assert!(self.width == self.height); det(&self.data, self.width) } pub fn adjugate(&self) -> Matrix { debug_assert!(self.width == self.height); let mut matrix = Matrix::new(self.height, self.width); if matrix.height == 2 { matrix[(0, 0)] = self[(1, 1)]; matrix[(0, 1)] = -self[(0, 1)]; matrix[(1, 0)] = -self[(1, 0)]; matrix[(1, 1)] = self[(0, 0)]; } else { for r in 0..matrix.height { for c in 0..matrix.width { let sign = if (r + c) % 2 == 0 { 1.0 } else { -1.0 }; matrix[(r, c)] = self.minor(r, c).determinant() * sign; } } } matrix } pub fn inverse(&self) -> Matrix { let mut matrix = Matrix::new(self.width, self.height); if self.height == self.width && self.height == 1 { matrix[(0, 0)] = 1.0 / self[(0, 0)]; } else { } matrix } } impl ops::Index<(usize, usize)> for Matrix { type Output = f64; fn index(&self, pos: (usize, usize)) -> &Self::Output { &self.data[(self.width * pos.0) + pos.1] } } impl ops::IndexMut<(usize, usize)> for Matrix { fn index_mut(&mut self, pos: (usize, usize)) -> &mut Self::Output { &mut self.data[(self.width * pos.0) + pos.1] } } impl<'a> ops::Mul<&'a Matrix> for f64 { type Output = Matrix; fn mul(self, rhs: &'a Matrix) -> Matrix { let mut matrix = Matrix::new(rhs.height, rhs.width); for r in 0..rhs.height { for c in 0..rhs.width { matrix[(r, c)] = self * rhs[(r, c)]; } } matrix } } impl<'a> ops::Mul<&'a Matrix> for Matrix { type Output = Matrix; fn mul(self, rhs: &'a Matrix) -> Matrix { let mut matrix = Matrix::new(self.height, rhs.width); for r in 0..matrix.height { for c in 0..matrix.width { let mut value = 0.0; for x in 0..self.width { value += self[(r, x)] * rhs[(x, c)]; } matrix[(r, c)] = value; } } matrix } } impl<'a> ops::Mul<&'a Matrix> for &'a Matrix { type Output = Matrix; fn mul(self, rhs: &'a Matrix) -> Matrix { let mut matrix = Matrix::new(self.height, rhs.width); for r in 0..matrix.height { for c in 0..matrix.width { let mut value = 0.0; for x in 0..self.width { value += self[(r, x)] * rhs[(x, c)]; } matrix[(r, c)] = value; } } matrix } } impl<'a> ops::Add<&'a Matrix> for &'a Matrix { type Output = Matrix; fn add(self, rhs: &'a Matrix) -> Matrix { let mut matrix = Matrix::new(self.height, self.width); for r in 0..matrix.height { for c in 0..matrix.width { matrix[(r, c)] = self[(r, c)] + rhs[(r, c)]; } } matrix } }