mini_matrix/matrix.rs
1//! # mini_matrix
2//!
3//! A mini linear algebra library implemented in Rust.
4
5use num::{Float, Num};
6use std::fmt::{Debug, Display};
7
8use std::ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub, SubAssign};
9use std::ops::{Deref, DerefMut, Index, IndexMut};
10
11use crate::Vector;
12
13/// A generic matrix type with `M` rows and `N` columns.
14///
15/// # Examples
16///
17/// ```
18/// use mini_matrix::Matrix;
19///
20/// let matrix = Matrix::<f64, 2, 2>::from([[1.0, 2.0], [3.0, 4.0]]);
21/// assert_eq!(matrix.size(), (2, 2));
22/// ```
23#[derive(Debug, Clone, Copy, PartialEq)]
24pub struct Matrix<T, const M: usize, const N: usize> {
25 /// The underlying storage for the matrix elements.
26 pub store: [[T; N]; M],
27}
28
29impl<T, const M: usize, const N: usize> Matrix<T, M, N>
30where
31 T: Copy + Default,
32{
33 /// Creates a new `Matrix` from the given 2D array.
34 ///
35 /// # Examples
36 ///
37 /// ```
38 /// use mini_matrix::Matrix;
39 ///
40 /// let matrix = Matrix::<i32, 2, 3>::from([[1, 2, 3], [4, 5, 6]]);
41 /// ```
42 pub fn from(data: [[T; N]; M]) -> Self {
43 Self { store: data }
44 }
45
46 /// Returns the dimensions of the matrix as a tuple `(rows, columns)`.
47 ///
48 /// # Examples
49 ///
50 /// ```
51 /// use mini_matrix::Matrix;
52 ///
53 /// let matrix = Matrix::<f64, 3, 4>::zero();
54 /// assert_eq!(matrix.size(), (3, 4));
55 /// ```
56 pub const fn size(&self) -> (usize, usize) {
57 (M, N)
58 }
59
60 /// Creates a new `Matrix` with all elements set to the default value of type `T`.
61 ///
62 /// # Examples
63 ///
64 /// ```
65 /// use mini_matrix::Matrix;
66 ///
67 /// let matrix = Matrix::<f64, 2, 2>::zero();
68 /// assert_eq!(matrix.store, [[0.0, 0.0], [0.0, 0.0]]);
69 /// ```
70 pub fn zero() -> Self {
71 Self {
72 store: [[T::default(); N]; M],
73 }
74 }
75
76 pub fn from_vecs(vecs: Vec<Vec<T>>) -> Self {
77 let mut store = [[T::default(); N]; M];
78 for (i, vec) in vecs.iter().enumerate() {
79 for (j, elem) in vec.iter().enumerate() {
80 store[i][j] = *elem;
81 }
82 }
83 Self { store }
84 }
85
86 #[allow(dead_code)]
87 fn map<F>(&self, mut f: F) -> Matrix<T, M, N>
88 where
89 F: FnMut(T) -> T,
90 {
91 let mut result = Matrix::<T, M, N>::zero();
92 for i in 0..M {
93 for j in 0..N {
94 result[(i, j)] = f(self[(i, j)]);
95 }
96 }
97 result
98 }
99}
100
101impl<T, const M: usize, const N: usize> Matrix<T, M, N>
102where
103 T: AddAssign + SubAssign + MulAssign + Copy + Default,
104{
105 /// Adds another matrix to this matrix in-place.
106 ///
107 /// # Examples
108 ///
109 /// ```
110 /// use mini_matrix::Matrix;
111 ///
112 /// let mut a = Matrix::<i32, 2, 2>::from([[1, 2], [3, 4]]);
113 /// let b = Matrix::<i32, 2, 2>::from([[5, 6], [7, 8]]);
114 /// a.add(&b);
115 /// assert_eq!(a.store, [[6, 8], [10, 12]]);
116 /// ```
117 pub fn add(&mut self, other: &Self) {
118 for (l_row, r_row) in self.store.iter_mut().zip(other.store.iter()) {
119 for (l, r) in l_row.iter_mut().zip(r_row.iter()) {
120 *l += *r;
121 }
122 }
123 }
124
125 /// Subtracts another matrix from this matrix in-place.
126 ///
127 /// # Examples
128 ///
129 /// ```
130 /// use mini_matrix::Matrix;
131 ///
132 /// let mut a = Matrix::<i32, 2, 2>::from([[5, 6], [7, 8]]);
133 /// let b = Matrix::<i32, 2, 2>::from([[1, 2], [3, 4]]);
134 /// a.sub(&b);
135 /// assert_eq!(a.store, [[4, 4], [4, 4]]);
136 /// ```
137 pub fn sub(&mut self, other: &Self) {
138 for (l_row, r_row) in self.store.iter_mut().zip(other.store.iter()) {
139 for (l, r) in l_row.iter_mut().zip(r_row.iter()) {
140 *l -= *r;
141 }
142 }
143 }
144
145 /// Multiplies this matrix by a scalar value in-place.
146 ///
147 /// # Examples
148 ///
149 /// ```
150 /// use mini_matrix::Matrix;
151 ///
152 /// let mut a = Matrix::<i32, 2, 2>::from([[1, 2], [3, 4]]);
153 /// a.scl(2);
154 /// assert_eq!(a.store, [[2, 4], [6, 8]]);
155 /// ```
156 pub fn scl(&mut self, scalar: T) {
157 for row in self.store.iter_mut() {
158 for elem in row.iter_mut() {
159 *elem *= scalar;
160 }
161 }
162 }
163}
164
165impl<T, const M: usize, const N: usize> IndexMut<(usize, usize)> for Matrix<T, M, N> {
166 /// Mutably indexes into the matrix, allowing modification of its elements.
167 ///
168 /// # Arguments
169 ///
170 /// * `index` - A tuple `(row, column)` specifying the element to access.
171 ///
172 /// # Examples
173 ///
174 /// ```
175 /// use mini_matrix::Matrix;
176 ///
177 /// let mut matrix = Matrix::<i32, 2, 2>::from([[1, 2], [3, 4]]);
178 /// matrix[(0, 1)] = 5;
179 /// assert_eq!(matrix.store, [[1, 5], [3, 4]]);
180 /// ```
181 fn index_mut(&mut self, index: (usize, usize)) -> &mut T {
182 &mut self.store[index.0][index.1]
183 }
184}
185
186impl<T, const M: usize, const N: usize> Index<(usize, usize)> for Matrix<T, M, N> {
187 type Output = T;
188
189 /// Immutably indexes into the matrix, allowing read access to its elements.
190 ///
191 /// # Arguments
192 ///
193 /// * `index` - A tuple `(row, column)` specifying the element to access.
194 ///
195 /// # Examples
196 ///
197 /// ```
198 /// use mini_matrix::Matrix;
199 ///
200 /// let matrix = Matrix::<i32, 2, 2>::from([[1, 2], [3, 4]]);
201 /// assert_eq!(matrix[(1, 0)], 3);
202 /// ```
203 fn index(&self, index: (usize, usize)) -> &Self::Output {
204 &self.store.index(index.0).index(index.1)
205 }
206}
207
208impl<T, const M: usize, const N: usize> Deref for Matrix<T, M, N> {
209 type Target = [[T; N]; M];
210
211 /// Dereferences the matrix, allowing it to be treated as a slice.
212 ///
213 /// # Note
214 ///
215 /// This implementation is currently unfinished.
216 ///
217 /// # Examples
218 ///
219 /// ```
220 /// use mini_matrix::Matrix;
221 ///
222 /// let matrix = Matrix::<i32, 2, 2>::from([[1, 2], [3, 4]]);
223 /// // Usage example will be available once implementation is complete
224 /// ```
225 fn deref(&self) -> &Self::Target {
226 &self.store
227 }
228}
229
230impl<T, const M: usize, const N: usize> DerefMut for Matrix<T, M, N> {
231 /// Mutably dereferences the matrix, allowing it to be treated as a mutable slice.
232 ///
233 /// # Note
234 ///
235 /// This implementation is currently unfinished.
236 ///
237 /// # Examples
238 ///
239 /// ```
240 /// use mini_matrix::Matrix;
241 ///
242 /// let mut matrix = Matrix::<i32, 2, 2>::from([[1, 2], [3, 4]]);
243 /// // Usage example will be available once implementation is complete
244 /// ```
245 fn deref_mut(&mut self) -> &mut Self::Target {
246 &mut self.store
247 }
248}
249
250impl<T, const M: usize, const N: usize> Add for Matrix<T, M, N>
251where
252 T: AddAssign + Copy + Num,
253{
254 type Output = Self;
255
256 /// Adds two matrices element-wise.
257 ///
258 /// # Examples
259 ///
260 /// ```
261 /// use mini_matrix::Matrix;
262 ///
263 /// let a = Matrix::<i32, 2, 2>::from([[1, 2], [3, 4]]);
264 /// let b = Matrix::<i32, 2, 2>::from([[5, 6], [7, 8]]);
265 /// let c = a + b;
266 /// assert_eq!(c.store, [[6, 8], [10, 12]]);
267 /// ```
268 fn add(self, rhs: Self) -> Self::Output {
269 let mut result = self.clone();
270 for (l_row, r_row) in result.store.iter_mut().zip(rhs.store.iter()) {
271 for (l, r) in l_row.iter_mut().zip(r_row.iter()) {
272 *l += *r;
273 }
274 }
275 result
276 }
277}
278
279impl<T, const M: usize, const N: usize> Sub for Matrix<T, M, N>
280where
281 T: SubAssign + Copy + Num,
282{
283 type Output = Self;
284
285 /// Subtracts one matrix from another element-wise.
286 ///
287 /// # Examples
288 ///
289 /// ```
290 /// use mini_matrix::Matrix;
291 ///
292 /// let a = Matrix::<i32, 2, 2>::from([[5, 6], [7, 8]]);
293 /// let b = Matrix::<i32, 2, 2>::from([[1, 2], [3, 4]]);
294 /// let c = a - b;
295 /// assert_eq!(c.store, [[4, 4], [4, 4]]);
296 /// ```
297 fn sub(self, rhs: Self) -> Self::Output {
298 let mut result = self.clone();
299 for (l_row, r_row) in result.store.iter_mut().zip(rhs.store.iter()) {
300 for (l, r) in l_row.iter_mut().zip(r_row.iter()) {
301 *l -= *r;
302 }
303 }
304 result
305 }
306}
307
308impl<T, const M: usize, const N: usize> Mul<T> for Matrix<T, M, N>
309where
310 T: MulAssign + Copy + Num,
311{
312 type Output = Self;
313
314 /// Multiplies a matrix by a scalar value.
315 ///
316 /// # Examples
317 ///
318 /// ```
319 /// use mini_matrix::Matrix;
320 ///
321 /// let a = Matrix::<i32, 2, 2>::from([[1, 2], [3, 4]]);
322 /// let b = a * 2;
323 /// assert_eq!(b.store, [[2, 4], [6, 8]]);
324 /// ```
325 fn mul(self, rhs: T) -> Self::Output {
326 let mut result = self.clone();
327 for row in result.store.iter_mut() {
328 for elem in row.iter_mut() {
329 *elem *= rhs;
330 }
331 }
332 result
333 }
334}
335
336impl<T, const M: usize, const N: usize> Mul<Vector<T, N>> for Matrix<T, M, N>
337where
338 T: MulAssign + AddAssign + Copy + Num + Default,
339{
340 type Output = Vector<T, M>;
341
342 /// Multiplies a matrix by a vector.
343 ///
344 /// # Examples
345 ///
346 /// ```
347 /// use mini_matrix::{Matrix, Vector};
348 ///
349 /// let a = Matrix::<i32, 2, 2>::from([[1, 2], [3, 4]]);
350 /// let v = Vector::<i32, 2>::from([5, 6]);
351 /// let result = a * v;
352 /// assert_eq!(result.store, [17, 39]);
353 /// ```
354 fn mul(self, rhs: Vector<T, N>) -> Self::Output {
355 let mut result = Vector::zero();
356 for (idx, row) in self.store.iter().enumerate() {
357 for (e1, e2) in row.iter().zip(rhs.store.iter()) {
358 result.store[idx] += *e1 * *e2;
359 }
360 }
361 result
362 }
363}
364
365impl<T, const M: usize, const N: usize> Mul<Matrix<T, N, N>> for Matrix<T, M, N>
366where
367 T: MulAssign + AddAssign + Copy + Num + Default,
368{
369 type Output = Self;
370
371 /// Multiplies two matrices.
372 ///
373 /// # Examples
374 ///
375 /// ```
376 /// use mini_matrix::Matrix;
377 ///
378 /// let a = Matrix::<i32, 2, 2>::from([[1, 2], [3, 4]]);
379 /// let b = Matrix::<i32, 2, 2>::from([[5, 6], [7, 8]]);
380 /// let c = a * b;
381 /// assert_eq!(c.store, [[17, 23], [39, 53]]);
382 /// ```
383 fn mul(self, rhs: Matrix<T, N, N>) -> Self::Output {
384 let mut result = Matrix::zero();
385 for (i, row) in self.store.iter().enumerate() {
386 for (j, col) in rhs.store.iter().enumerate() {
387 for (e1, e2) in row.iter().zip(col.iter()) {
388 result.store[i][j] += *e1 * *e2;
389 }
390 }
391 }
392 result
393 }
394}
395
396impl<T, const M: usize, const N: usize> Neg for Matrix<T, M, N>
397where
398 T: Neg<Output = T> + Copy,
399{
400 type Output = Self;
401
402 /// Negates all elements of the matrix.
403 ///
404 /// # Examples
405 ///
406 /// ```
407 /// use mini_matrix::Matrix;
408 ///
409 /// let a = Matrix::<i32, 2, 2>::from([[1, -2], [-3, 4]]);
410 /// let b = -a;
411 /// assert_eq!(b.store, [[-1, 2], [3, -4]]);
412 /// ````
413 fn neg(self) -> Self::Output {
414 let mut result = self.clone();
415 for row in result.store.iter_mut() {
416 for elem in row.iter_mut() {
417 *elem = -*elem;
418 }
419 }
420 result
421 }
422}
423
424impl<T, const M: usize, const N: usize> Display for Matrix<T, M, N>
425where
426 T: AddAssign + SubAssign + MulAssign + Copy + std::fmt::Display,
427{
428 /// Formats the matrix for display.
429 ///
430 /// Each row of the matrix is displayed on a new line, with elements separated by commas.
431 /// Elements are formatted with one decimal place precision.
432 ///
433 /// # Examples
434 ///
435 /// ```
436 /// use mini_matrix::Matrix;
437 ///
438 /// let a = Matrix::<f32, 2, 2>::from([[1.0, 2.5], [3.7, 4.2]]);
439 /// println!("{}", a);
440 /// // Output:
441 /// // // [1.0, 2.5]
442 /// // // [3.7, 4.2]
443 /// ```
444 fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
445 for (i, row) in self.store.iter().enumerate() {
446 if i > 0 {
447 writeln!(f)?;
448 }
449 write!(f, "// [")?;
450 for (j, val) in row.iter().enumerate() {
451 if j > 0 {
452 write!(f, ", ")?;
453 }
454 write!(f, "{:.1}", val)?;
455 }
456 write!(f, "]")?;
457 }
458 write!(f, "\n")?;
459 Ok(())
460 }
461}
462
463impl<T, const M: usize, const N: usize> Matrix<T, M, N>
464where
465 T: Num + Copy + AddAssign + Default,
466{
467 /// Multiplies the matrix by a vector.
468 ///
469 /// # Arguments
470 /// * `vec` - The vector to multiply with the matrix.
471 ///
472 /// # Returns
473 /// The resulting vector of the multiplication.
474 pub fn mul_vec(&mut self, vec: &Vector<T, N>) -> Vector<T, N> {
475 let mut result = Vector::zero();
476 for (idx, row) in self.store.iter_mut().enumerate() {
477 for (e1, e2) in row.iter_mut().zip(vec.store.iter()) {
478 result[idx] += *e1 * *e2;
479 }
480 }
481 result
482 }
483
484 /// Multiplies the matrix by another matrix.
485 ///
486 /// # Arguments
487 /// * `mat` - The matrix to multiply with.
488 ///
489 /// # Returns
490 /// The resulting matrix of the multiplication.
491 pub fn mul_mat(&mut self, mat: &Matrix<T, M, N>) -> Matrix<T, M, N> {
492 let mut result = Matrix::zero();
493 for i in 0..M {
494 for j in 0..N {
495 for k in 0..N {
496 result[(i, j)] += self[(i, k)] * mat[(k, j)];
497 }
498 }
499 }
500 result
501 }
502}
503
504impl<T, const M: usize, const N: usize> Matrix<T, M, N>
505where
506 T: Num + Copy + AddAssign + Default,
507{
508 /// Calculates the trace of the matrix.
509 ///
510 /// The trace is defined as the sum of the elements on the main diagonal.
511 ///
512 /// # Panics
513 ///
514 /// Panics if the matrix is not square (i.e., if M != N).
515 ///
516 /// # Examples
517 ///
518 /// ```
519 /// use mini_matrix::Matrix;
520 ///
521 /// let mut a = Matrix::<i32, 3, 3>::from([[1, 2, 3], [4, 5, 6], [7, 8, 9]]);
522 /// assert_eq!(a.trace(), 15);
523 /// ```
524 pub fn trace(&self) -> T {
525 assert!(M == N, "Matrix must be square to calculate trace");
526
527 let mut result = T::default();
528 for i in 0..M {
529 result += self[(i, i)];
530 }
531 result
532 }
533}
534
535/* ********************************************* */
536/* Exercise 09 - Transpose */
537/* ********************************************* */
538impl<T, const M: usize, const N: usize> Matrix<T, M, N>
539where
540 T: Copy + Default,
541{
542 /// Computes the transpose of the matrix.
543 ///
544 /// # Examples
545 ///
546 /// ```
547 /// use mini_matrix::Matrix;
548 ///
549 /// let mut a = Matrix::<i32, 2, 3>::from([[1, 2, 3], [4, 5, 6]]);
550 /// let b = a.transpose();
551 /// assert_eq!(b.store, [[1, 4], [2, 5], [3, 6]]);
552 /// ```
553 pub fn transpose(&mut self) -> Matrix<T, N, M> {
554 let mut result = Matrix::zero();
555 for i in 0..M {
556 for j in 0..N {
557 result[(j, i)] = self[(i, j)];
558 }
559 }
560 result
561 }
562}
563
564/* ********************************************** */
565/* Exercise XX - Identity */
566/* ********************************************** */
567impl<T, const M: usize, const N: usize> Matrix<T, M, N>
568where
569 T: Copy + Default + Num,
570{
571 /// Creates an identity matrix.
572 ///
573 /// # Examples
574 ///
575 /// ```
576 /// use mini_matrix::Matrix;
577 ///
578 /// let a = Matrix::<i32, 3, 3>::identity();
579 /// assert_eq!(a.store, [[1, 0, 0], [0, 1, 0], [0, 0, 1]]);
580 /// ```
581 pub fn identity() -> Matrix<T, M, N> {
582 let mut result = Matrix::zero();
583 for i in 0..M {
584 result[(i, i)] = T::one();
585 }
586 result
587 }
588}
589
590/* ************************************************* */
591/* Exercise 12 - Row Echelon */
592/* ************************************************* */
593
594impl<T, const M: usize, const N: usize> Matrix<T, M, N>
595where
596 T: Copy + Default + Mul<Output = T> + PartialEq + Num + Div<Output = T> + Sub<Output = T>,
597{
598 /// Converts a given matrix to its Reduced Row-Echelon Form (RREF).
599 ///
600 /// Row-echelon form is a specific form of a matrix that is particularly useful for solving
601 /// systems of linear equations and for understanding the properties of the matrix. To explain
602 /// it simply and in detail, let's go through what row-echelon form is, how to achieve it, and
603 /// an example to illustrate the process.
604 ///
605 /// # Row-Echelon Form
606 ///
607 /// A matrix is in row-echelon form if it satisfies the following conditions:
608 ///
609 /// 1. **Leading Entry**: The leading entry (first non-zero number from the left) in each row is 1.
610 /// This is called the pivot.
611 /// 2. **Zeros Below**: The pivot in each row is to the right of the pivot in the row above, and
612 /// all entries below a pivot are zeros.
613 /// 3. **Rows of Zeros**: Any rows consisting entirely of zeros are at the bottom of the matrix.
614 ///
615 /// # Reduced Row-Echelon Form
616 ///
617 /// A matrix is in reduced row-echelon form (RREF) if it additionally satisfies:
618 ///
619 /// 4. **Leading Entries**: Each leading 1 is the only non-zero entry in its column.
620 ///
621 /// # Achieving Row-Echelon Form
622 ///
623 /// To convert a matrix into row-echelon form, we use a process called Gaussian elimination.
624 /// This involves performing row operations:
625 ///
626 /// 1. **Row swapping**: Swapping the positions of two rows.
627 /// 2. **Row multiplication**: Multiplying all entries of a row by a non-zero scalar.
628 /// 3. **Row addition**: Adding or subtracting the multiples of one row to another row.
629 ///
630 ///
631 /// Let's consider an example with a `3 x 3` matrix:
632 ///
633 ///
634 /// A = [
635 /// [1, 2, 3],
636 /// [4, 5, 6],
637 /// [7, 8, 9]
638 /// ]
639 ///
640 /// ## Step-by-Step Process
641 ///
642 /// 1. **Starting Matrix**:
643 ///
644 /// [
645 /// [1, 2, 3],
646 /// [4, 5, 6],
647 /// [7, 8, 9]
648 /// ]
649 ///
650 /// 2. **Make the Pivot of Row 1 (already 1)**:
651 ///
652 /// The first leading entry is already 1.
653 ///
654 /// 3. **Eliminate Below Pivot in Column 1**:
655 ///
656 /// Subtract 4 times the first row from the second row:
657 /// ```markdown
658 /// R2 = R2 - 4R1
659 /// [
660 /// [1, 2, 3],
661 /// [0, -3, -6],
662 /// [7, 8, 9]
663 /// ]
664 /// ```
665 /// Subtract 7 times the first row from the third row:
666 /// ```markdown
667 /// R3 = R3 - 7R1
668 /// [
669 /// [1, 2, 3],
670 /// [0, -3, -6],
671 /// [0, -6, -12]
672 /// ]
673 /// ```
674 /// 4. **Make the Pivot of Row 2**:
675 ///
676 /// Divide the second row by -3 to make the pivot 1:
677 /// ```markdown
678 /// R2 = (1 / -3) * R2
679 /// [
680 /// [1, 2, 3],
681 /// [0, 1, 2],
682 /// [0, -6, -12]
683 /// ]
684 ///
685 /// 5. **Eliminate Below Pivot in Column 2**:
686 ///
687 /// Add 6 times the second row to the third row:
688 /// ```markdown
689 /// R3 = R3 + 6R2
690 /// [
691 /// [1, 2, 3],
692 /// [0, 1, 2],
693 /// [0, 0, 0]
694 /// ]
695 /// ```
696 /// Now, the matrix is in row-echelon form.
697 ///
698 /// # Examples
699 ///
700 /// ```
701 /// use mini_matrix::Matrix;
702 ///
703 /// let a = Matrix::<f64, 3, 4>::from([
704 /// [1.0, 2.0, 3.0, 4.0],
705 /// [5.0, 6.0, 7.0, 8.0],
706 /// [9.0, 10.0, 11.0, 12.0]
707 /// ]);
708 /// let b = a.row_echelon();
709 /// // Check the result (approximate due to floating-point arithmetic)
710 /// ```
711 pub fn row_echelon(&self) -> Matrix<T, M, N> {
712 let mut result = self.clone();
713 // let mut matrix_out = result.store;
714 let mut pivot = 0;
715 let row_count = M;
716 let column_count = N;
717
718 'outer: for r in 0..row_count {
719 if column_count <= pivot {
720 break;
721 }
722 let mut i = r;
723 while result[(i, pivot)] == T::default() {
724 i = i + 1;
725 if i == row_count {
726 i = r;
727 pivot = pivot + 1;
728 if column_count == pivot {
729 pivot = pivot - 1;
730 break 'outer;
731 }
732 }
733 }
734 for j in 0..row_count {
735 let temp = result[(r, j)];
736 result[(r, j)] = result[(i, j)];
737 result[(i, j)] = temp;
738 }
739 let divisor = result[(r, pivot)];
740 if divisor != T::default() {
741 for j in 0..column_count {
742 result[(r, j)] = result[(r, j)] / divisor;
743 }
744 }
745 for j in 0..row_count {
746 if j != r {
747 let hold = result[(j, pivot)];
748 for k in 0..column_count {
749 result[(j, k)] = result[(j, k)] - (hold * result[(r, k)]);
750 }
751 }
752 }
753 pivot = pivot + 1;
754 }
755 result
756 }
757}
758
759/************************************************ * */
760/* Exercise 12 - Determinant */
761/************************************************ */
762
763impl<T, const M: usize, const N: usize> Matrix<T, M, N>
764where
765 T: Copy + Default + Mul + Num + Neg<Output = T> + AddAssign + Debug,
766{
767 /// Computes the determinant of the matrix.
768 ///
769 /// # Determinant in General
770 ///
771 /// The determinant is a scalar value that can be computed from the elements of a square matrix.
772 /// It provides important properties about the matrix and the linear transformation it represents.
773 /// In general, the determinant represents the scaling factor of the volume when the matrix is
774 /// used as a linear transformation. It can be positive, negative, or zero, each with different
775 /// implications:
776 ///
777 /// - **\(\det(A) = 0\)**:
778 /// - The matrix `A` is **singular** and does not have an inverse.
779 /// - The columns (or rows) of `A` are linearly dependent.
780 /// - The transformation collapses the space into a lower-dimensional subspace.
781 /// - Geometrically, it indicates that the volume of the transformed space is 0.
782 ///
783 /// - **\(\det(A) > 0\)**:
784 /// - The matrix `A` is **non-singular** and has an inverse.
785 /// - The transformation preserves the orientation of the space.
786 /// - Geometrically, it indicates a positive scaling factor of the volume.
787 ///
788 /// - **\(\det(A) < 0\)**:
789 /// - The matrix `A` is **non-singular** and has an inverse.
790 /// - The transformation reverses the orientation of the space.
791 /// - Geometrically, it indicates a negative scaling factor of the volume.
792 ///
793 /// # Example
794 ///
795 /// Consider a `2 x 2` matrix:
796 ///
797 /// ```text
798 /// A = [
799 /// [1, 2],
800 /// [3, 4]
801 /// ]
802 /// ```
803 ///
804 /// The determinant is:
805 ///
806 /// ```text
807 /// det(A) = 1 * 4 - 2 * 3 = 4 - 6 = -2
808 /// ```
809 ///
810 /// This indicates that the transformation represented by `A` scales areas by a factor of 2 and
811 /// reverses their orientation.
812 ///
813 /// # Panics
814 ///
815 /// Panics if the matrix size is larger than `4 x 4`.
816 ///
817 /// # Returns
818 ///
819 /// The determinant of the matrix.
820 ///
821 /// # Examples
822 ///
823 /// ```
824 /// use mini_matrix::Matrix;
825 ///
826 /// let a = Matrix::<i32, 3, 3>::from([[1, 2, 3], [4, 5, 6], [7, 8, 9]]);
827 /// assert_eq!(a.determinant(), 0);
828 /// ```
829 pub fn determinant(&self) -> T {
830 match M {
831 1 => self[(0, 0)],
832 2 => self[(0, 0)] * self[(1, 1)] - self[(0, 1)] * self[(1, 0)],
833 3 => self.determinant_3x3(),
834 4 => (0..4)
835 .map(|i| {
836 let sign = if i % 2 == 0 { T::one() } else { -T::one() };
837 let cofactor = self.get_cofactor(0, i);
838 sign * self[(0, i)] * cofactor.determinant_3x3()
839 })
840 .fold(T::default(), |acc, x| acc + x),
841 _ => panic!("Determinant not implemented for matrices larger than 4x4"),
842 }
843 }
844
845 fn determinant_3x3(&self) -> T {
846 self[(0, 0)] * (self[(1, 1)] * self[(2, 2)] - self[(1, 2)] * self[(2, 1)])
847 - self[(0, 1)] * (self[(1, 0)] * self[(2, 2)] - self[(1, 2)] * self[(2, 0)])
848 + self[(0, 2)] * (self[(1, 0)] * self[(2, 1)] - self[(1, 1)] * self[(2, 0)])
849 }
850
851 fn get_cofactor(&self, row: usize, col: usize) -> Matrix<T, 3, 3> {
852 let mut cofactor_matrix = Matrix::<T, 3, 3>::zero();
853 let mut row_index = 0;
854
855 for r in 0..M {
856 if r == row {
857 continue;
858 }
859 let mut col_index = 0;
860 for c in 0..N {
861 if c == col {
862 continue;
863 }
864 cofactor_matrix[(row_index, col_index)] = self[(r, c)];
865 col_index += 1;
866 }
867 row_index += 1;
868 }
869 cofactor_matrix
870 }
871}
872
873/********************************************* */
874/* Exercise 12 - Inverse */
875/********************************************* */
876
877impl<T, const M: usize, const N: usize> Matrix<T, M, N>
878where
879 T: Copy + Default + Mul + Num + Neg<Output = T> + AddAssign + Debug + Float,
880{
881 /// Calculates the inverse of the matrix.
882 ///
883 /// This method supports matrices up to 3x3 in size.
884 ///
885 /// # Returns
886 ///
887 /// Returns `Ok(Matrix)` if the inverse exists, or an `Err` with a descriptive message if not.
888 ///
889 /// # Examples
890 ///
891 /// ```
892 /// use mini_matrix::Matrix;
893 ///
894 /// let a = Matrix::<f64, 2, 2>::from([[1.0, 2.0], [3.0, 4.0]]);
895 /// let inv = a.inverse().unwrap();
896 /// // Check the result (approximate due to floating-point arithmetic)
897 /// ```
898 pub fn inverse(&self) -> Result<Self, &'static str> {
899 if M != N {
900 return Err("Matrix must be square to calculate inverse");
901 }
902
903 let det = self.determinant();
904
905 if det == T::zero() {
906 return Err("Matrix is singular and has no inverse");
907 }
908
909 let mut inv = Matrix::<T, N, M>::zero();
910 for i in 0..M {
911 for j in 0..N {
912 let coffactor = match M {
913 2 => self.cofactor1x1(i, j).determinant(),
914 3 => self.cofactor2x2(i, j).determinant(),
915 _ => return Err("Inverse not implemented for matrices larger than 3x3"),
916 };
917 let base: i32 = -1;
918 inv[(i, j)] = (coffactor * T::from(base.pow((i + j) as u32)).unwrap()) / det;
919 }
920 }
921 let inv = inv.transpose();
922 Ok(inv)
923 }
924}
925
926/***************************************** */
927/* Exercise 13 - Rank */
928/***************************************** */
929
930impl<T, const M: usize, const N: usize> Matrix<T, M, N>
931where
932 T: Copy + Default + Mul + Num + Neg<Output = T> + AddAssign + PartialEq,
933{
934 /// Calculates the rank of the matrix.
935 ///
936 /// The rank is determined by computing the row echelon form and counting non-zero rows.
937 ///
938 /// # Examples
939 ///
940 /// ```
941 /// use mini_matrix::Matrix;
942 ///
943 /// let a = Matrix::<i32, 3, 3>::from([[1, 2, 3], [4, 5, 6], [7, 8, 9]]);
944 /// assert_eq!(a.rank(), 2);
945 /// ```
946 pub fn rank(&self) -> usize {
947 let mut rank = M;
948 let row_echelon = self.row_echelon();
949 for i in 0..M {
950 let mut is_zero = true;
951 for j in 0..N {
952 if row_echelon[(i, j)] != T::default() {
953 is_zero = false;
954 break;
955 }
956 }
957 if is_zero {
958 rank -= 1;
959 }
960 }
961 rank
962 }
963}