devela/num/geom/linear/matrix/methods.rs
1// devela::num::geom::alg::matrix::methods
2//
3//! Methods for matrices.
4//
5// TODO
6// - rank
7// - is_square
8// - is ... (see paper)
9// - impl ops
10
11use crate::{iif, ExtNumConst, Matrix};
12
13impl<const R: usize, const C: usize, const LEN: usize, const MAX_LEN_DET: usize, T: Copy>
14 Matrix<T, R, C, LEN, true, MAX_LEN_DET>
15{
16 /// Returns the transposed matrix (C × R).
17 pub const fn transpose(&self) -> Matrix<T, C, R, LEN, true, MAX_LEN_DET> {
18 let mut transposed_data = [self.data[0]; LEN];
19 let mut i = 0;
20 while i < R {
21 let mut j = 0;
22 while j < C {
23 transposed_data[j * R + i] = self.data[i * C + j];
24 j += 1;
25 }
26 i += 1;
27 }
28 Matrix { data: transposed_data }
29 }
30}
31
32/* numeric methods */
33
34macro_rules! impl_matrix {
35 () => {
36 impl_matrix![int: i8, i16, i32, i64, i128, isize];
37 impl_matrix![float: f32, f64];
38 };
39 (float: $($T:ty),+) => { $(
40 impl_matrix![@both: $T];
41 // impl_matrix![@float: $T];
42 )+ };
43 (int: $($T:ty),+) => { $(
44 impl_matrix![@both: $T];
45 // impl_matrix![@int: $T];
46 )+ };
47
48 (@both: $T:ty) => {
49 // methods for row-major matrices
50 impl<const R: usize, const C: usize, const LEN: usize, const MAX_LEN_DET: usize>
51 Matrix<$T, R, C, LEN, true, MAX_LEN_DET> {
52
53 /* misc */
54
55 /// Returns the identity matrix, or `None` if the matrix is not square.
56 pub const fn identity() -> Option<Self> {
57 if R == C { Some(Self::identity_unchecked()) } else { None }
58 }
59
60 /// Creates an NxN identity matrix.
61 ///
62 /// # Panics
63 /// Panics if the matrix is not square.
64 pub const fn identity_unchecked() -> Self {
65 let data = [<$T>::NUM_ZERO.unwrap(); LEN];
66 let mut matrix = Matrix { data };
67 let mut i = 0;
68 while i < R {
69 *matrix.at_mut(i, i) = <$T>::NUM_ONE;
70 i += 1;
71 }
72 matrix
73 }
74
75 /* ops */
76
77 /// Returns the element-wise sum of two matrices.
78 pub const fn add(&self, other: &Self) -> Self {
79 let mut result = [self.data[0]; LEN];
80 let mut i = 0;
81 while i < LEN {
82 result[i] = self.data[i] + other.data[i];
83 i += 1;
84 }
85 Self { data: result }
86 }
87
88 /// Returns the matrix scaled by a scalar value.
89 pub const fn scale(&self, scalar: $T) -> Self {
90 let mut result = [self.data[0]; LEN];
91 let mut i = 0;
92 while i < LEN {
93 result[i] = self.data[i] * scalar;
94 i += 1;
95 }
96 Self { data: result }
97 }
98
99 /// Computes the matrix product (R × C) * (C × C2) = (R × C2) = LEN2.
100 pub const fn mul<const C2: usize, const LEN2: usize>(
101 &self, other: &Matrix<$T, C, C2, LEN2, true, MAX_LEN_DET>,
102 ) -> Matrix<$T, R, C2, LEN2, true, MAX_LEN_DET> {
103 let mut result = [<$T>::NUM_ZERO.unwrap(); LEN2];
104 let mut i = 0;
105 while i < R {
106 let mut j = 0;
107 while j < C2 {
108 let mut sum = self.data[i * C] * other.data[j]; // 1st element mul
109 let mut k = 1;
110 while k < C {
111 sum += self.data[i * C + k] * other.data[k * C2 + j];
112 k += 1;
113 }
114 result[i * C2 + j] = sum;
115 j += 1;
116 }
117 i += 1;
118 }
119 Matrix { data: result }
120 }
121
122 /* determinant */
123
124 /// Returns the determinant if the matrix is square, or `None` otherwise.
125 ///
126 /// # Panics
127 /// - If called on a non-square matrix in debug mode, it will panic.
128 pub const fn determinant(&self) -> Option<$T> {
129 iif![R == C; Some(self.determinant_unchecked()); None]
130 }
131
132 /// Returns the determinant without checking matrix squareness.
133 ///
134 /// # Panics
135 /// - Panics in debug mode if `R != C` (non-square matrix).
136 /// - May panic on overflow for integer types.
137 pub const fn determinant_unchecked(&self) -> $T {
138 debug_assert![R == C, "a non-square matrix has no determinant"];
139 let s = self;
140 match R {
141 1 => s.at(0, 0),
142 2 => s.at(0, 0) * s.at(1, 1) - s.at(0, 1) * s.at(1, 0),
143 3 => { s.at(0, 0) * (s.at(1, 1) * s.at(2, 2) - s.at(1, 2) * s.at(2, 1))
144 - s.at(0, 1) * (s.at(1, 0) * s.at(2, 2) - s.at(1, 2) * s.at(2, 0))
145 + s.at(0, 2) * (s.at(1, 0) * s.at(2, 1) - s.at(1, 1) * s.at(2, 0)) }
146 _ => {
147 debug_assert![R*C <= MAX_LEN_DET, "MAX_LEN_DET is too small"];
148 let mut buffer = [<$T>::NUM_ZERO.unwrap(); MAX_LEN_DET];
149 Self::submatrix_determinant(R, &s.data, &mut buffer)
150 }
151 }
152 }
153
154 /// Extracts a `(D-1)x(D-1)` submatrix by removing the given row and column.
155 ///
156 /// Arguments
157 /// - `n`: Matrix dimension.
158 /// - `row`, `col`: The row and column to exclude.
159 /// - `matrix`: Source matrix in row-major order.
160 /// - `buffer`: Target buffer to store the submatrix.
161 ///
162 /// # Panics
163 /// Panics in debug mode if `buffer.len() < (n-1)*(n-1)`.
164 pub const fn submatrix(
165 n: usize,
166 row: usize,
167 col: usize,
168 matrix: &[$T],
169 buffer: &mut [$T]
170 ) {
171 debug_assert!(buffer.len() >= (n-1)*(n-1), "Buffer is too small");
172 let (mut idx, mut r) = (0, 0);
173 while r < n {
174 iif![r == row; { r += 1; continue; }];
175 let mut c = 0;
176 while c < n {
177 iif![c == col; { c += 1; continue; }];
178 buffer[idx] = matrix[r * n + c];
179 idx += 1;
180 c += 1;
181 }
182 r += 1;
183 }
184 }
185
186 /// Computes the determinant of a square matrix using Laplace expansion.
187 ///
188 /// This method is **only valid for square matrices** of size `dim × dim`
189 /// and is intended for cases where `dim > 3`.
190 ///
191 /// How it works:
192 /// - Expands along the first row using minors (submatrices).
193 /// - Recursively computes determinants of `(dim-1)×(dim-1)` matrices.
194 ///
195 /// Notes:
196 /// - Uses a temporary `buffer` to avoid repeated allocations.
197 /// - The `buffer` must be at least `(dim-1)^2 + (dim-2)^2` elements long.
198 /// - `MAX_LEN_DET` defines the upper bound for buffer size.
199 /// - It has `O(N!) factorial time complexity due to recursive expansion.
200 ///
201 /// # Panics
202 /// - Panics in debug mode if `buffer.len()` is insufficient.
203 /// - Panics if matrix is not square (should never happen when used internally).
204 pub const fn submatrix_determinant(dim: usize, matrix: &[$T], buffer: &mut [$T]) -> $T {
205 iif![dim == 1; return matrix[0]];
206 iif![dim == 2; return matrix[0] * matrix[3] - matrix[1] * matrix[2]];
207
208 let required_size = (dim - 1) * (dim - 1) + (dim - 2) * (dim - 2);
209 debug_assert!(buffer.len() >= required_size, "buffer is too small");
210
211 let (mut determinant, mut col_idx) = (<$T>::NUM_ZERO.unwrap(), 0);
212 while col_idx < dim {
213 let minor_len = (dim - 1) * (dim - 1);
214 // 1. Split buffer into current level & next (for recursion)
215 let (minor_matrix, next_minor_buffer) = buffer.split_at_mut(minor_len);
216 // 2. Extract the minor submatrix
217 Self::submatrix(dim, 0, col_idx, matrix, minor_matrix);
218 // 3. Compute determinant recursively
219 let sub_det = Self::submatrix_determinant(dim - 1, minor_matrix, next_minor_buffer);
220 determinant += Self::parity_sign(col_idx) * matrix[col_idx] * sub_det;
221 col_idx += 1;
222 }
223 determinant
224 }
225
226 /* utility methods */
227
228 /// Returns an immutable reference to the element at (`row`, `col`).
229 #[inline(always)]
230 pub const fn at(&self, row: usize, col: usize) -> $T {
231 self.data[row * C + col]
232 }
233 /// Returns a shared reference to the element at (`row`, `col`).
234 #[inline(always)]
235 pub const fn at_ref(&self, row: usize, col: usize) -> &$T {
236 &self.data[row * C + col]
237 }
238 /// Returns an exclusive reference to the element at (`row`, `col`).
239 #[inline(always)]
240 pub const fn at_mut(&mut self, row: usize, col: usize) -> &mut $T {
241 &mut self.data[row * C + col]
242 }
243
244 /// Returns alternating ±1 based on the column index for determinant expansion.
245 ///
246 /// Returns `1` for even indices, and `-1` for odd indices.
247 #[inline(always)]
248 const fn parity_sign(i: usize) -> $T {
249 iif![i % 2 == 0; <$T>::NUM_ONE; <$T>::NUM_NEG_ONE.unwrap()]
250 }
251 }
252 };
253 // (@float: $T:ty) => {};
254 // (@int: $T:ty) => {}
255}
256impl_matrix![];