abels_complex/complex/f64/
rectangular.rs

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

abels_complex/complex/f64/
rectangular.rs
