|
@@ -1,7 +1,9 @@
|
|
|
use std::ops::{Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign};
|
|
|
use std::ptr::swap;
|
|
|
|
|
|
-use crate::{Float, Numeric};
|
|
|
+use crate::{Float, Numeric, Primitive};
|
|
|
+#[cfg(feature = "angular")]
|
|
|
+use crate::Angular;
|
|
|
|
|
|
/// Struct representing a dense matrix.
|
|
|
#[repr(transparent)]
|
|
@@ -159,9 +161,11 @@ impl<T: Numeric + Float + PartialOrd, const DIMENSION: usize> SquareMatrix<T, DI
|
|
|
///
|
|
|
/// This function is unsafe, because even when the result is false,
|
|
|
/// the matrix may already have been modified.
|
|
|
- pub unsafe fn lup_decompose(&mut self, pivot: &mut RowVector<usize, DIMENSION>, tolerance: T) -> bool {
|
|
|
- for i in 0..DIMENSION {
|
|
|
- pivot[(0, i)] = i;
|
|
|
+ pub unsafe fn lup_decompose(&mut self, pivot: &mut Vec<usize>, tolerance: T) -> bool {
|
|
|
+ pivot.clear();
|
|
|
+ pivot.reserve(DIMENSION + 1);
|
|
|
+ for i in 0..(DIMENSION + 1) {
|
|
|
+ pivot.push(i);
|
|
|
}
|
|
|
|
|
|
let zero = T::whole(0);
|
|
@@ -183,13 +187,13 @@ impl<T: Numeric + Float + PartialOrd, const DIMENSION: usize> SquareMatrix<T, DI
|
|
|
}
|
|
|
|
|
|
if max_i != i {
|
|
|
- let p_i = (&mut pivot.data[0][i]) as *mut usize;
|
|
|
- let p_i_max = (&mut pivot.data[0][max_i]) as *mut usize;
|
|
|
- swap(p_i, p_i_max);
|
|
|
+ pivot.swap(i, max_i);
|
|
|
|
|
|
let p_a: *mut [T; DIMENSION] = (&mut self.data[i]) as *mut [T; DIMENSION];
|
|
|
let p_a_max: *mut [T; DIMENSION] = (&mut self.data[max_i]) as *mut [T; DIMENSION];
|
|
|
swap(p_a, p_a_max);
|
|
|
+
|
|
|
+ pivot[DIMENSION] += 1;
|
|
|
}
|
|
|
|
|
|
for j in (i + 1)..DIMENSION {
|
|
@@ -207,9 +211,28 @@ impl<T: Numeric + Float + PartialOrd, const DIMENSION: usize> SquareMatrix<T, DI
|
|
|
true
|
|
|
}
|
|
|
|
|
|
+ /// Calculate the determinant of the input matrix.
|
|
|
+ pub fn determinant(mut self, tolerance: T) -> Option<T> {
|
|
|
+ let mut pivot = Vec::new();
|
|
|
+ if !unsafe { self.lup_decompose(&mut pivot, tolerance) } {
|
|
|
+ return None;
|
|
|
+ }
|
|
|
+
|
|
|
+ let mut result = self[(0, 0)];
|
|
|
+ for i in 1..DIMENSION {
|
|
|
+ result *= self[(i, i)];
|
|
|
+ }
|
|
|
+
|
|
|
+ if (pivot[DIMENSION] - DIMENSION) % 2 == 0 {
|
|
|
+ Some(result)
|
|
|
+ } else {
|
|
|
+ Some(-result)
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
/// Calculate the inverse of the input matrix.
|
|
|
pub fn inverted(mut self, tolerance: T) -> Option<Self> {
|
|
|
- let mut pivot: RowVector<usize, DIMENSION> = RowVector::new();
|
|
|
+ let mut pivot = Vec::new();
|
|
|
if !unsafe { self.lup_decompose(&mut pivot, tolerance) } {
|
|
|
return None;
|
|
|
}
|
|
@@ -221,7 +244,7 @@ impl<T: Numeric + Float + PartialOrd, const DIMENSION: usize> SquareMatrix<T, DI
|
|
|
|
|
|
for j in 0..DIMENSION {
|
|
|
for i in 0..DIMENSION {
|
|
|
- result[(i, j)] = if pivot[(0, i)] == j { one } else { zero };
|
|
|
+ result[(i, j)] = if pivot[i] == j { one } else { zero };
|
|
|
|
|
|
for k in 0..i {
|
|
|
let factor = result[(k, j)];
|
|
@@ -373,6 +396,61 @@ impl<T: Numeric, const ROWS: usize, const COMMON: usize, const COLUMNS: usize> M
|
|
|
}
|
|
|
}
|
|
|
|
|
|
+/// Apply a transformation matrix to a column vector.
|
|
|
+impl<T: Numeric + Float + Primitive> Mul<ColumnVector<T, 3>> for SquareMatrix<T, 4> {
|
|
|
+ type Output = ColumnVector<T, 3>;
|
|
|
+
|
|
|
+ fn mul(self, rhs: ColumnVector<T, 3>) -> Self::Output {
|
|
|
+ let x_old = rhs[(0, 0)];
|
|
|
+ let y_old = rhs[(1, 0)];
|
|
|
+ let z_old = rhs[(2, 0)];
|
|
|
+
|
|
|
+ cvec![
|
|
|
+ self[(0,0)] * x_old + self[(0, 1)] * y_old + self[(0, 2)] * z_old + self[(0, 3)],
|
|
|
+ self[(1,0)] * x_old + self[(1, 1)] * y_old + self[(1, 2)] * z_old + self[(1, 3)],
|
|
|
+ self[(2,0)] * x_old + self[(2, 1)] * y_old + self[(2, 2)] * z_old + self[(2, 3)],
|
|
|
+ ]
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+impl<T: Numeric + Float + Primitive> SquareMatrix<T, 4> {
|
|
|
+ /// Create a 3D translation matrix.
|
|
|
+ pub fn translation(x: T, y: T, z: T) -> Self {
|
|
|
+ let mut result = SquareMatrix::identity();
|
|
|
+ result[(0, 3)] = x;
|
|
|
+ result[(1, 3)] = y;
|
|
|
+ result[(2, 3)] = z;
|
|
|
+ result
|
|
|
+ }
|
|
|
+
|
|
|
+ /// Create a 3D scale matrix.
|
|
|
+ pub fn scale(x: T, y: T, z: T) -> Self {
|
|
|
+ let one = unsafe { T::whole(1) };
|
|
|
+ SquareMatrix::diagonal(&[x, y, z, one])
|
|
|
+ }
|
|
|
+
|
|
|
+ /// Create a 3D rotation matrix.
|
|
|
+ #[cfg(feature = "angular")]
|
|
|
+ pub fn rotation<A: Angular<T>>(angle: A, x: T, y: T, z: T) -> Self {
|
|
|
+ let zero = unsafe { T::whole(0) };
|
|
|
+ let one = unsafe { T::whole(1) };
|
|
|
+
|
|
|
+ let rad_angle = angle.radiant();
|
|
|
+ let c = rad_angle.cos();
|
|
|
+ let s = rad_angle.sin();
|
|
|
+ let omc = one - c;
|
|
|
+
|
|
|
+ Self {
|
|
|
+ data: [
|
|
|
+ [x * x * omc + c, x * y * omc - z * s, x * z * omc + y * s, zero],
|
|
|
+ [y * x * omc + z * s, y * y * omc + c, y * z * omc - x * s, zero],
|
|
|
+ [x * z * omc - y * s, y * z * omc + x * s, z * z * omc + c, zero],
|
|
|
+ [zero, zero, zero, one],
|
|
|
+ ],
|
|
|
+ }
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
#[cfg(test)]
|
|
|
mod tests {
|
|
|
use crate::{ColumnVector, Complex, cplx, Float, GenericMatrix, Numeric, RowVector, SquareMatrix};
|