abels_complex/complex/f64/
polar.rs

1use core::f64::consts::{FRAC_PI_2, LN_2, LN_10, PI, TAU};
2use core::fmt;
3use core::ops::*;
4use core::write;
5#[cfg(feature = "libm")]
6use num_traits::real::Real;
7
8type Rectangular = crate::complex::rectangular::Complex<FT>;
9type Polar = crate::complex::polar::ComplexPolar<FT>;
10type FT = f64;
11
12impl Polar {
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 I: Self = Self::new(1.0, FRAC_PI_2);
16    pub const NEG_ONE: Self = Self::new(1.0, PI);
17    pub const NEG_I: Self = Self::new(1.0, -FRAC_PI_2);
18
19    /// Computes the conjugate.
20    pub const fn conjugate(self) -> Self {
21        Self::new(self.abs, -self.arg)
22    }
23
24    /// Computes the real component.
25    pub fn re(self) -> FT {
26        self.abs * self.arg.cos()
27    }
28
29    /// Computes the imaginary component.
30    pub fn im(self) -> FT {
31        self.abs * self.arg.sin()
32    }
33
34    /// Computes the squared absolute value.
35    pub const fn abs_sq(self) -> FT {
36        self.abs * self.abs
37    }
38
39    /// Computes the reciprocal.
40    pub fn recip(self) -> Self {
41        Self::new(self.abs.recip(), -self.arg)
42    }
43
44    /// Computes the principle square root.
45    pub fn sqrt(self) -> Self {
46        Self::new(self.abs.sqrt(), self.arg / 2.0)
47    }
48
49    /// Convert to rectangular form.
50    pub fn to_rectangular(self) -> Rectangular {
51        let (sin, cos) = self.arg.sin_cos();
52        self.abs * Rectangular::new(cos, sin)
53    }
54
55    /// Computes `e^self` where `e` is the base of the natural logarithm.
56    pub fn exp(self) -> Self {
57        self.to_rectangular().exp()
58    }
59
60    /// Computes the principle natural logarithm.
61    pub fn ln(self) -> Rectangular {
62        Rectangular::new(self.abs.ln(), self.arg)
63    }
64
65    /// Computes the principle logarithm in base 2.
66    pub fn log2(self) -> Rectangular {
67        self.ln() / LN_2
68    }
69
70    /// Computes the principle logarithm in base 10.
71    pub fn log10(self) -> Rectangular {
72        self.ln() / LN_10
73    }
74
75    /// Raises `self` to a floating point power.
76    pub fn powf(self, x: FT) -> Self {
77        if x < 0.0 && self.abs == 0.0 {
78            return Self::ZERO;
79        }
80        Self::new(self.abs.powf(x), self.arg * x)
81    }
82
83    /// Raises `self` to an integer power.
84    pub fn powi(self, n: i32) -> Self {
85        if n < 0 && self.abs == 0.0 {
86            return Self::ZERO;
87        }
88        Self::new(self.abs.powi(n), self.arg * n as FT)
89    }
90
91    /// Normalizes the absolute value and the argument into the range `[0, ∞)` and `(-π, +π]` respectively.
92    pub fn normalize(mut self) -> Self {
93        #[cfg(feature = "libm")]
94        {
95            self.arg = num_traits::Euclid::rem_euclid(&self.arg, &TAU);
96        }
97        #[cfg(not(feature = "libm"))]
98        {
99            self.arg = self.arg.rem_euclid(TAU);
100        }
101        if self.abs < 0.0 {
102            self.abs = -self.abs;
103            if self.arg <= 0.0 {
104                self.arg += PI;
105            } else {
106                self.arg -= PI;
107            }
108        } else {
109            if self.arg > PI {
110                self.arg -= TAU;
111            } else if self.arg <= -PI {
112                self.arg += TAU;
113            }
114        }
115        self
116    }
117}
118
119impl Mul for Polar {
120    type Output = Self;
121    fn mul(mut self, other: Self) -> Self {
122        self *= other;
123        self
124    }
125}
126
127impl Mul<FT> for Polar {
128    type Output = Self;
129    fn mul(mut self, re: FT) -> Self::Output {
130        self *= re;
131        self
132    }
133}
134
135impl Mul<Polar> for FT {
136    type Output = Polar;
137    fn mul(self, mut other: Self::Output) -> Self::Output {
138        other *= self;
139        other
140    }
141}
142
143impl MulAssign for Polar {
144    fn mul_assign(&mut self, other: Self) {
145        self.abs *= other.abs;
146        self.arg += other.arg;
147    }
148}
149
150impl MulAssign<FT> for Polar {
151    fn mul_assign(&mut self, re: FT) {
152        self.abs *= re;
153    }
154}
155
156impl Div for Polar {
157    type Output = Self;
158    fn div(mut self, other: Self) -> Self {
159        self /= other;
160        self
161    }
162}
163
164impl Div<FT> for Polar {
165    type Output = Self;
166    fn div(mut self, re: FT) -> Self {
167        self /= re;
168        self
169    }
170}
171
172impl Div<Polar> for FT {
173    type Output = Polar;
174    fn div(self, other: Self::Output) -> Self::Output {
175        self * other.recip()
176    }
177}
178
179impl DivAssign for Polar {
180    fn div_assign(&mut self, other: Self) {
181        *self *= other.recip();
182    }
183}
184
185impl DivAssign<FT> for Polar {
186    fn div_assign(&mut self, re: FT) {
187        self.abs /= re;
188    }
189}
190
191impl Neg for Polar {
192    type Output = Self;
193    fn neg(mut self) -> Self {
194        self.abs = -self.abs;
195        self
196    }
197}
198
199impl fmt::Display for Polar {
200    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
201        fn fmt_x(f: &mut fmt::Formatter, x: FT) -> fmt::Result {
202            if let Some(p) = f.precision() {
203                write!(f, "{x:.*}", p)
204            } else {
205                write!(f, "{x}")
206            }
207        }
208        let pi_radians = self.arg / PI;
209        fmt_x(f, self.abs)?;
210        if pi_radians == 0.0 || self.abs == 0.0 {
211            Ok(())
212        } else if pi_radians == 1.0 {
213            write!(f, "e^iπ")
214        } else {
215            write!(f, "e^")?;
216            fmt_x(f, pi_radians)?;
217            write!(f, "iπ")
218        }
219    }
220}
221
222#[cfg(feature = "rand")]
223impl rand::distr::Distribution<Polar> for rand::distr::StandardUniform {
224    fn sample<R: rand::Rng + ?Sized>(&self, rng: &mut R) -> Polar {
225        Polar::new(self.sample(rng), rng.random_range((-PI).next_up()..=PI))
226    }
227}
228
229#[cfg(feature = "approx")]
230use approx::{AbsDiffEq, RelativeEq, UlpsEq};
231
232#[cfg(feature = "approx")]
233impl AbsDiffEq for Polar {
234    type Epsilon = <FT as AbsDiffEq>::Epsilon;
235    fn default_epsilon() -> Self::Epsilon {
236        FT::default_epsilon()
237    }
238    fn abs_diff_eq(&self, other: &Self, epsilon: Self::Epsilon) -> bool {
239        FT::abs_diff_eq(&self.abs, &other.abs, epsilon)
240            && FT::abs_diff_eq(&self.arg, &other.arg, epsilon)
241    }
242}
243
244#[cfg(feature = "approx")]
245impl RelativeEq for Polar {
246    fn default_max_relative() -> Self::Epsilon {
247        FT::default_max_relative()
248    }
249    fn relative_eq(
250        &self,
251        other: &Self,
252        epsilon: Self::Epsilon,
253        max_relative: Self::Epsilon,
254    ) -> bool {
255        FT::relative_eq(&self.abs, &other.abs, epsilon, max_relative)
256            && FT::relative_eq(&self.arg, &other.arg, epsilon, max_relative)
257    }
258}
259
260#[cfg(feature = "approx")]
261impl UlpsEq for Polar {
262    fn default_max_ulps() -> u32 {
263        FT::default_max_ulps()
264    }
265    fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool {
266        FT::ulps_eq(&self.abs, &other.abs, epsilon, max_ulps)
267            && FT::ulps_eq(&self.arg, &other.arg, epsilon, max_ulps)
268    }
269}
270
271impl From<FT> for Polar {
272    fn from(value: FT) -> Self {
273        Self::new(value, 0.0)
274    }
275}
276
277#[cfg(test)]
278mod tests {
279    use super::*;
280    use approx::*;
281    use core::f64::consts::{E, FRAC_PI_2, PI};
282    use rand::{
283        Rng, SeedableRng,
284        distr::{Distribution, StandardUniform, Uniform, uniform::*},
285        rngs::StdRng,
286    };
287
288    const NUM_SAMPLES: usize = 100;
289
290    fn random_samples<T>() -> impl core::iter::Iterator<Item = T>
291    where
292        StandardUniform: Distribution<T>,
293    {
294        StdRng::seed_from_u64(21)
295            .sample_iter(StandardUniform)
296            .take(NUM_SAMPLES)
297    }
298
299    fn uniform_samples<T>(low: T, high: T) -> impl core::iter::Iterator<Item = T>
300    where
301        T: SampleUniform,
302    {
303        StdRng::seed_from_u64(21)
304            .sample_iter(Uniform::new(low, high).unwrap())
305            .take(NUM_SAMPLES)
306    }
307
308    #[test]
309    fn multiplication() {
310        for z0 in random_samples::<Polar>() {
311            for z1 in random_samples::<Polar>() {
312                let z = z0 * z1;
313                assert_eq!(z.abs, z0.abs * z1.abs);
314                assert_eq!(z.arg, z0.arg + z1.arg);
315
316                let z = z0 * z1.re();
317                assert_eq!(z.normalize().abs, z0.abs * z1.re().abs());
318                assert_eq!(z.arg, z0.arg);
319
320                let z = z0.re() * z1;
321                assert_eq!(z.normalize().abs, z0.re().abs() * z1.abs);
322                assert_eq!(z.arg, z1.arg);
323
324                let mut z = z0;
325                z *= z1;
326                assert_eq!(z, z0 * z1);
327
328                let mut z = z0;
329                z *= z1.re();
330                assert_eq!(z, z0 * z1.re());
331            }
332            assert_eq!(z0 * Polar::ONE, z0);
333            assert_eq!((z0 * Polar::ZERO).abs, 0.0);
334            assert_eq!((z0 * 0.0).abs, 0.0);
335        }
336    }
337
338    #[test]
339    fn division() {
340        for z0 in random_samples::<Polar>() {
341            for z1 in random_samples::<Polar>() {
342                let z = z0 / z1;
343                assert_ulps_eq!(z.abs, z0.abs / z1.abs);
344                assert_ulps_eq!(z.arg, z0.arg - z1.arg);
345
346                let z = z0 / z1.re();
347                assert_ulps_eq!(z.normalize().abs, z0.abs / z1.re().abs());
348                assert_ulps_eq!(z.arg, z0.arg);
349
350                let z = z0.re() / z1;
351                assert_ulps_eq!(z.normalize().abs, z0.re().abs() / z1.abs);
352                assert_ulps_eq!(z.arg, -z1.arg);
353
354                let mut z = z0;
355                z /= z1;
356                assert_eq!(z, z0 / z1);
357
358                let mut z = z0;
359                z /= z1.re();
360                assert_eq!(z, z0 / z1.re());
361            }
362            assert_ulps_eq!(z0 / z0, Polar::ONE);
363            assert_eq!((Polar::ZERO / z0).abs, 0.0);
364        }
365    }
366
367    #[test]
368    fn negation() {
369        assert_eq!((-Polar::ONE).normalize(), Polar::NEG_ONE);
370        assert_eq!((-Polar::I).normalize(), Polar::NEG_I);
371        assert_eq!((-Polar::NEG_ONE).normalize(), Polar::ONE);
372        assert_ulps_eq!((-Polar::NEG_I).normalize(), Polar::I);
373    }
374
375    #[test]
376    fn reciprocal() {
377        for z in random_samples::<Polar>() {
378            assert_ulps_eq!(z.recip(), 1.0 / z);
379            assert_ulps_eq!(z * z.recip(), Polar::ONE);
380        }
381        assert_eq!(Polar::ONE.recip(), Polar::ONE);
382        assert_eq!(Polar::I.recip(), Polar::NEG_I);
383        assert_eq!(Polar::NEG_ONE.recip().normalize(), Polar::NEG_ONE);
384        assert_eq!(Polar::NEG_I.recip(), Polar::I);
385    }
386
387    #[test]
388    fn sqrt() {
389        for z in random_samples::<Polar>() {
390            assert_eq!(z.sqrt().abs, z.abs.sqrt());
391            assert_eq!(z.sqrt().arg, z.arg / 2.0);
392        }
393        assert_eq!(Polar::ONE.sqrt(), Polar::ONE);
394        assert_eq!(Polar::NEG_ONE.sqrt(), Polar::I);
395        assert_eq!(Polar::ONE.sqrt(), Polar::ONE);
396    }
397
398    #[test]
399    fn abs() {
400        for z in random_samples::<Polar>() {
401            assert_eq!(z.abs_sq(), z.abs * z.abs);
402        }
403        assert_eq!(Polar::ONE.abs, 1.0);
404        assert_eq!(Polar::I.abs, 1.0);
405        assert_eq!(Polar::NEG_ONE.abs, 1.0);
406        assert_eq!(Polar::NEG_I.abs, 1.0);
407    }
408
409    #[test]
410    fn conjugate() {
411        for z in random_samples::<Polar>() {
412            assert_eq!(z.conjugate().re(), z.re());
413            assert_eq!(z.conjugate().im(), -z.im());
414            assert_eq!(z.conjugate().conjugate(), z);
415        }
416        assert_ulps_eq!(Polar::ONE.conjugate().normalize(), Polar::ONE);
417        assert_ulps_eq!(Polar::I.conjugate().normalize(), Polar::NEG_I);
418        assert_ulps_eq!(Polar::NEG_ONE.conjugate().normalize(), Polar::NEG_ONE);
419        assert_ulps_eq!(Polar::NEG_I.conjugate().normalize(), Polar::I);
420    }
421
422    #[test]
423    fn arg() {
424        assert_eq!(Polar::ONE.arg, 0.0);
425        assert_eq!(Polar::I.arg, FRAC_PI_2);
426        assert_eq!(Polar::NEG_ONE.arg, PI);
427        assert_eq!(Polar::NEG_I.arg, -FRAC_PI_2);
428    }
429
430    #[test]
431    fn exp() {
432        for z in random_samples::<Polar>() {
433            assert_eq!(z.exp().abs, z.re().exp());
434            assert_eq!(z.exp().arg, z.im());
435            assert_ulps_eq!(z.exp().ln(), z.to_rectangular());
436        }
437        assert_ulps_eq!(Polar::ONE.exp(), Polar::new(E, 0.0));
438        assert_ulps_eq!(Polar::I.exp(), Polar::new(1.0, 1.0));
439        assert_ulps_eq!(Polar::NEG_ONE.exp(), Polar::new(E.recip(), 0.0));
440        assert_ulps_eq!(Polar::NEG_I.exp(), Polar::new(1.0, -1.0));
441    }
442
443    #[test]
444    fn log() {
445        for z in random_samples::<Polar>() {
446            assert_eq!(z.ln().re, z.abs.ln());
447            assert_eq!(z.ln().im, z.arg);
448            assert_ulps_eq!(z.ln().exp(), z);
449        }
450        assert_eq!(Polar::ONE.ln(), Rectangular::ZERO);
451        assert_eq!(Polar::I.ln(), Rectangular::I * FRAC_PI_2);
452        assert_eq!(Polar::NEG_ONE.ln(), Rectangular::I * PI);
453        assert_eq!(Polar::NEG_I.ln(), Rectangular::I * -FRAC_PI_2);
454
455        assert_ulps_eq!(Polar::new(E, 0.0).ln(), Rectangular::ONE);
456        assert_ulps_eq!(Polar::new(2.0, 0.0).log2(), Rectangular::ONE);
457        assert_ulps_eq!(Polar::new(10.0, 0.0).log10(), Rectangular::ONE);
458    }
459
460    #[test]
461    fn powi() {
462        for z in random_samples::<Polar>() {
463            assert_eq!(z.powi(0), Polar::ONE);
464            assert_eq!(z.powi(1), z);
465            for n in random_samples::<i32>() {
466                assert_eq!(z.powi(n).abs, z.abs.powi(n));
467                assert_eq!(z.powi(n).arg, z.arg * n as FT);
468            }
469        }
470        for n in random_samples::<i32>() {
471            assert_eq!(Polar::ZERO.powi(n.abs()), Polar::ZERO);
472            assert_eq!(Polar::ONE.powi(n), Polar::ONE);
473        }
474    }
475
476    #[test]
477    fn powf() {
478        for z in random_samples::<Polar>() {
479            assert_eq!(z.powf(0.0), Polar::ONE);
480            assert_eq!(z.powf(1.0), z);
481            for n in random_samples::<i32>() {
482                let x = n as FT * 0.01;
483                assert_eq!(z.powf(x).abs, z.abs.powf(x));
484                assert_eq!(z.powf(x).arg, z.arg * x);
485            }
486        }
487        for n in random_samples::<i32>() {
488            let x = n as FT * 0.01;
489            assert_eq!(Polar::ZERO.powf(x.abs()), Polar::ZERO);
490            assert_eq!(Polar::ONE.powf(x), Polar::ONE);
491        }
492    }
493
494    #[test]
495    fn normalize() {
496        for z in random_samples::<Polar>() {
497            for n in uniform_samples::<i32>(-99, 99) {
498                let w = Polar::new(z.abs, z.arg + n as FT * TAU);
499                assert_ulps_eq!(z, w.normalize(), epsilon = 2000.0 * FT::EPSILON);
500
501                assert_ulps_eq!(
502                    Polar::new(-z.abs, z.arg).normalize(),
503                    Polar::new(z.abs, z.arg + PI).normalize(),
504                    epsilon = 2000.0 * FT::EPSILON
505                );
506            }
507        }
508    }
509}
abels_complex/complex/f64/polar.rs

abels_complex/complex/f64/
polar.rs
