vector_victor/
math.rs

1// This Source Code Form is subject to the terms of the Mozilla Public
2// License, v. 2.0. If a copy of the MPL was not distributed with this
3// file, You can obtain one at http://mozilla.org/MPL/2.0/.
4
5use crate::{Matrix, Vector};
6use num_traits::real::Real;
7use num_traits::{Inv, Num, NumOps, One, Pow, Signed, Zero};
8use std::iter::{zip, Product, Sum};
9use std::ops::{Add, Mul};
10
11/// Operations for column vectors
12impl<T: Copy, const N: usize> Vector<T, N> {
13    /** Compute the dot product of two vectors, otherwise known as the scalar product.
14
15    This is the sum of the elementwise product, or in math terms
16
17    $vec(a) * vec(b) = sum_(i=1)^n a_i b_i = a_1 b_1 + a_2 b_2 + ... + a_n b_n$
18
19    for example, $\[\[1],\[2],\[3]] * \[\[4],\[5],\[6]] = (1 * 4) + (2 * 5) + (3 * 6) = 32$
20
21    For vectors in euclidean space, this has the property that it is equal to the magnitudes of
22    the vectors times the cosine of the angle between them.
23
24    $vec(a) * vec(b) = |vec(a)| |vec(b)| cos(theta)$
25
26    this also gives it the special property that the dot product of a vector and itself is the
27    square of its magnitude. You may recognize the 2D version as the
28    [pythagorean theorem](https://en.wikipedia.org/wiki/Pythagorean_theorem).
29
30    see [dot product](https://en.wikipedia.org/wiki/Dot_product) on Wikipedia for more
31    information. */
32    pub fn dot<R>(&self, rhs: &R) -> T
33    where
34        for<'s> &'s Self: Mul<&'s R, Output = Self>,
35        T: Sum<T>,
36    {
37        (self * rhs).elements().cloned().sum()
38    }
39
40    pub fn sqrmag(&self) -> T
41    where
42        for<'s> &'s Self: Mul<&'s Self, Output = Self>,
43        T: Sum<T>,
44    {
45        self.dot(self)
46    }
47
48    pub fn mag(&self) -> T
49    where
50        T: Sum<T> + Mul<T> + Real,
51    {
52        self.sqrmag().sqrt()
53    }
54
55    pub fn normalized(&self) -> Option<Self>
56    where
57        T: Sum<T> + Mul<T> + Real,
58    {
59        match self.mag() {
60            mag if mag.abs() < T::epsilon() => None,
61            mag => Some(self / mag),
62        }
63    }
64}
65
66/// Cross product operations for column vectors in $RR^3$
67impl<T: Copy> Vector<T, 3> {
68    pub fn cross_r<R: Copy>(&self, rhs: &Vector<R, 3>) -> Self
69    where
70        T: NumOps<R> + NumOps,
71    {
72        Self::vec([
73            (self[1] * rhs[2]) - (self[2] * rhs[1]),
74            (self[2] * rhs[0]) - (self[0] * rhs[2]),
75            (self[0] * rhs[1]) - (self[1] * rhs[0]),
76        ])
77    }
78
79    pub fn cross_l<R: Copy>(&self, rhs: &Vector<R, 3>) -> Vector<R, 3>
80    where
81        R: NumOps<T> + NumOps,
82    {
83        rhs.cross_r(self)
84    }
85}
86
87/// Operations for Matrices
88impl<T: Copy, const M: usize, const N: usize> Matrix<T, M, N> {
89    pub fn mmul<R: Copy, const P: usize>(&self, rhs: &Matrix<R, N, P>) -> Matrix<T, M, P>
90    where
91        T: Default + NumOps<R> + Sum,
92    {
93        let mut result: Matrix<T, M, P> = Default::default();
94
95        for (m, a) in self.rows().enumerate() {
96            for (n, b) in rhs.cols().enumerate() {
97                result[(m, n)] = a.dot(&b)
98            }
99        }
100
101        return result;
102    }
103
104    /// Computes the absolute value of each element of the matrix
105    pub fn abs(&self) -> Self
106    where
107        T: Signed + Default,
108    {
109        self.elements().map(|&x| x.abs()).collect()
110    }
111
112    /// Computes the sign of each element of the matrix
113    pub fn signum(&self) -> Self
114    where
115        T: Signed + Default,
116    {
117        self.elements().map(|&x| x.signum()).collect()
118    }
119
120    /// Raises every element to the power of rhs, where rhs is either a scalar or a matrix of exponents
121    pub fn pow<R, O>(self, rhs: R) -> O
122    where
123        Self: Pow<R, Output = O>,
124    {
125        Pow::pow(self, rhs)
126    }
127}
128
129// Sum up matrices
130impl<T: Copy, const M: usize, const N: usize> Sum for Matrix<T, M, N>
131where
132    Self: Zero + Add<Output = Self>,
133{
134    fn sum<I: Iterator<Item = Self>>(iter: I) -> Self {
135        iter.fold(Self::zero(), Self::add)
136    }
137}
138
139// Product of matrices
140impl<T: Copy, const M: usize, const N: usize> Product for Matrix<T, M, N>
141where
142    Self: One + Mul<Output = Self>,
143{
144    fn product<I: Iterator<Item = Self>>(iter: I) -> Self {
145        iter.fold(Self::one(), Self::mul)
146    }
147}
148
149/// Inverse trait. Note that this is the elementwise inverse, not the matrix inverse.
150/// For the inverse matrix see [`LUDecomposable::inv()`](crate::decompose::LUDecompose::inv())
151impl<T: Copy + Inv<Output = T> + Default, const M: usize, const N: usize> Inv for Matrix<T, M, N> {
152    type Output = Self;
153
154    fn inv(self) -> Self::Output {
155        self.elements().map(|t| t.inv()).collect()
156    }
157}
158
159/// Pow for $Matrix^{scalar}$
160impl<T, R, O, const M: usize, const N: usize> Pow<R> for Matrix<T, M, N>
161where
162    T: Copy + Pow<R, Output = O>,
163    R: Copy + Num,
164    O: Copy + Default,
165{
166    type Output = Matrix<O, M, N>;
167
168    fn pow(self, rhs: R) -> Self::Output {
169        self.elements().map(|&x| x.pow(rhs)).collect()
170    }
171}
172
173/// Pow for $Matrix^{Matrix}$
174impl<T, R, O, const M: usize, const N: usize> Pow<Matrix<R, M, N>> for Matrix<T, M, N>
175where
176    T: Copy + Pow<R, Output = O>,
177    R: Copy,
178    O: Copy + Default,
179{
180    type Output = Matrix<O, M, N>;
181
182    fn pow(self, rhs: Matrix<R, M, N>) -> Self::Output {
183        zip(self.elements(), rhs.elements())
184            .map(|(x, &r)| x.pow(r))
185            .collect()
186    }
187}