matrix/
complex.rs

1use std::{
2    fmt::{Debug, Display},
3    iter::Sum,
4    ops::{
5        Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign,
6    },
7};
8
9use crate::{
10    matrix::Transpose,
11    scalar::{Lerp, MulAdd, Scalar, Sqrt},
12    utils::EPSILON,
13    vector::Angle,
14    Dot, Matrix, Vector, V,
15};
16
17#[derive(Copy, Clone, Default, PartialEq, PartialOrd)]
18pub struct Complex {
19    pub x: f64,
20    pub y: f64,
21}
22
23pub trait Conj {
24    fn conj(&self) -> Complex;
25}
26
27#[macro_export]
28macro_rules! C {
29    ($r:expr, $i:expr) => {
30        Complex::from([$r, $i])
31    };
32}
33
34impl Debug for Complex {
35    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
36        write!(f, "({} + {}i)", self.x, self.y)
37    }
38}
39
40impl Display for Complex {
41    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
42        write!(f, "({} + {}i)", self.x, self.y)
43    }
44}
45
46impl From<[f64; 2]> for Complex {
47    fn from(value: [f64; 2]) -> Self {
48        Complex {
49            x: value[0],
50            y: value[1],
51        }
52    }
53}
54
55impl Scalar for Complex {
56    type AbsOutput = f64;
57
58    fn abs(&self) -> Self::AbsOutput {
59        (self.x.powi(2) + self.y.powi(2)).powf(0.5)
60    }
61
62    fn inv(self) -> Self {
63        Complex {
64            x: self.x / (self.x.powi(2) + self.y.powi(2)),
65            y: -self.y / (self.x.powi(2) + self.y.powi(2)),
66        }
67    }
68
69    fn one() -> Self {
70        C!(1., 0.)
71    }
72
73    type TanOutput = Complex;
74    fn tan(self) -> Self::TanOutput {
75        Complex {
76            x: f64::sin(2. * self.x)
77                / (f64::cos(2. * self.x) + f64::cosh(2. * self.y)),
78            y: f64::sinh(2. * self.y)
79                / (f64::cos(2. * self.x) + f64::cosh(2. * self.y)),
80        }
81    }
82
83    type SinOutput = Complex;
84    fn sin(self) -> Self::SinOutput {
85        Complex {
86            x: f64::sin(self.x) * f64::cosh(self.y),
87            y: f64::cos(self.x) * f64::sinh(self.y),
88        }
89    }
90
91    type CosOutput = Complex;
92    fn cos(self) -> Self::CosOutput {
93        Complex {
94            x: f64::cos(self.x) * f64::cosh(self.y),
95            y: f64::sin(self.x) * f64::sinh(self.y),
96        }
97    }
98
99    fn is_non_zero(&self) -> bool {
100        self.x.abs() > EPSILON || self.y.abs() > EPSILON
101    }
102}
103
104impl Sqrt for Complex {
105    fn sqrt(self: Self) -> Self {
106        Complex {
107            x: ((self.abs() + self.x) / 2.).sqrt(),
108            y: self.y / self.y.abs() * ((self.abs() - self.x) / 2.).sqrt(),
109        }
110    }
111}
112
113impl Lerp for Complex {
114    fn lerp(u: Self, v: Self, t: f32) -> Self {
115        match t {
116            0. => u,
117            1. => v,
118            p => (v - u) * p as f64 + u,
119        }
120    }
121}
122
123impl Sum for Complex {
124    fn sum<I: Iterator<Item = Self>>(iter: I) -> Self {
125        iter.fold(Complex::default(), |a, b| a + b)
126    }
127}
128
129impl Neg for Complex {
130    type Output = Self;
131    fn neg(self) -> Self::Output {
132        Complex {
133            x: self.x.neg(),
134            y: self.y.neg(),
135        }
136    }
137}
138
139impl Add for Complex {
140    type Output = Self;
141    fn add(self, rhs: Self) -> Self::Output {
142        Complex {
143            x: self.x + rhs.x,
144            y: self.y + rhs.y,
145        }
146    }
147}
148
149impl AddAssign for Complex {
150    fn add_assign(&mut self, rhs: Self) {
151        self.x += rhs.x;
152        self.y += rhs.y;
153    }
154}
155
156impl AddAssign<&Complex> for Complex {
157    fn add_assign(&mut self, rhs: &Complex) {
158        self.x += rhs.x;
159        self.y += rhs.y;
160    }
161}
162
163impl Sub for Complex {
164    type Output = Self;
165    fn sub(self, rhs: Self) -> Self::Output {
166        Complex {
167            x: self.x - rhs.x,
168            y: self.y - rhs.y,
169        }
170    }
171}
172
173impl SubAssign for Complex {
174    fn sub_assign(&mut self, rhs: Self) {
175        self.x -= rhs.x;
176        self.y -= rhs.y;
177    }
178}
179
180impl SubAssign<&Complex> for Complex {
181    fn sub_assign(&mut self, rhs: &Complex) {
182        self.x -= rhs.x;
183        self.y -= rhs.y;
184    }
185}
186
187impl Mul for Complex {
188    type Output = Self;
189    fn mul(self, rhs: Self) -> Self::Output {
190        Complex {
191            x: self.x * rhs.x - self.y * rhs.y,
192            y: self.x * rhs.y + self.y * rhs.x,
193        }
194    }
195}
196
197impl MulAssign for Complex {
198    fn mul_assign(&mut self, rhs: Self) {
199        let x = self.x * rhs.x - self.y * rhs.y;
200        self.y = self.x * rhs.y + self.y * rhs.x;
201        self.x = x;
202    }
203}
204
205impl MulAssign<&Complex> for Complex {
206    fn mul_assign(&mut self, rhs: &Complex) {
207        let x = self.x * rhs.x - self.y * rhs.y;
208        self.y = self.x * rhs.y + self.y * rhs.x;
209        self.x = x;
210    }
211}
212
213impl MulAdd<Complex, Complex> for Complex {
214    fn mul_add(self, a: &Self, b: &Self) -> Self {
215        Complex {
216            x: self.x.mul_add(a.x, self.y.mul_add(-a.y, b.x)),
217            y: self.x.mul_add(a.y, self.y.mul_add(a.x, b.y)),
218        }
219    }
220}
221
222impl Mul<f64> for Complex {
223    type Output = Complex;
224    fn mul(self, rhs: f64) -> Self::Output {
225        Complex {
226            x: self.x * rhs,
227            y: self.y * rhs,
228        }
229    }
230}
231
232impl MulAssign<f64> for Complex {
233    fn mul_assign(&mut self, rhs: f64) {
234        self.x *= rhs;
235        self.y *= rhs;
236    }
237}
238
239impl MulAdd<f64, Complex> for Complex {
240    fn mul_add(self, a: &f64, b: &Self) -> Self {
241        Complex {
242            x: self.x.mul_add(*a, b.x),
243            y: self.y.mul_add(*a, b.y),
244        }
245    }
246}
247
248impl Div for Complex {
249    type Output = Self;
250    fn div(self, rhs: Self) -> Self {
251        Complex {
252            x: (self.x * rhs.x + self.y * rhs.y)
253                / (rhs.x.powi(2) + rhs.y.powi(2)),
254            y: (self.y * rhs.x - self.x * rhs.y)
255                / (rhs.x.powi(2) + rhs.y.powi(2)),
256        }
257    }
258}
259
260impl DivAssign<&Complex> for Complex {
261    fn div_assign(&mut self, rhs: &Complex) {
262        let x =
263            (self.x * rhs.x + self.y * rhs.y) / (rhs.x.powi(2) + rhs.y.powi(2));
264        self.y =
265            (self.y * rhs.x - self.x * rhs.y) / (rhs.x.powi(2) + rhs.y.powi(2));
266        self.x = x;
267    }
268}
269
270impl DivAssign for Complex {
271    fn div_assign(&mut self, rhs: Self) {
272        let x =
273            (self.x * rhs.x + self.y * rhs.y) / (rhs.x.powi(2) + rhs.y.powi(2));
274        self.y =
275            (self.y * rhs.x - self.x * rhs.y) / (rhs.x.powi(2) + rhs.y.powi(2));
276        self.x = x;
277    }
278}
279
280impl Dot<Complex> for Vector<Complex> {
281    fn dot(&self, v: &Vector<Complex>) -> Complex {
282        assert_eq!(v.size(), self.size(), "vectors must be the same size");
283
284        self * v
285    }
286}
287
288impl Mul<&Vector<Complex>> for &Vector<Complex> {
289    type Output = Complex;
290
291    fn mul(self, rhs: &Vector<Complex>) -> Self::Output {
292        self * &rhs._d
293    }
294}
295
296impl Mul<&Vec<Complex>> for &Vector<Complex> {
297    type Output = Complex;
298
299    fn mul(self, rhs: &Vec<Complex>) -> Self::Output {
300        self * &rhs[..]
301    }
302}
303
304impl Conj for Complex {
305    fn conj(&self) -> Complex {
306        Complex {
307            x: self.x,
308            y: -self.y,
309        }
310    }
311}
312
313impl Mul<&[Complex]> for &Vector<Complex> {
314    type Output = Complex;
315
316    fn mul(self, rhs: &[Complex]) -> Self::Output {
317        let mut sum = Complex::default();
318
319        for (a, b) in self._d.iter().zip(rhs) {
320            sum = a.mul_add(&b.conj(), &sum);
321        }
322
323        sum
324    }
325}
326
327impl Dot<Complex> for [Complex] {
328    fn dot(&self, v: &Vector<Complex>) -> Complex {
329        assert_eq!(v.size(), self.len(), "vectors must be the same size");
330
331        self * v
332    }
333}
334
335impl Mul<&Vector<Complex>> for &[Complex] {
336    type Output = Complex;
337
338    fn mul(self, rhs: &Vector<Complex>) -> Self::Output {
339        let mut sum = Complex::default();
340
341        for (a, b) in self.iter().zip(rhs._d.iter()) {
342            sum = a.mul_add(&b.conj(), &sum);
343        }
344
345        sum
346    }
347}
348impl MulAdd<Complex, Vector<Complex>> for Vector<Complex> {
349    fn mul_add(self, a: &Complex, b: &Vector<Complex>) -> Self {
350        assert!(self.size() == b.size(), "vectors must be the same size");
351
352        let mut vec = Vec::with_capacity(self.size());
353
354        for i in 0..self.size() {
355            vec.push(self[i].mul_add(a, &b[i]))
356        }
357
358        V!(vec)
359    }
360}
361
362impl MulAdd<f64, Vector<Complex>> for Vector<Complex> {
363    fn mul_add(self, a: &f64, b: &Vector<Complex>) -> Self {
364        assert!(self.size() == b.size(), "vectors must be the same size");
365
366        let mut vec = Vec::with_capacity(self.size());
367
368        for i in 0..self.size() {
369            vec.push(self[i].mul_add(a, &b[i]))
370        }
371
372        V!(vec)
373    }
374}
375
376impl Angle for Vector<Complex> {
377    type Output = f64;
378    fn angle_cos(u: &Vector<Complex>, v: &Vector<Complex>) -> Self::Output {
379        u.dot(v).x / (u.norm() * v.norm())
380    }
381}
382
383impl Transpose<Complex> for Matrix<Complex> {
384    fn transpose(&self) -> Matrix<Complex> {
385        let mut vec = Vec::with_capacity(self.rows * self.cols);
386        for i in 0..self.cols {
387            for j in 0..self.rows {
388                vec.push(self[j][i].conj());
389            }
390        }
391
392        Matrix {
393            rows: self.cols,
394            cols: self.rows,
395            _d: vec,
396        }
397    }
398}
399
400impl Lerp for Vector<Complex> {
401    fn lerp(u: Self, v: Self, t: f32) -> Self {
402        match t {
403            0. => u,
404            1. => v,
405            p => {
406                let mut vec = Vec::with_capacity(u.size());
407
408                for i in 0..u.size() {
409                    vec.push((v[i] - u[i]).mul_add(&(p as f64), &u[i]))
410                }
411
412                V!(vec)
413            }
414        }
415    }
416}
417
418impl Lerp for Matrix<Complex> {
419    fn lerp(u: Self, v: Self, t: f32) -> Self {
420        match t {
421            0. => u,
422            1. => v,
423            p => {
424                let mut vec = Vec::with_capacity(u._d.len());
425
426                for i in 0..u._d.len() {
427                    vec.push((v._d[i] - u._d[i]).mul_add(&(p as f64), &u._d[i]))
428                }
429
430                Matrix {
431                    _d: vec,
432                    cols: u.cols,
433                    rows: u.rows,
434                }
435            }
436        }
437    }
438}
matrix/complex.rs

matrix/
complex.rs
