new: Added det_for_field method

This commit is contained in:
Sayantan Santra 2023-05-25 22:59:01 -05:00
parent 19e410322e
commit 3f8373164b
Signed by: SinTan1729
GPG key ID: EB3E68BFBA25C85F
2 changed files with 67 additions and 4 deletions

View file

@ -13,7 +13,7 @@ use num::{
}; };
use std::{ use std::{
fmt::{self, Debug, Display, Formatter}, fmt::{self, Debug, Display, Formatter},
ops::{Add, Mul, Sub}, ops::{Add, Div, Mul, Sub},
result::Result, result::Result,
}; };
@ -129,7 +129,8 @@ impl<T: Mul + Add + Sub> Matrix<T> {
/// Return the determinant of a square matrix. This method additionally requires [`Zero`], /// Return the determinant of a square matrix. This method additionally requires [`Zero`],
/// [`One`] and [`Copy`] traits. Also, we need that the [`Mul`] and [`Add`] operations /// [`One`] and [`Copy`] traits. Also, we need that the [`Mul`] and [`Add`] operations
/// return the same type `T`. /// return the same type `T`. This uses basic recursive algorithm using cofactor-minor.
/// See [`det_in_field`](Self::det_in_field()) for faster determinant calculation in fields.
/// It'll throw an error if the provided matrix isn't square. /// It'll throw an error if the provided matrix isn't square.
/// # Example /// # Example
/// ``` /// ```
@ -146,7 +147,7 @@ impl<T: Mul + Add + Sub> Matrix<T> {
{ {
if self.is_square() { if self.is_square() {
// It's a recursive algorithm using minors. // It's a recursive algorithm using minors.
// TODO: Implement a faster algorithm. Maybe use row reduction for fields. // TODO: Implement a faster algorithm.
let out = if self.width() == 1 { let out = if self.width() == 1 {
self.entries[0][0] self.entries[0][0]
} else { } else {
@ -168,6 +169,61 @@ impl<T: Mul + Add + Sub> Matrix<T> {
} }
} }
/// Return the determinant of a square matrix over a field i.e. needs [`One`] and [`Div`] traits.
/// See [`det`](Self::det()) for determinants in rings.
/// This method uses row reduction as is much faster.
/// It'll throw an error if the provided matrix isn't square.
/// # Example
/// ```
/// use matrix_basic::Matrix;
/// let m = Matrix::from(vec![vec![1,2],vec![3,4]]).unwrap();
/// assert_eq!(m.det(),Ok(-2));
/// ```
pub fn det_in_field(&self) -> Result<T, &'static str>
where
T: Copy,
T: Mul<Output = T>,
T: Sub<Output = T>,
T: Zero,
T: One,
T: PartialEq,
T: Div<Output = T>,
{
if self.is_square() {
// Cloning is necessary as we'll be doing row operations on it.
let mut rows = self.entries.clone();
let mut multiplier = T::one();
for i in 0..self.height() {
// First check if the row has diagonal element 0, if yes, then swap.
if rows[i][i] == T::zero() {
let mut zero_column = true;
for j in (i + 1)..self.height() {
if rows[j][i] != T::zero() {
rows.swap(i, j);
multiplier = T::zero() - multiplier;
zero_column = false;
break;
}
}
if zero_column {
return Ok(T::zero());
}
}
for j in (i + 1)..self.height() {
for k in (i + 1)..self.width() {
rows[j][k] = rows[j][k] - rows[i][k] * rows[j][i] / rows[i][i];
}
}
}
for (i, row) in rows.iter().enumerate() {
multiplier = multiplier * row[i];
}
Ok(multiplier)
} else {
Err("Provided matrix isn't square.")
}
}
/// Creates a zero matrix of a given size. /// Creates a zero matrix of a given size.
pub fn zero(height: usize, width: usize) -> Self pub fn zero(height: usize, width: usize) -> Self
where where
@ -213,7 +269,7 @@ impl<T: Debug + Mul + Add + Sub> Display for Matrix<T> {
} }
impl<T: Mul<Output = T> + Add + Sub + Copy + Zero> Mul for Matrix<T> { impl<T: Mul<Output = T> + Add + Sub + Copy + Zero> Mul for Matrix<T> {
// TODO: Implement a faster algorithm. Maybe use row reduction for fields. // TODO: Implement a faster algorithm.
type Output = Self; type Output = Self;
fn mul(self, other: Self) -> Self { fn mul(self, other: Self) -> Self {
let width = self.width(); let width = self.width();

View file

@ -22,7 +22,14 @@ fn add_sub_test() {
fn det_test() { fn det_test() {
let a = Matrix::from(vec![vec![1, 2, 0], vec![0, 3, 5], vec![0, 0, 10]]).unwrap(); let a = Matrix::from(vec![vec![1, 2, 0], vec![0, 3, 5], vec![0, 0, 10]]).unwrap();
let b = Matrix::from(vec![vec![1, 2, 0], vec![0, 3, 5]]).unwrap(); let b = Matrix::from(vec![vec![1, 2, 0], vec![0, 3, 5]]).unwrap();
let c = Matrix::from(vec![
vec![0.0, 0.0, 10.0],
vec![0.0, 3.0, 5.0],
vec![1.0, 2.0, 0.0],
])
.unwrap();
assert_eq!(a.det(), Ok(30)); assert_eq!(a.det(), Ok(30));
assert_eq!(c.det_in_field(), Ok(-30.0));
assert!(b.det().is_err()); assert!(b.det().is_err());
} }