qm/
m.rs

1//! matrix
2//!
3
4pub mod m3;
5pub mod m4;
6
7use num::Float;
8
9use crate::v::{v3::Vector3, v4::Vector4};
10
11/// cofactor
12pub fn cofactor<F: Float + std::fmt::Debug>(
13  m: &Vec<Vec<F>>, i: usize, j: usize) -> F {
14/*
15  assert_eq!(m.len(), m[0].len());
16  assert!(m.len() >= 2);
17*/
18  if m.len() == 2 {
19    return [[m[1][1]], [-m[1][0]], [-m[0][1]], [m[0][0]]][i][j];
20  }
21  let d = det(
22    &m.iter().enumerate().flat_map(|(ri, r)|
23      if i == ri { vec![] } // to be skipped by flatten
24      else {
25        vec![r.iter().enumerate().flat_map(|(cj, &c)|
26          if j == cj { vec![] } // to be skipped by flatten
27          else { vec![c] }
28        ).collect::<Vec<_>>()]
29      }
30    ).collect::<Vec<_>>()
31  );
32  if (i + j) % 2 == 0 { d } else { -d }
33}
34
35/// transpose
36pub fn transpose<F: Float + std::fmt::Debug>(m: &Vec<Vec<F>>) -> Vec<Vec<F>> {
37  (0..m[0].len()).into_iter().map(|j|
38    (0..m.len()).into_iter().map(|i|
39      m[i][j]
40    ).collect::<Vec<_>>()
41  ).collect::<Vec<_>>()
42}
43
44/// det
45pub fn det<F: Float + std::fmt::Debug>(m: &Vec<Vec<F>>) -> F {
46/*
47  assert_eq!(m.len(), m[0].len());
48  assert!(m.len() >= 2);
49*/
50  if m.len() == 2 { return m[0][0] * m[1][1] - m[0][1] * m[1][0]; } // to fast
51  let mut d = <F>::from(0).unwrap();
52  for (cj, &c) in m[0].iter().enumerate() {
53    d = d + c * cofactor(m, 0, cj); // += ops::AddAssign
54  }
55  d
56}
57
58/// inv
59/// - TODO: recursive cofactor (slow), change to LU decomposition to be faster
60/// - p: prec (assume det = 0)
61pub fn inv<F: Float + std::fmt::Debug>(m: &Vec<Vec<F>>, p: F) ->
62  Option<Vec<Vec<F>>> {
63  assert_eq!(m.len(), m[0].len());
64  assert!(m.len() >= 2);
65  let d = det(m);
66  if crate::prec_eq_f(d, p, <F>::from(0).unwrap()) { return None; }
67  Some(transpose(&m.iter().enumerate().map(|(ri, r)|
68    r.iter().enumerate().map(|(cj, _c)|
69      cofactor(m, ri, cj) / d
70    ).collect::<Vec<_>>()
71  ).collect::<Vec<_>>()))
72}
73
74/// TMatrix
75pub trait TMatrix<F: Float + std::fmt::Debug> {
76  /// constructor col major from v3 (move)
77  fn colmajor3(_m: Vec<Vector3<F>>) -> Self where Self: Sized { panic!("cm3") }
78  /// constructor row major from v3 (move)
79  fn rowmajor3(_m: Vec<Vector3<F>>) -> Self where Self: Sized { panic!("rm3") }
80  /// constructor col major from v4 (move)
81  fn colmajor4(_m: Vec<Vector4<F>>) -> Self where Self: Sized { panic!("cm4") }
82  /// constructor row major from v4 (move)
83  fn rowmajor4(_m: Vec<Vector4<F>>) -> Self where Self: Sized { panic!("rm4") }
84  /// constructor col major
85  fn col_major(m: &Vec<Vec<F>>) -> Self;
86  /// constructor row major
87  fn row_major(m: &Vec<Vec<F>>) -> Self;
88  /// constructor row major
89  fn new(m: &Vec<Vec<F>>) -> Self;
90  /// constructor
91  fn identity() -> Self;
92  /// check equal with precision
93  fn prec_eq(&self, e: F, m: &impl TMatrix<F>) -> bool;
94  /// like as slice v3
95  fn mev3(&self) -> &[Vector3<F>] { panic!("mev3") }
96  /// like as slice v4
97  fn mev4(&self) -> &[Vector4<F>] { panic!("mev4") }
98  /// m dot self
99  fn dot_m(&self, m: &impl TMatrix<F>) -> Self;
100  /// row to v3
101  fn rowv3(&self, _j: usize) -> Vector3<F> { panic!("rowv3") }
102  /// col to v3
103  fn colv3(&self, _i: usize) -> Vector3<F> { panic!("colv3") }
104  /// row to v4
105  fn rowv4(&self, _j: usize) -> Vector4<F> { panic!("rowv4") }
106  /// col to v4
107  fn colv4(&self, _i: usize) -> Vector4<F> { panic!("colv4") }
108  /// to_vec
109  fn to_vec(&self) -> Vec<Vec<F>>;
110  /// transpose
111  fn transpose(&self) -> Self where Self: Sized {
112    Self::col_major(&self.to_vec())
113  }
114  /// det
115  fn det(&self) -> F {
116    crate::m::det(&self.to_vec())
117  }
118  /// inv
119  /// - p: prec (assume det = 0)
120  fn inv(&self, p: F) -> Option<Self> where Self: Sized {
121    match crate::m::inv(&self.to_vec(), p) {
122    None => None,
123    Some(m) => Some(Self::new(&m))
124    }
125  }
126}
qm/m.rs

qm/
m.rs
