abels_complex/
polar.rs

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

abels_complex/
polar.rs
