abels_complex/
polar.rs

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

abels_complex/
polar.rs
