use super::*;
#[repr(transparent)]
pub struct Matrix<T, const N: usize, const M: usize>(pub(crate) [Vector<T, { N }>; M]);
impl<T, const N: usize, const M: usize> Matrix<T, { N }, { M }> {
pub fn swap_columns(&mut self, a: usize, b: usize) {
unsafe { core::ptr::swap(&mut self.0[a], &mut self.0[b]) };
}
pub fn swap_rows(&mut self, a: usize, b: usize) {
for v in self.0.iter_mut() {
unsafe { core::ptr::swap(&mut v[a], &mut v[b]) };
}
}
pub fn swap_elements(&mut self, (acol, arow): (usize, usize), (bcol, brow): (usize, usize)) {
unsafe { core::ptr::swap(&mut self[acol][arow], &mut self[bcol][brow]) };
}
pub fn column_iter<'a>(&'a self) -> impl Iterator<Item = &'a Vector<T, { N }>> {
self.0.iter()
}
pub fn column_iter_mut<'a>(&'a mut self) -> impl Iterator<Item = &'a mut Vector<T, { N }>> {
self.0.iter_mut()
}
pub fn row_iter<'a>(&'a self) -> impl Iterator<Item = RowView<'a, T, { N }, { M }>> {
RowIter {
row: 0,
matrix: self,
}
}
pub fn row_iter_mut<'a>(&'a mut self) -> impl Iterator<Item = RowViewMut<'a, T, { N }, { M }>> {
RowIterMut {
row: 0,
matrix: self,
phantom: PhantomData,
}
}
pub fn map<Out, F>(self, mut f: F) -> Matrix<Out, { N }, { M }>
where
F: FnMut(T) -> Out,
{
let mut from = MaybeUninit::new(self);
let mut to = MaybeUninit::<Matrix<Out, { N }, { M }>>::uninit();
let fromp: *mut MaybeUninit<Vector<T, { N }>> = unsafe { mem::transmute(&mut from) };
let top: *mut Vector<Out, { N }> = unsafe { mem::transmute(&mut to) };
for i in 0..M {
unsafe {
let fromp: *mut MaybeUninit<T> = mem::transmute(fromp.add(i));
let top: *mut Out = mem::transmute(top.add(i));
for j in 0..N {
top.add(j)
.write(f(fromp.add(j).replace(MaybeUninit::uninit()).assume_init()));
}
}
}
unsafe { to.assume_init() }
}
pub fn transpose(self) -> Matrix<T, { M }, { N }> {
let mut from = MaybeUninit::new(self);
let mut trans = MaybeUninit::<[Vector<T, { M }>; N]>::uninit();
let fromp: *mut Vector<MaybeUninit<T>, { N }> = unsafe { mem::transmute(&mut from) };
let transp: *mut Vector<T, { M }> = unsafe { mem::transmute(&mut trans) };
for j in 0..N {
let mut row = MaybeUninit::<[T; M]>::uninit();
let rowp: *mut T = unsafe { mem::transmute(&mut row) };
for k in 0..M {
unsafe {
let fromp: *mut MaybeUninit<T> = mem::transmute(fromp.add(k));
rowp.add(k)
.write(fromp.add(j).replace(MaybeUninit::uninit()).assume_init());
}
}
let row = Vector::<T, { M }>::from(unsafe { row.assume_init() });
unsafe {
transp.add(j).write(row);
}
}
Matrix::<T, { M }, { N }>(unsafe { trans.assume_init() })
}
}
impl<T, const N: usize> Matrix<T, { N }, { N }>
where
T: Clone,
{
pub fn diagonal(&self) -> Vector<T, { N }> {
let mut diag = MaybeUninit::<[T; N]>::uninit();
let diagp: *mut T = unsafe { mem::transmute(&mut diag) };
for i in 0..N {
unsafe {
diagp.add(i).write(self.0[i].0[i].clone());
}
}
Vector::<T, { N }>(unsafe { diag.assume_init() })
}
}
impl<T, const N: usize> Matrix<T, { N }, { N }>
where
T: Clone + PartialOrd + Product + Real + One + Zero,
T: Neg<Output = T>,
T: Add<T, Output = T> + Sub<T, Output = T>,
T: Mul<T, Output = T> + Div<T, Output = T>,
Self: Add<Self>,
Self: Sub<Self>,
Self: Mul<Self>,
Self: Mul<Vector<T, { N }>, Output = Vector<T, { N }>>,
{
pub fn lu(mut self) -> Option<LU<T, { N }>> {
let mut p = Permutation::<{ N }>::unit();
for i in 0..N {
let mut max_a = T::zero();
let mut imax = i;
for k in i..N {
let abs = self[i][k].clone().abs();
if abs > max_a {
max_a = abs;
imax = k;
}
}
if max_a.is_zero() {
return None;
}
if imax != i {
p.swap(i, imax);
self.swap_rows(i, imax);
}
for j in i + 1..N {
self[(j, i)] = self[(j, i)].clone() / self[(i, i)].clone();
for k in i + 1..N {
self[(j, k)] =
self[(j, k)].clone() - self[(j, i)].clone() * self[(i, k)].clone();
}
}
}
Some(LU(p, self))
}
pub fn determinant(&self) -> T {
self.clone().lu().map_or(T::zero(), |x| x.determinant())
}
pub fn invert(self) -> Option<Self> {
self.lu().map(|x| x.invert())
}
}
impl<T, const N: usize, const M: usize> From<[Vector<T, { N }>; M]>
for Matrix<T, { N }, { M }>
{
fn from(array: [Vector<T, { N }>; M]) -> Self {
Matrix::<T, { N }, { M }>(array)
}
}
impl<T, const N: usize, const M: usize> From<[[T; N]; M]> for Matrix<T, { N }, { M }> {
fn from(array: [[T; N]; M]) -> Self {
let mut array = MaybeUninit::<[[T; N]; M]>::new(array);
let mut vec_array: MaybeUninit<[Vector<T, { N }>; M]> = MaybeUninit::uninit();
let arrayp: *mut MaybeUninit<[T; N]> = unsafe { mem::transmute(&mut array) };
let vec_arrayp: *mut Vector<T, { N }> = unsafe { mem::transmute(&mut vec_array) };
for i in 0..M {
unsafe {
vec_arrayp.add(i).write(Vector::<T, { N }>(
arrayp.add(i).replace(MaybeUninit::uninit()).assume_init(),
));
}
}
Matrix::<T, { N }, { M }>(unsafe { vec_array.assume_init() })
}
}
impl<T> From<Quaternion<T>> for Matrix3<T>
where
T: Add + Mul + Sub + Real + One + Copy + Clone,
{
fn from(quat: Quaternion<T>) -> Self {
let x2 = quat.v.x() + quat.v.x();
let y2 = quat.v.y() + quat.v.y();
let z2 = quat.v.z() + quat.v.z();
let xx2 = x2 * quat.v.x();
let xy2 = x2 * quat.v.y();
let xz2 = x2 * quat.v.z();
let yy2 = y2 * quat.v.y();
let yz2 = y2 * quat.v.z();
let zz2 = z2 * quat.v.z();
let sy2 = y2 * quat.s;
let sz2 = z2 * quat.s;
let sx2 = x2 * quat.s;
matrix![
[T::one() - yy2 - zz2, xy2 + sz2, xz2 - sy2],
[xy2 - sz2, T::one() - xx2 - zz2, yz2 + sx2],
[xz2 + sy2, yz2 - sx2, T::one() - xx2 - yy2],
]
}
}
pub type Matrix2<T> = Matrix<T, 2, 2>;
pub type Matrix3<T> = Matrix<T, 3, 3>;
pub type Matrix4<T> = Matrix<T, 4, 4>;
#[inline]
#[doc(hidden)]
pub fn new_matrix<T: Clone, const N: usize, const M: usize>(
rows: [[T; M]; N],
) -> Matrix<T, { N }, { M }> {
Matrix::<T, { M }, { N }>::from(rows).transpose()
}
#[macro_export]
macro_rules! matrix {
( $item:expr ) => {
$crate::new_matrix([
[ $item ]
])
};
( $($rows:expr),* $(,)? ) => {
$crate::new_matrix([
$($rows),*
])
};
}
impl<T, const N: usize, const M: usize> Clone for Matrix<T, { N }, { M }>
where
T: Clone,
{
fn clone(&self) -> Self {
Matrix::<T, { N }, { M }>(self.0.clone())
}
}
impl<T, const N: usize, const M: usize> Copy for Matrix<T, { N }, { M }> where T: Copy {}
impl<T, const N: usize, const M: usize> Deref for Matrix<T, { N }, { M }> {
type Target = [Vector<T, { N }>; M];
fn deref(&self) -> &Self::Target {
&self.0
}
}
impl<T, const N: usize, const M: usize> DerefMut for Matrix<T, { N }, { M }> {
fn deref_mut(&mut self) -> &mut Self::Target {
&mut self.0
}
}
impl<T, const N: usize, const M: usize> Hash for Matrix<T, { N }, { M }>
where
T: Hash,
{
fn hash<H: Hasher>(&self, state: &mut H) {
for i in 0..M {
self.0[i].hash(state);
}
}
}
impl<T, const N: usize, const M: usize> FromIterator<T> for Matrix<T, { N }, { M }> {
fn from_iter<I>(iter: I) -> Self
where
I: IntoIterator<Item = T>,
{
let mut iter = iter.into_iter();
let mut new = MaybeUninit::<[Vector<T, { N }>; M]>::uninit();
let newp: *mut Vector<T, { N }> = unsafe { mem::transmute(&mut new) };
for i in 0..M {
let mut newv = MaybeUninit::<Vector<T, { N }>>::uninit();
let newvp: *mut T = unsafe { mem::transmute(&mut newv) };
for j in 0..N {
if let Some(next) = iter.next() {
unsafe { newvp.add(j).write(next) };
} else {
panic!(
"too few items in iterator to create Matrix<_, {}, {}>",
N, M
);
}
}
unsafe {
newp.add(i)
.write(mem::replace(&mut newv, MaybeUninit::uninit()).assume_init());
}
}
if iter.next().is_some() {
panic!(
"too many items in iterator to create Matrix<_, {}, {}>",
N, M
);
}
Matrix::<T, { N }, { M }>(unsafe { new.assume_init() })
}
}
impl<T, const N: usize, const M: usize> FromIterator<Vector<T, { N }>> for Matrix<T, { N }, { M }> {
fn from_iter<I>(iter: I) -> Self
where
I: IntoIterator<Item = Vector<T, { N }>>,
{
let mut iter = iter.into_iter();
let mut new = MaybeUninit::<[Vector<T, { N }>; M]>::uninit();
let newp: *mut Vector<T, { N }> = unsafe { mem::transmute(&mut new) };
for i in 0..M {
if let Some(v) = iter.next() {
unsafe {
newp.add(i).write(v);
}
} else {
panic!(
"too few items in iterator to create Matrix<_, {}, {}>",
N, M
);
}
}
Matrix::<T, { N }, { M }>(unsafe { new.assume_init() })
}
}
impl<T, const N: usize, const M: usize> IntoIterator for Matrix<T, { N }, { M }> {
type Item = Vector<T, { N }>;
type IntoIter = ArrayIter<Vector<T, { N }>, { M }>;
fn into_iter(self) -> Self::IntoIter {
let Matrix(array) = self;
ArrayIter {
array: MaybeUninit::new(array),
pos: 0,
}
}
}
impl<T, const N: usize, const M: usize> Zero for Matrix<T, { N }, { M }>
where
T: Zero,
Vector<T, { N }>: Zero,
{
fn zero() -> Self {
let mut zero_mat = MaybeUninit::<[Vector<T, { N }>; M]>::uninit();
let matp: *mut Vector<T, { N }> = unsafe { mem::transmute(&mut zero_mat) };
for i in 0..M {
unsafe {
matp.add(i).write(Vector::<T, { N }>::zero());
}
}
Matrix::<T, { N }, { M }>(unsafe { zero_mat.assume_init() })
}
fn is_zero(&self) -> bool {
for i in 0..M {
if !self.0[i].is_zero() {
return false;
}
}
true
}
}
impl<T, const N: usize> One for Matrix<T, { N }, { N }>
where
T: Zero + One + Clone,
Self: PartialEq<Self>,
{
fn one() -> Self {
let mut unit_mat = MaybeUninit::<[Vector<T, { N }>; N]>::uninit();
let matp: *mut Vector<T, { N }> = unsafe { mem::transmute(&mut unit_mat) };
for i in 0..N {
let mut unit_vec = MaybeUninit::<Vector<T, { N }>>::uninit();
let vecp: *mut T = unsafe { mem::transmute(&mut unit_vec) };
for j in 0..i {
unsafe {
vecp.add(j).write(<T as Zero>::zero());
}
}
unsafe {
vecp.add(i).write(<T as One>::one());
}
for j in (i + 1)..N {
unsafe {
vecp.add(j).write(<T as Zero>::zero());
}
}
unsafe {
matp.add(i).write(unit_vec.assume_init());
}
}
Matrix::<T, { N }, { N }>(unsafe { unit_mat.assume_init() })
}
fn is_one(&self) -> bool {
self == &<Self as One>::one()
}
}
impl<T, const N: usize, const M: usize> Index<usize> for Matrix<T, { N }, { M }> {
type Output = Vector<T, { N }>;
fn index(&self, column: usize) -> &Self::Output {
&self.0[column]
}
}
impl<T, const N: usize, const M: usize> IndexMut<usize> for Matrix<T, { N }, { M }> {
fn index_mut(&mut self, column: usize) -> &mut Self::Output {
&mut self.0[column]
}
}
impl<T, const N: usize, const M: usize> Index<(usize, usize)> for Matrix<T, { N }, { M }> {
type Output = T;
fn index(&self, (row, column): (usize, usize)) -> &Self::Output {
&self.0[column][row]
}
}
impl<T, const N: usize, const M: usize> IndexMut<(usize, usize)> for Matrix<T, { N }, { M }> {
fn index_mut(&mut self, (row, column): (usize, usize)) -> &mut Self::Output {
&mut self.0[column][row]
}
}
impl<A, B, RHS, const N: usize, const M: usize> PartialEq<RHS> for Matrix<A, { N }, { M }>
where
RHS: Deref<Target = [Vector<B, { N }>; M]>,
A: PartialEq<B>,
{
fn eq(&self, other: &RHS) -> bool {
for (a, b) in self.0.iter().zip(other.deref().iter()) {
if !a.eq(b) {
return false;
}
}
true
}
}
impl<T, const N: usize, const M: usize> fmt::Debug for Matrix<T, { N }, { M }>
where
T: fmt::Debug,
{
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
write!(f, "Matrix [ ")?;
for i in 0..N {
write!(f, "[ ")?;
for j in 0..M {
write!(f, "{:?} ", self.0[j].0[i])?;
}
write!(f, "] ")?;
}
write!(f, "]")
}
}
impl<A, B, const N: usize, const M: usize> Add<Matrix<B, { N }, { M }>> for Matrix<A, { N }, { M }>
where
A: Add<B>,
{
type Output = Matrix<<A as Add<B>>::Output, { N }, { M }>;
fn add(self, rhs: Matrix<B, { N }, { M }>) -> Self::Output {
let mut mat = MaybeUninit::<[Vector<<A as Add<B>>::Output, { N }>; M]>::uninit();
let mut lhs = MaybeUninit::new(self);
let mut rhs = MaybeUninit::new(rhs);
let matp: *mut Vector<<A as Add<B>>::Output, { N }> = unsafe { mem::transmute(&mut mat) };
let lhsp: *mut MaybeUninit<Vector<A, { N }>> = unsafe { mem::transmute(&mut lhs) };
let rhsp: *mut MaybeUninit<Vector<B, { N }>> = unsafe { mem::transmute(&mut rhs) };
for i in 0..M {
unsafe {
matp.add(i).write(
lhsp.add(i).replace(MaybeUninit::uninit()).assume_init()
+ rhsp.add(i).replace(MaybeUninit::uninit()).assume_init(),
);
}
}
Matrix::<<A as Add<B>>::Output, { N }, { M }>(unsafe { mat.assume_init() })
}
}
impl<A, B, const N: usize, const M: usize> AddAssign<Matrix<B, { N }, { M }>>
for Matrix<A, { N }, { M }>
where
A: AddAssign<B>,
{
fn add_assign(&mut self, rhs: Matrix<B, { N }, { M }>) {
let mut rhs = MaybeUninit::new(rhs);
let rhsp: *mut MaybeUninit<Vector<B, { N }>> = unsafe { mem::transmute(&mut rhs) };
for i in 0..M {
self.0[i] += unsafe { rhsp.add(i).replace(MaybeUninit::uninit()).assume_init() };
}
}
}
impl<A, B, const N: usize, const M: usize> Sub<Matrix<B, { N }, { M }>> for Matrix<A, { N }, { M }>
where
A: Sub<B>,
{
type Output = Matrix<<A as Sub<B>>::Output, { N }, { M }>;
fn sub(self, rhs: Matrix<B, { N }, { M }>) -> Self::Output {
let mut mat = MaybeUninit::<[Vector<<A as Sub<B>>::Output, { N }>; M]>::uninit();
let mut lhs = MaybeUninit::new(self);
let mut rhs = MaybeUninit::new(rhs);
let matp: *mut Vector<<A as Sub<B>>::Output, { N }> = unsafe { mem::transmute(&mut mat) };
let lhsp: *mut MaybeUninit<Vector<A, { N }>> = unsafe { mem::transmute(&mut lhs) };
let rhsp: *mut MaybeUninit<Vector<B, { N }>> = unsafe { mem::transmute(&mut rhs) };
for i in 0..M {
unsafe {
matp.add(i).write(
lhsp.add(i).replace(MaybeUninit::uninit()).assume_init()
- rhsp.add(i).replace(MaybeUninit::uninit()).assume_init(),
);
}
}
Matrix::<<A as Sub<B>>::Output, { N }, { M }>(unsafe { mat.assume_init() })
}
}
impl<A, B, const N: usize, const M: usize> SubAssign<Matrix<B, { N }, { M }>>
for Matrix<A, { N }, { M }>
where
A: SubAssign<B>,
{
fn sub_assign(&mut self, rhs: Matrix<B, { N }, { M }>) {
let mut rhs = MaybeUninit::new(rhs);
let rhsp: *mut MaybeUninit<Vector<B, { N }>> = unsafe { mem::transmute(&mut rhs) };
for i in 0..M {
self.0[i] -= unsafe { rhsp.add(i).replace(MaybeUninit::uninit()).assume_init() };
}
}
}
impl<T, const N: usize, const M: usize> Neg for Matrix<T, { N }, { M }>
where
T: Neg,
{
type Output = Matrix<<T as Neg>::Output, { N }, { M }>;
fn neg(self) -> Self::Output {
let mut from = MaybeUninit::new(self);
let mut mat = MaybeUninit::<[Vector<<T as Neg>::Output, { N }>; M]>::uninit();
let fromp: *mut MaybeUninit<Vector<T, { N }>> = unsafe { mem::transmute(&mut from) };
let matp: *mut Vector<<T as Neg>::Output, { N }> = unsafe { mem::transmute(&mut mat) };
for i in 0..M {
unsafe {
matp.add(i).write(
fromp
.add(i)
.replace(MaybeUninit::uninit())
.assume_init()
.neg(),
);
}
}
Matrix::<<T as Neg>::Output, { N }, { M }>(unsafe { mat.assume_init() })
}
}
impl<T, const N: usize, const M: usize, const P: usize> Mul<Matrix<T, { M }, { P }>>
for Matrix<T, { N }, { M }>
where
T: Add<T, Output = T> + Mul<T, Output = T> + Clone,
Vector<T, { M }>: InnerSpace,
{
type Output = Matrix<<Vector<T, { M }> as VectorSpace>::Scalar, { N }, { P }>;
fn mul(self, rhs: Matrix<T, { M }, { P }>) -> Self::Output {
let mut mat = MaybeUninit::<
[Vector<<Vector<T, { M }> as VectorSpace>::Scalar, { N }>; P],
>::uninit();
let matp: *mut Vector<<Vector<T, { M }> as VectorSpace>::Scalar, { N }> =
unsafe { mem::transmute(&mut mat) };
for i in 0..P {
let mut column =
MaybeUninit::<[<Vector<T, { M }> as VectorSpace>::Scalar; N]>::uninit();
let columnp: *mut <Vector<T, { M }> as VectorSpace>::Scalar =
unsafe { mem::transmute(&mut column) };
for j in 0..N {
let mut row = MaybeUninit::<[T; M]>::uninit();
let rowp: *mut T = unsafe { mem::transmute(&mut row) };
for k in 0..M {
unsafe {
rowp.add(k).write(self.0[k].0[j].clone());
}
}
let row = Vector::<T, { M }>::from(unsafe { row.assume_init() });
unsafe {
columnp.add(j).write(row.dot(rhs.0[i].clone()));
}
}
let column = Vector::<<Vector<T, { M }> as VectorSpace>::Scalar, { N }>(unsafe {
column.assume_init()
});
unsafe {
matp.add(i).write(column);
}
}
Matrix::<<Vector<T, { M }> as VectorSpace>::Scalar, { N }, { P }>(unsafe {
mat.assume_init()
})
}
}
impl<T, const N: usize, const M: usize> Mul<Vector<T, { M }>> for Matrix<T, { N }, { M }>
where
T: Add<T, Output = T> + Mul<T, Output = T> + Clone,
Vector<T, { M }>: InnerSpace,
{
type Output = Vector<<Vector<T, { M }> as VectorSpace>::Scalar, { N }>;
fn mul(self, rhs: Vector<T, { M }>) -> Self::Output {
let mut column =
MaybeUninit::<[<Vector<T, { M }> as VectorSpace>::Scalar; N]>::uninit();
let columnp: *mut <Vector<T, { M }> as VectorSpace>::Scalar =
unsafe { mem::transmute(&mut column) };
for j in 0..N {
let mut row = MaybeUninit::<[T; M]>::uninit();
let rowp: *mut T = unsafe { mem::transmute(&mut row) };
for k in 0..M {
unsafe {
rowp.add(k).write(self.0[k].0[j].clone());
}
}
let row = Vector::<T, { M }>::from(unsafe { row.assume_init() });
unsafe {
columnp.add(j).write(row.dot(rhs.clone()));
}
}
Vector::<<Vector<T, { M }> as VectorSpace>::Scalar, { N }>(unsafe { column.assume_init() })
}
}
impl<T, const N: usize, const M: usize> Mul<T> for Matrix<T, { N }, { M }>
where
T: Mul<T, Output = T> + Clone,
{
type Output = Matrix<T, { N }, { M }>;
fn mul(self, scalar: T) -> Self::Output {
let mut mat = MaybeUninit::<[Vector<T, { N }>; M]>::uninit();
let matp: *mut Vector<T, { N }> = unsafe { mem::transmute(&mut mat) };
for i in 0..M {
unsafe {
matp.add(i).write(self.0[i].clone() * scalar.clone());
}
}
Matrix::<T, { N }, { M }>(unsafe { mat.assume_init() })
}
}
impl<const N: usize, const M: usize> Mul<Matrix<f32, { N }, { M }>> for f32 {
type Output = Matrix<f32, { N }, { M }>;
fn mul(self, mat: Matrix<f32, { N }, { M }>) -> Self::Output {
mat.map(|x| x * self)
}
}
impl<const N: usize, const M: usize> Mul<Matrix<f64, { N }, { M }>> for f64 {
type Output = Matrix<f64, { N }, { M }>;
fn mul(self, mat: Matrix<f64, { N }, { M }>) -> Self::Output {
mat.map(|x| x * self)
}
}
#[derive(Copy, Clone)]
pub struct Permutation<const N: usize> {
arr: [usize; N],
num_swaps: usize,
}
impl<const N: usize> fmt::Debug for Permutation<{ N }> {
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
write!(f, "[ ")?;
for i in 0..N {
write!(f, "{:?} ", self.arr[i])?;
}
write!(f, "] ")
}
}
impl<RHS, const N: usize> PartialEq<RHS> for Permutation<{ N }>
where
RHS: Deref<Target = [usize; N]>,
{
fn eq(&self, other: &RHS) -> bool {
for (a, b) in self.arr.iter().zip(other.deref().iter()) {
if !a.eq(b) {
return false;
}
}
true
}
}
impl<const N: usize> Deref for Permutation<{ N }> {
type Target = [usize; N];
fn deref(&self) -> &Self::Target {
&self.arr
}
}
impl<const N: usize> DerefMut for Permutation<{ N }> {
fn deref_mut(&mut self) -> &mut Self::Target {
&mut self.arr
}
}
impl<const N: usize> Permutation<{ N }> {
pub fn unit() -> Permutation<{ N }> {
let mut arr: [MaybeUninit<usize>; N] = MaybeUninit::uninit_array();
let arr = unsafe {
for i in 0..N {
arr[i] = MaybeUninit::new(i);
}
transmute_copy::<_, _>(&arr)
};
Permutation { arr, num_swaps: 0 }
}
pub fn swap(&mut self, a: usize, b: usize) {
self.num_swaps += 1;
self.arr.swap(a, b);
}
pub fn num_swaps(&self) -> usize {
self.num_swaps
}
}
impl<T, const N: usize> Mul<Vector<T, { N }>> for Permutation<{ N }>
where
T: Clone,
{
type Output = Vector<T, { N }>;
fn mul(self, rhs: Vector<T, { N }>) -> Self::Output {
Vector::from_iter((0..N).map(|i| rhs[self[i]].clone()))
}
}
#[derive(Copy, Clone, Debug)]
pub struct LU<T, const N: usize>(Permutation<{ N }>, Matrix<T, { N }, { N }>);
impl<T, const N: usize> Index<(usize, usize)> for LU<T, { N }> {
type Output = T;
fn index(&self, (row, column): (usize, usize)) -> &Self::Output {
&self.1[(row, column)]
}
}
impl<T, const N: usize> LU<T, { N }>
where
T: Clone
+ PartialEq
+ One
+ Zero
+ Product
+ Neg<Output = T>
+ Sub<T, Output = T>
+ Mul<T, Output = T>
+ Div<T, Output = T>,
{
pub fn p(&self) -> &Permutation<{ N }> {
&self.0
}
pub fn solve(&self, b: Vector<T, { N }>) -> Vector<T, { N }> {
let mut x = self.0.clone() * b;
for i in 0..N {
for k in 0..i {
x[i] = x[i].clone() - self[(i, k)].clone() * x[k].clone();
}
}
for i in (0..N).rev() {
for k in i + 1..N {
x[i] = x[i].clone() - self[(i, k)].clone() * x[k].clone();
}
x[i] = x[i].clone() / self[(i, i)].clone();
}
x
}
pub fn determinant(&self) -> T {
let det: T = self.1.diagonal().into_iter().product();
if self.0.num_swaps % 2 == 1 {
-det
} else {
det
}
}
pub fn invert(self) -> Matrix<T, { N }, { N }> {
Matrix::<T, { N }, { N }>::one()
.into_iter()
.map(|col| self.solve(col))
.collect()
}
}
#[cfg(feature = "rand")]
impl<T, const N: usize, const M: usize> Distribution<Matrix<T, { N }, { M }>> for Standard
where
Standard: Distribution<Vector<T, { N }>>,
{
fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Matrix<T, { N }, { M }> {
let mut rand = MaybeUninit::<[Vector<T, { N }>; { M }]>::uninit();
let randp: *mut Vector<T, { N }> = unsafe { mem::transmute(&mut rand) };
for i in 0..M {
unsafe {
randp.add(i).write(self.sample(rng));
}
}
Matrix::<T, { N }, { M }>(unsafe { rand.assume_init() })
}
}
#[cfg(feature = "serde")]
impl<T, const N: usize, const M: usize> Serialize for Matrix<T, { N }, { M }>
where
Vector<T, { N }>: Serialize,
{
fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
where
S: Serializer,
{
let mut seq = serializer.serialize_tuple(M)?;
for i in 0..M {
seq.serialize_element(&self.0[i])?;
}
seq.end()
}
}
#[cfg(feature = "serde")]
impl<'de, T, const N: usize, const M: usize> Deserialize<'de> for Matrix<T, { N }, { M }>
where
T: Deserialize<'de>,
{
fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>
where
D: Deserializer<'de>,
{
deserializer
.deserialize_tuple(N, ArrayVisitor::<[Vector<T, { N }>; { M }]>::new())
.map(Matrix)
}
}
macro_rules! into_mint_column_matrix {
($mint_name:ident, $rows:expr, $cols:expr $( , ($col_name:ident, $col_idx:expr ) )+) => {
#[cfg(feature = "mint")]
impl<T: Copy> Into<mint::$mint_name<T>> for Matrix<T, {$rows}, {$cols}> {
fn into(self) -> mint::$mint_name<T> {
mint::$mint_name {
$(
$col_name: self.0[$col_idx].into(),
)*
}
}
}
}
}
into_mint_column_matrix!(ColumnMatrix2, 2, 2, (x, 0), (y, 1));
into_mint_column_matrix!(ColumnMatrix3, 3, 3, (x, 0), (y, 1), (z, 2));
into_mint_column_matrix!(ColumnMatrix4, 4, 4, (x, 0), (y, 1), (z, 2), (w, 3));
into_mint_column_matrix!(ColumnMatrix2x3, 2, 3, (x, 0), (y, 1), (z, 2));
into_mint_column_matrix!(ColumnMatrix2x4, 2, 4, (x, 0), (y, 1), (z, 2), (w, 3));
into_mint_column_matrix!(ColumnMatrix3x2, 3, 2, (x, 0), (y, 1));
into_mint_column_matrix!(ColumnMatrix3x4, 3, 4, (x, 0), (y, 1), (z, 2), (w, 3));
into_mint_column_matrix!(ColumnMatrix4x2, 4, 2, (x, 0), (y, 1));
into_mint_column_matrix!(ColumnMatrix4x3, 4, 3, (x, 0), (y, 1), (z, 2));
macro_rules! from_mint_column_matrix {
($mint_name:ident, $rows:expr, $cols:expr, $($component:ident),+) => {
#[cfg(feature = "mint")]
impl<T> From<mint::$mint_name<T>> for Matrix<T, {$rows}, {$cols}> {
fn from(m: mint::$mint_name<T>) -> Self {
Self([
$(
Vector::<T, {$rows}>::from(m.$component),
)*
])
}
}
}
}
from_mint_column_matrix!(ColumnMatrix2, 2, 2, x, y);
from_mint_column_matrix!(ColumnMatrix3, 3, 3, x, y, z);
from_mint_column_matrix!(ColumnMatrix4, 4, 4, x, y, z, w);
from_mint_column_matrix!(ColumnMatrix2x3, 2, 3, x, y, z);
from_mint_column_matrix!(ColumnMatrix2x4, 2, 4, x, y, z, w);
from_mint_column_matrix!(ColumnMatrix3x2, 3, 2, x, y);
from_mint_column_matrix!(ColumnMatrix3x4, 3, 4, x, y, z, w);
from_mint_column_matrix!(ColumnMatrix4x2, 4, 2, x, y);
from_mint_column_matrix!(ColumnMatrix4x3, 4, 3, x, y, z);
macro_rules! into_mint_row_matrix {
($mint_name:ident, $rows:expr, $cols:expr $( , ($col_name:ident, $col_idx:expr ) )+) => {
#[cfg(feature = "mint")]
impl<T: Copy> Into<mint::$mint_name<T>> for Matrix<T, {$rows}, {$cols}> {
fn into(self) -> mint::$mint_name<T> {
let transposed = self.transpose();
mint::$mint_name {
$(
$col_name: transposed.0[$col_idx].into(),
)*
}
}
}
}
}
into_mint_row_matrix!(RowMatrix2, 2, 2, (x, 0), (y, 1));
into_mint_row_matrix!(RowMatrix3, 3, 3, (x, 0), (y, 1), (z, 2));
into_mint_row_matrix!(RowMatrix4, 4, 4, (x, 0), (y, 1), (z, 2), (w, 3));
into_mint_row_matrix!(RowMatrix2x3, 2, 3, (x, 0), (y, 1));
into_mint_row_matrix!(RowMatrix2x4, 2, 4, (x, 0), (y, 1));
into_mint_row_matrix!(RowMatrix3x2, 3, 2, (x, 0), (y, 1), (z, 2));
into_mint_row_matrix!(RowMatrix3x4, 3, 4, (x, 0), (y, 1), (z, 2));
into_mint_row_matrix!(RowMatrix4x2, 4, 2, (x, 0), (y, 1), (z, 2), (w, 3));
into_mint_row_matrix!(RowMatrix4x3, 4, 3, (x, 0), (y, 1), (z, 2), (w, 3));
macro_rules! from_mint_row_matrix {
($mint_name:ident, $rows:expr, $cols:expr, $($component:ident),+) => {
#[cfg(feature = "mint")]
impl<T> From<mint::$mint_name<T>> for Matrix<T, {$rows}, {$cols}> {
fn from(m: mint::$mint_name<T>) -> Self {
Matrix::<T, {$cols}, {$rows}>([
$(
Vector::<T, {$cols}>::from(m.$component),
)*
]).transpose()
}
}
}
}
from_mint_row_matrix!(RowMatrix2, 2, 2, x, y);
from_mint_row_matrix!(RowMatrix3, 3, 3, x, y, z);
from_mint_row_matrix!(RowMatrix4, 4, 4, x, y, z, w);
from_mint_row_matrix!(RowMatrix2x3, 2, 3, x, y);
from_mint_row_matrix!(RowMatrix2x4, 2, 4, x, y);
from_mint_row_matrix!(RowMatrix3x2, 3, 2, x, y, z);
from_mint_row_matrix!(RowMatrix3x4, 3, 4, x, y, z);
from_mint_row_matrix!(RowMatrix4x2, 4, 2, x, y, z, w);
from_mint_row_matrix!(RowMatrix4x3, 4, 3, x, y, z, w);