abels_complex/complex/f32/
rectangular.rs

1use core::f32::consts::{LN_2, LN_10};
2use core::fmt;
3use core::ops::*;
4use core::write;
5#[cfg(feature = "libm")]
6use num_traits::real::Real;
7
8type Polar = crate::complex::polar::ComplexPolar<FT>;
9type Rectangular = crate::complex::rectangular::Complex<FT>;
10type FT = f32;
11
12impl Rectangular {
13    pub const ZERO: Self = Self::new(0.0, 0.0);
14    pub const ONE: Self = Self::new(1.0, 0.0);
15    pub const NEG_ONE: Self = Self::new(-1.0, 0.0);
16    pub const I: Self = Self::new(0.0, 1.0);
17    pub const NEG_I: Self = Self::new(0.0, -1.0);
18
19    /// Computes the conjugate.
20    pub const fn conjugate(self) -> Self {
21        Self::new(self.re, -self.im)
22    }
23
24    /// Computes the absolute value.
25    pub fn abs(self) -> FT {
26        self.abs_sq().sqrt()
27    }
28
29    /// Computes the squared absolute value.
30    ///
31    /// This is faster than `abs()` as it avoids a square root operation.
32    pub fn abs_sq(self) -> FT {
33        self.re * self.re + self.im * self.im
34    }
35
36    /// Computes the argument in the range `(-π, +π]`.
37    pub fn arg(self) -> FT {
38        self.im.atan2(self.re)
39    }
40
41    /// Computes the reciprocal.
42    pub fn recip(self) -> Self {
43        self.conjugate() / self.abs_sq()
44    }
45
46    /// Computes the principle square root.
47    pub fn sqrt(self) -> Self {
48        let abs = self.abs();
49        Self::new(
50            (0.5 * (abs + self.re)).sqrt(),
51            (0.5 * (abs - self.re)).sqrt().copysign(self.im),
52        )
53    }
54
55    /// Convert to polar form.
56    pub fn to_polar(self) -> Polar {
57        Polar::new(self.abs(), self.arg())
58    }
59
60    /// Computes `e^self` where `e` is the base of the natural logarithm.
61    pub fn exp(self) -> Polar {
62        Polar::new(self.re.exp(), self.im)
63    }
64
65    /// Computes the principle natural logarithm.
66    pub fn ln(self) -> Self {
67        self.to_polar().ln()
68    }
69
70    /// Computes the principle logarithm in base 2.
71    pub fn log2(self) -> Self {
72        self.ln() / LN_2
73    }
74
75    /// Computes the principle logarithm in base 10.
76    pub fn log10(self) -> Self {
77        self.ln() / LN_10
78    }
79
80    /// Raises `self` to an integer power.
81    pub fn powi(self, n: i32) -> Polar {
82        self.to_polar().powi(n)
83    }
84
85    /// Raises `self` to a floating point power.
86    pub fn powf(self, x: FT) -> Polar {
87        self.to_polar().powf(x)
88    }
89
90    /// Computes the euclidian distance between two points.
91    pub fn distance(self, other: Self) -> FT {
92        (self - other).abs()
93    }
94
95    /// Computes the squared euclidian distance between two points.
96    pub fn distance_squared(self, other: Self) -> FT {
97        (self - other).abs_sq()
98    }
99
100    /// Casts to a glam::Vec2.
101    #[cfg(feature = "glam")]
102    pub fn as_vec2(self) -> glam::Vec2 {
103        glam::vec2(self.re, self.im)
104    }
105}
106
107impl Add for Rectangular {
108    type Output = Self;
109    fn add(mut self, other: Self) -> Self::Output {
110        self += other;
111        self
112    }
113}
114
115impl Add<FT> for Rectangular {
116    type Output = Self;
117    fn add(mut self, re: FT) -> Self::Output {
118        self += re;
119        self
120    }
121}
122
123impl Add<Rectangular> for FT {
124    type Output = Rectangular;
125    fn add(self, mut z: Self::Output) -> Self::Output {
126        z += self;
127        z
128    }
129}
130
131impl AddAssign for Rectangular {
132    fn add_assign(&mut self, other: Self) {
133        self.re += other.re;
134        self.im += other.im;
135    }
136}
137
138impl AddAssign<FT> for Rectangular {
139    fn add_assign(&mut self, re: FT) {
140        self.re += re;
141    }
142}
143
144impl Sub for Rectangular {
145    type Output = Self;
146    fn sub(mut self, other: Self) -> Self::Output {
147        self -= other;
148        self
149    }
150}
151
152impl Sub<FT> for Rectangular {
153    type Output = Self;
154    fn sub(mut self, re: FT) -> Self::Output {
155        self -= re;
156        self
157    }
158}
159
160impl Sub<Rectangular> for FT {
161    type Output = Rectangular;
162    fn sub(self, z: Self::Output) -> Self::Output {
163        Rectangular::new(self - z.re, -z.im)
164    }
165}
166
167impl SubAssign for Rectangular {
168    fn sub_assign(&mut self, other: Self) {
169        self.re -= other.re;
170        self.im -= other.im;
171    }
172}
173
174impl SubAssign<FT> for Rectangular {
175    fn sub_assign(&mut self, re: FT) {
176        self.re -= re;
177    }
178}
179
180impl Mul for Rectangular {
181    type Output = Self;
182    fn mul(mut self, other: Self) -> Self {
183        self *= other;
184        self
185    }
186}
187
188impl Mul<FT> for Rectangular {
189    type Output = Self;
190    fn mul(mut self, re: FT) -> Self {
191        self *= re;
192        self
193    }
194}
195
196impl Mul<Rectangular> for FT {
197    type Output = Rectangular;
198    fn mul(self, mut other: Self::Output) -> Self::Output {
199        other *= self;
200        other
201    }
202}
203
204impl MulAssign for Rectangular {
205    fn mul_assign(&mut self, other: Self) {
206        let re = self.re * other.re - self.im * other.im;
207        self.im = self.re * other.im + self.im * other.re;
208        self.re = re;
209    }
210}
211
212impl MulAssign<FT> for Rectangular {
213    fn mul_assign(&mut self, re: FT) {
214        self.re *= re;
215        self.im *= re;
216    }
217}
218
219impl Div for Rectangular {
220    type Output = Self;
221    fn div(self, other: Self) -> Self::Output {
222        self * other.recip()
223    }
224}
225
226impl Div<FT> for Rectangular {
227    type Output = Self;
228    fn div(mut self, re: FT) -> Self::Output {
229        self /= re;
230        self
231    }
232}
233
234impl Div<Rectangular> for FT {
235    type Output = Rectangular;
236    fn div(self, other: Self::Output) -> Self::Output {
237        self * other.recip()
238    }
239}
240
241impl DivAssign for Rectangular {
242    fn div_assign(&mut self, other: Self) {
243        *self = *self / other;
244    }
245}
246
247impl DivAssign<FT> for Rectangular {
248    fn div_assign(&mut self, re: FT) {
249        self.re /= re;
250        self.im /= re;
251    }
252}
253
254impl Neg for Rectangular {
255    type Output = Self;
256    fn neg(self) -> Self::Output {
257        Self::new(-self.re, -self.im)
258    }
259}
260
261impl fmt::Display for Rectangular {
262    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
263        fn fmt_x(f: &mut fmt::Formatter, x: FT, sign_plus: bool) -> fmt::Result {
264            match (f.precision(), sign_plus) {
265                (None, false) => write!(f, "{}", x),
266                (None, true) => write!(f, "{:+}", x),
267                (Some(p), false) => write!(f, "{:.*}", p, x),
268                (Some(p), true) => write!(f, "{:+.*}", p, x),
269            }
270        }
271        match (self.re, self.im, f.sign_plus()) {
272            (re, 0.0, sp) => fmt_x(f, re, sp),
273            (0.0, 1.0, false) => write!(f, "i"),
274            (0.0, 1.0, true) => write!(f, "+i"),
275            (0.0, -1.0, _) => write!(f, "-i"),
276            (0.0, im, sp) => {
277                fmt_x(f, im, sp)?;
278                write!(f, "i")
279            }
280            (re, 1.0, sp) => {
281                fmt_x(f, re, sp)?;
282                write!(f, "+i")
283            }
284            (re, -1.0, sp) => {
285                fmt_x(f, re, sp)?;
286                write!(f, "-i")
287            }
288            (re, im, sp) => {
289                fmt_x(f, re, sp)?;
290                fmt_x(f, im, true)?;
291                write!(f, "i")
292            }
293        }
294    }
295}
296
297#[cfg(feature = "rand")]
298impl rand::distr::Distribution<Rectangular> for rand::distr::StandardUniform {
299    fn sample<R: rand::Rng + ?Sized>(&self, rng: &mut R) -> Rectangular {
300        rng.sample::<Polar, _>(self).to_rectangular()
301    }
302}
303
304#[cfg(feature = "approx")]
305use approx::{AbsDiffEq, RelativeEq, UlpsEq};
306
307#[cfg(feature = "approx")]
308impl AbsDiffEq for Rectangular {
309    type Epsilon = <FT as AbsDiffEq>::Epsilon;
310    fn default_epsilon() -> Self::Epsilon {
311        FT::default_epsilon()
312    }
313    fn abs_diff_eq(&self, other: &Self, epsilon: Self::Epsilon) -> bool {
314        FT::abs_diff_eq(&self.re, &other.re, epsilon)
315            && FT::abs_diff_eq(&self.im, &other.im, epsilon)
316    }
317}
318
319#[cfg(feature = "approx")]
320impl RelativeEq for Rectangular {
321    fn default_max_relative() -> Self::Epsilon {
322        FT::default_max_relative()
323    }
324    fn relative_eq(
325        &self,
326        other: &Self,
327        epsilon: Self::Epsilon,
328        max_relative: Self::Epsilon,
329    ) -> bool {
330        FT::relative_eq(&self.re, &other.re, epsilon, max_relative)
331            && FT::relative_eq(&self.im, &other.im, epsilon, max_relative)
332    }
333}
334
335#[cfg(feature = "approx")]
336impl UlpsEq for Rectangular {
337    fn default_max_ulps() -> u32 {
338        FT::default_max_ulps()
339    }
340    fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool {
341        FT::ulps_eq(&self.re, &other.re, epsilon, max_ulps)
342            && FT::ulps_eq(&self.im, &other.im, epsilon, max_ulps)
343    }
344}
345
346impl From<FT> for Rectangular {
347    fn from(value: FT) -> Self {
348        Self::new(value, 0.0)
349    }
350}
351
352#[cfg(feature = "glam")]
353impl From<glam::Vec2> for Rectangular {
354    fn from(v: glam::Vec2) -> Self {
355        Self::new(v.x, v.y)
356    }
357}
358
359#[cfg(feature = "glam")]
360impl From<Rectangular> for glam::Vec2 {
361    fn from(z: Rectangular) -> Self {
362        z.as_vec2()
363    }
364}
365
366#[cfg(test)]
367mod tests {
368    use super::*;
369    use approx::*;
370    use core::f32::consts::{E, FRAC_PI_2, PI, SQRT_2};
371    use core::iter::Iterator;
372    use rand::{
373        Rng, SeedableRng,
374        distr::{Distribution, StandardUniform},
375        rngs::StdRng,
376    };
377
378    const NUM_SAMPLES: usize = 100;
379
380    fn random_samples<T>() -> impl Iterator<Item = T>
381    where
382        StandardUniform: Distribution<T>,
383    {
384        StdRng::seed_from_u64(21)
385            .sample_iter(StandardUniform)
386            .take(NUM_SAMPLES)
387    }
388
389    #[test]
390    fn addition() {
391        for z0 in random_samples::<Rectangular>() {
392            for z1 in random_samples::<Rectangular>() {
393                let z = z0 + z1;
394                assert_eq!(z.re, z0.re + z1.re);
395                assert_eq!(z.im, z0.im + z1.im);
396
397                let z = z0 + z1.re;
398                assert_eq!(z.re, z0.re + z1.re);
399                assert_eq!(z.im, z0.im);
400
401                let z = z0.re + z1;
402                assert_eq!(z.re, z0.re + z1.re);
403                assert_eq!(z.im, z1.im);
404
405                let mut z = z0;
406                z += z1;
407                assert_eq!(z, z0 + z1);
408
409                let mut z = z0;
410                z += z1.re;
411                assert_eq!(z, z0 + z1.re);
412            }
413            assert_eq!(z0 + Rectangular::ZERO, z0);
414        }
415    }
416
417    #[test]
418    fn subtraction() {
419        for z0 in random_samples::<Rectangular>() {
420            for z1 in random_samples::<Rectangular>() {
421                let z = z0 - z1;
422                assert_eq!(z.re, z0.re - z1.re);
423                assert_eq!(z.im, z0.im - z1.im);
424
425                let z = z0 - z1.re;
426                assert_eq!(z.re, z0.re - z1.re);
427                assert_eq!(z.im, z0.im);
428
429                let z = z0.re - z1;
430                assert_eq!(z.re, z0.re - z1.re);
431                assert_eq!(z.im, -z1.im);
432
433                let mut z = z0;
434                z -= z1;
435                assert_eq!(z, z0 - z1);
436
437                let mut z = z0;
438                z -= z1.re;
439                assert_eq!(z, z0 - z1.re);
440            }
441            assert_eq!(z0 - z0, Rectangular::ZERO);
442            assert_eq!(z0 - Rectangular::ZERO, z0);
443        }
444    }
445
446    #[test]
447    fn multiplication() {
448        for z0 in random_samples::<Rectangular>() {
449            for z1 in random_samples::<Rectangular>() {
450                let z = z0 * z1;
451                assert_ulps_eq!(z.abs(), z0.abs() * z1.abs());
452                assert_ulps_eq!(
453                    z.arg().sin(),
454                    (z0.arg() + z1.arg()).sin(),
455                    epsilon = 4.0 * FT::EPSILON
456                );
457
458                let z = z0 * z1.re;
459                assert_eq!(z.re, z0.re * z1.re);
460                assert_eq!(z.im, z0.im * z1.re);
461
462                let z = z0.re * z1;
463                assert_eq!(z.re, z0.re * z1.re);
464                assert_eq!(z.im, z0.re * z1.im);
465
466                let mut z = z0;
467                z *= z1;
468                assert_eq!(z, z0 * z1);
469
470                let mut z = z0;
471                z *= z1.re;
472                assert_eq!(z, z0 * z1.re);
473            }
474            assert_eq!(z0 * Rectangular::ONE, z0);
475            assert_eq!(z0 * Rectangular::ZERO, Rectangular::ZERO);
476            assert_eq!(z0 * 0.0, Rectangular::ZERO);
477        }
478    }
479
480    #[test]
481    fn division() {
482        for z0 in random_samples::<Rectangular>() {
483            for z1 in random_samples::<Rectangular>() {
484                let z = z0 / z1;
485                assert_relative_eq!(
486                    z.abs(),
487                    z0.abs() / z1.abs(),
488                    max_relative = 3.0 * FT::EPSILON
489                );
490                assert_ulps_eq!(
491                    z.arg().sin(),
492                    (z0.arg() - z1.arg()).sin(),
493                    epsilon = 4.0 * FT::EPSILON
494                );
495
496                let z = z0 / z1.re;
497                assert_eq!(z.re, z0.re / z1.re);
498                assert_eq!(z.im, z0.im / z1.re);
499
500                let z = z0.re / z1;
501                assert_ulps_eq!(z.abs(), z0.re.abs() / z1.abs());
502                assert_ulps_eq!(
503                    z.arg().sin(),
504                    (-z0.re.signum() * z1.arg()).sin(),
505                    epsilon = 2.0 * FT::EPSILON
506                );
507
508                let mut z = z0;
509                z /= z1;
510                assert_eq!(z, z0 / z1);
511
512                let mut z = z0;
513                z /= z1.re;
514                assert_eq!(z, z0 / z1.re);
515            }
516            assert_ulps_eq!(z0 / z0, Rectangular::ONE);
517            assert_eq!(Rectangular::ZERO / z0, Rectangular::ZERO);
518        }
519    }
520
521    #[test]
522    fn negation() {
523        for z in random_samples::<Rectangular>() {
524            assert_eq!(-z, Rectangular::new(-z.re, -z.im));
525        }
526        assert_eq!(-Rectangular::ONE, Rectangular::NEG_ONE);
527        assert_eq!(-Rectangular::I, Rectangular::NEG_I);
528        assert_eq!(-Rectangular::NEG_ONE, Rectangular::ONE);
529        assert_eq!(-Rectangular::NEG_I, Rectangular::I);
530    }
531
532    #[test]
533    fn reciprocal() {
534        for z in random_samples::<Rectangular>() {
535            assert_eq!(z.recip(), 1.0 / z);
536            assert_ulps_eq!(z * z.recip(), Rectangular::ONE);
537        }
538        assert_eq!(Rectangular::ONE.recip(), Rectangular::ONE);
539        assert_eq!(Rectangular::I.recip(), Rectangular::NEG_I);
540        assert_eq!(Rectangular::NEG_ONE.recip(), Rectangular::NEG_ONE);
541        assert_eq!(Rectangular::NEG_I.recip(), Rectangular::I);
542    }
543
544    #[test]
545    fn sqrt() {
546        for z in random_samples::<Rectangular>() {
547            assert_ulps_eq!(z.sqrt().abs(), z.abs().sqrt());
548            assert_ulps_eq!(
549                z.sqrt().arg(),
550                z.arg() / 2.0,
551                epsilon = 1400.0 * FT::EPSILON
552            );
553        }
554        assert_eq!(Rectangular::ONE.sqrt(), Rectangular::ONE);
555        assert_eq!(Rectangular::NEG_ONE.sqrt(), Rectangular::I);
556        assert_eq!(
557            Rectangular::new(0.0, 2.0).sqrt(),
558            Rectangular::new(1.0, 1.0)
559        );
560        assert_eq!(
561            Rectangular::new(0.0, -2.0).sqrt(),
562            Rectangular::new(1.0, -1.0)
563        );
564    }
565
566    #[test]
567    fn abs() {
568        for z in random_samples::<Rectangular>() {
569            assert_ulps_eq!(z.abs_sq(), z.abs() * z.abs());
570            assert_eq!(z.abs_sq(), z.re * z.re + z.im * z.im);
571        }
572        assert_eq!(Rectangular::ONE.abs(), 1.0);
573        assert_eq!(Rectangular::I.abs(), 1.0);
574        assert_eq!(Rectangular::NEG_ONE.abs(), 1.0);
575        assert_eq!(Rectangular::NEG_I.abs(), 1.0);
576        assert_eq!(Rectangular::new(1.0, 1.0).abs(), SQRT_2);
577        assert_eq!(Rectangular::new(-1.0, 1.0).abs(), SQRT_2);
578        assert_eq!(Rectangular::new(-1.0, -1.0).abs(), SQRT_2);
579        assert_eq!(Rectangular::new(1.0, -1.0).abs(), SQRT_2);
580    }
581
582    #[test]
583    fn conjugate() {
584        for z in random_samples::<Rectangular>() {
585            assert_eq!(z.conjugate().re, z.re);
586            assert_eq!(z.conjugate().im, -z.im);
587            assert_eq!(z.conjugate().conjugate(), z);
588        }
589        assert_eq!(Rectangular::ONE.conjugate(), Rectangular::ONE);
590        assert_eq!(Rectangular::I.conjugate(), Rectangular::NEG_I);
591        assert_eq!(Rectangular::NEG_ONE.conjugate(), Rectangular::NEG_ONE);
592        assert_eq!(Rectangular::NEG_I.conjugate(), Rectangular::I);
593    }
594
595    #[test]
596    fn arg() {
597        assert_eq!(Rectangular::ONE.arg(), 0.0);
598        assert_eq!(Rectangular::I.arg(), FRAC_PI_2);
599        assert_eq!(Rectangular::NEG_ONE.arg(), PI);
600        assert_eq!(Rectangular::NEG_I.arg(), -FRAC_PI_2);
601    }
602
603    #[test]
604    fn exp() {
605        for z in random_samples::<Rectangular>() {
606            assert_eq!(z.exp().abs, z.re.exp());
607            assert_eq!(z.exp().arg, z.im);
608            assert_ulps_eq!(z.exp().ln(), z);
609        }
610        assert_eq!(Rectangular::ONE.exp(), Polar::new(E, 0.0));
611        assert_eq!(Rectangular::I.exp(), Polar::new(1.0, 1.0));
612        assert_eq!(Rectangular::NEG_ONE.exp(), Polar::new(E.recip(), 0.0));
613        assert_eq!(Rectangular::NEG_I.exp(), Polar::new(1.0, -1.0));
614    }
615
616    #[test]
617    fn log() {
618        for z in random_samples::<Rectangular>() {
619            assert_eq!(z.ln().re, z.abs().ln());
620            assert_eq!(z.ln().im, z.arg());
621            assert_ulps_eq!(z.ln().exp(), z.to_polar());
622        }
623        assert_eq!(Rectangular::ONE.ln(), Rectangular::ZERO);
624        assert_eq!(Rectangular::I.ln(), Rectangular::I * FRAC_PI_2);
625        assert_eq!(Rectangular::NEG_ONE.ln(), Rectangular::I * PI);
626        assert_eq!(Rectangular::NEG_I.ln(), Rectangular::I * -FRAC_PI_2);
627
628        assert_ulps_eq!(Rectangular::new(E, 0.0).ln(), Rectangular::ONE);
629        assert_ulps_eq!(Rectangular::new(2.0, 0.0).log2(), Rectangular::ONE);
630        assert_ulps_eq!(Rectangular::new(10.0, 0.0).log10(), Rectangular::ONE);
631    }
632
633    #[test]
634    fn powi() {
635        for z in random_samples::<Rectangular>() {
636            assert_eq!(z.powi(0), Polar::ONE);
637            assert_eq!(z.powi(1), z.to_polar());
638            for n in random_samples::<i32>() {
639                assert_eq!(z.powi(n).abs, z.abs().powi(n));
640                assert_eq!(z.powi(n).arg, z.arg() * n as FT);
641            }
642        }
643        for n in random_samples::<i32>() {
644            assert_eq!(Rectangular::ZERO.powi(n.abs()), Polar::ZERO);
645            assert_eq!(Rectangular::ONE.powi(n), Polar::ONE);
646        }
647    }
648
649    #[test]
650    fn powf() {
651        for z in random_samples::<Rectangular>() {
652            assert_eq!(z.powf(0.0), Polar::ONE);
653            assert_eq!(z.powf(1.0), z.to_polar());
654            for n in random_samples::<i32>() {
655                let x = n as FT * 0.01;
656                assert_eq!(z.powf(x).abs, z.abs().powf(x));
657                assert_eq!(z.powf(x).arg, z.arg() * x);
658            }
659        }
660        for n in random_samples::<i32>() {
661            let x = n as FT * 0.01;
662            assert_eq!(Rectangular::ZERO.powf(x.abs()), Polar::ZERO);
663            assert_eq!(Rectangular::ONE.powf(x), Polar::ONE);
664        }
665    }
666}
abels_complex/complex/f32/rectangular.rs

abels_complex/complex/f32/
rectangular.rs
