abels_complex/complex/f64/
polar.rs

1use core::f64::consts::{FRAC_PI_2, PI};
2use core::ops::*;
3
4type Polar = crate::complex::polar::ComplexPolar<FT>;
5type FT = f64;
6
7impl Mul<Polar> for FT {
8    type Output = Polar;
9    fn mul(self, mut other: Self::Output) -> Self::Output {
10        other *= self;
11        other
12    }
13}
14
15impl Div<Polar> for FT {
16    type Output = Polar;
17    fn div(self, other: Self::Output) -> Self::Output {
18        self * other.recip()
19    }
20}
21
22impl Polar {
23    pub const I: Self = Self::new(1.0, FRAC_PI_2);
24    pub const NEG_ONE: Self = Self::new(1.0, PI);
25    pub const NEG_I: Self = Self::new(1.0, -FRAC_PI_2);
26}
27
28#[cfg(feature = "rand")]
29impl rand::distr::Distribution<Polar> for rand::distr::StandardUniform {
30    fn sample<R: rand::Rng + ?Sized>(&self, rng: &mut R) -> Polar {
31        Polar::new(self.sample(rng), rng.random_range((-PI).next_up()..=PI))
32    }
33}
34
35#[cfg(test)]
36mod tests {
37    use super::*;
38    use approx::*;
39    use core::f64::consts::{E, TAU};
40    use rand::{
41        Rng, SeedableRng,
42        distr::{Distribution, StandardUniform, Uniform, uniform::*},
43        rngs::StdRng,
44    };
45
46    type Rectangular = crate::complex::rectangular::Complex<FT>;
47
48    const NUM_SAMPLES: usize = 100;
49    const SEED: u64 = 21;
50
51    fn random_samples<T>() -> impl core::iter::Iterator<Item = T>
52    where
53        StandardUniform: Distribution<T>,
54    {
55        StdRng::seed_from_u64(SEED)
56            .sample_iter(StandardUniform)
57            .take(NUM_SAMPLES)
58    }
59
60    fn uniform_samples<T>(low: T, high: T) -> impl core::iter::Iterator<Item = T>
61    where
62        T: SampleUniform,
63    {
64        StdRng::seed_from_u64(SEED)
65            .sample_iter(Uniform::new(low, high).unwrap())
66            .take(NUM_SAMPLES)
67    }
68
69    #[test]
70    fn multiplication() {
71        for z0 in random_samples::<Polar>() {
72            for z1 in random_samples::<Polar>() {
73                let z = z0 * z1;
74                assert_eq!(z.abs, z0.abs * z1.abs);
75                assert_eq!(z.arg, z0.arg + z1.arg);
76
77                let z = z0 * z1.re();
78                assert_eq!(z.normalize().abs, z0.abs * z1.re().abs());
79                assert_eq!(z.arg, z0.arg);
80
81                let z = z0.re() * z1;
82                assert_eq!(z.normalize().abs, z0.re().abs() * z1.abs);
83                assert_eq!(z.arg, z1.arg);
84
85                let mut z = z0;
86                z *= z1;
87                assert_eq!(z, z0 * z1);
88
89                let mut z = z0;
90                z *= z1.re();
91                assert_eq!(z, z0 * z1.re());
92            }
93            assert_eq!(z0 * Polar::ONE, z0);
94            assert_eq!((z0 * Polar::ZERO).abs, 0.0);
95            assert_eq!((z0 * 0.0).abs, 0.0);
96        }
97    }
98
99    #[test]
100    fn division() {
101        for z0 in random_samples::<Polar>() {
102            for z1 in random_samples::<Polar>() {
103                let z = z0 / z1;
104                assert_ulps_eq!(z.abs, z0.abs / z1.abs);
105                assert_ulps_eq!(z.arg, z0.arg - z1.arg);
106
107                let z = z0 / z1.re();
108                assert_ulps_eq!(z.normalize().abs, z0.abs / z1.re().abs());
109                assert_ulps_eq!(z.arg, z0.arg);
110
111                let z = z0.re() / z1;
112                assert_ulps_eq!(z.normalize().abs, z0.re().abs() / z1.abs);
113                assert_ulps_eq!(z.arg, -z1.arg);
114
115                let mut z = z0;
116                z /= z1;
117                assert_eq!(z, z0 / z1);
118
119                let mut z = z0;
120                z /= z1.re();
121                assert_eq!(z, z0 / z1.re());
122            }
123            assert_ulps_eq!(z0 / z0, Polar::ONE);
124            assert_eq!((Polar::ZERO / z0).abs, 0.0);
125        }
126    }
127
128    #[test]
129    fn negation() {
130        assert_eq!((-Polar::ONE).normalize(), Polar::NEG_ONE);
131        assert_eq!((-Polar::I).normalize(), Polar::NEG_I);
132        assert_eq!((-Polar::NEG_ONE).normalize(), Polar::ONE);
133        assert_ulps_eq!((-Polar::NEG_I).normalize(), Polar::I);
134    }
135
136    #[test]
137    fn reciprocal() {
138        for z in random_samples::<Polar>() {
139            assert_ulps_eq!(z.recip(), 1.0 / z);
140            assert_ulps_eq!(z * z.recip(), Polar::ONE);
141        }
142        assert_eq!(Polar::ONE.recip(), Polar::ONE);
143        assert_eq!(Polar::I.recip(), Polar::NEG_I);
144        assert_eq!(Polar::NEG_ONE.recip().normalize(), Polar::NEG_ONE);
145        assert_eq!(Polar::NEG_I.recip(), Polar::I);
146    }
147
148    #[test]
149    fn sqrt() {
150        for z in random_samples::<Polar>() {
151            assert_eq!(z.sqrt().abs, z.abs.sqrt());
152            assert_eq!(z.sqrt().arg, z.arg / 2.0);
153        }
154        assert_eq!(Polar::ONE.sqrt(), Polar::ONE);
155        assert_eq!(Polar::NEG_ONE.sqrt(), Polar::I);
156        assert_eq!(Polar::ONE.sqrt(), Polar::ONE);
157    }
158
159    #[test]
160    fn abs() {
161        for z in random_samples::<Polar>() {
162            assert_eq!(z.abs_sq(), z.abs * z.abs);
163        }
164        assert_eq!(Polar::ONE.abs, 1.0);
165        assert_eq!(Polar::I.abs, 1.0);
166        assert_eq!(Polar::NEG_ONE.abs, 1.0);
167        assert_eq!(Polar::NEG_I.abs, 1.0);
168    }
169
170    #[test]
171    fn conjugate() {
172        for z in random_samples::<Polar>() {
173            assert_eq!(z.conjugate().re(), z.re());
174            assert_eq!(z.conjugate().im(), -z.im());
175            assert_eq!(z.conjugate().conjugate(), z);
176        }
177        assert_ulps_eq!(Polar::ONE.conjugate().normalize(), Polar::ONE);
178        assert_ulps_eq!(Polar::I.conjugate().normalize(), Polar::NEG_I);
179        assert_ulps_eq!(Polar::NEG_ONE.conjugate().normalize(), Polar::NEG_ONE);
180        assert_ulps_eq!(Polar::NEG_I.conjugate().normalize(), Polar::I);
181    }
182
183    #[test]
184    fn arg() {
185        assert_eq!(Polar::ONE.arg, 0.0);
186        assert_eq!(Polar::I.arg, FRAC_PI_2);
187        assert_eq!(Polar::NEG_ONE.arg, PI);
188        assert_eq!(Polar::NEG_I.arg, -FRAC_PI_2);
189    }
190
191    #[test]
192    fn exp() {
193        for z in random_samples::<Polar>() {
194            assert_eq!(z.exp().abs, z.re().exp());
195            assert_eq!(z.exp().arg, z.im());
196            assert_ulps_eq!(z.exp().ln(), z.to_rectangular());
197        }
198        assert_ulps_eq!(Polar::ONE.exp(), Polar::new(E, 0.0));
199        assert_ulps_eq!(Polar::I.exp(), Polar::new(1.0, 1.0));
200        assert_ulps_eq!(Polar::NEG_ONE.exp(), Polar::new(E.recip(), 0.0));
201        assert_ulps_eq!(Polar::NEG_I.exp(), Polar::new(1.0, -1.0));
202    }
203
204    #[test]
205    fn log() {
206        for z in random_samples::<Polar>() {
207            assert_eq!(z.ln().re, z.abs.ln());
208            assert_eq!(z.ln().im, z.arg);
209            assert_ulps_eq!(z.ln().exp(), z);
210        }
211        assert_eq!(Polar::ONE.ln(), Rectangular::ZERO);
212        assert_eq!(Polar::I.ln(), Rectangular::I * FRAC_PI_2);
213        assert_eq!(Polar::NEG_ONE.ln(), Rectangular::I * PI);
214        assert_eq!(Polar::NEG_I.ln(), Rectangular::I * -FRAC_PI_2);
215
216        assert_ulps_eq!(Polar::new(E, 0.0).ln(), Rectangular::ONE);
217        assert_ulps_eq!(Polar::new(2.0, 0.0).log2(), Rectangular::ONE);
218        assert_ulps_eq!(Polar::new(10.0, 0.0).log10(), Rectangular::ONE);
219    }
220
221    #[test]
222    fn powi() {
223        for z in random_samples::<Polar>() {
224            assert_eq!(z.powi(0), Polar::ONE);
225            assert_eq!(z.powi(1), z);
226            for n in random_samples::<i32>() {
227                assert_eq!(z.powi(n).abs, z.abs.powi(n));
228                assert_eq!(z.powi(n).arg, z.arg * n as FT);
229            }
230        }
231        for n in random_samples::<i32>() {
232            assert_eq!(Polar::ZERO.powi(n.abs()), Polar::ZERO);
233            assert_eq!(Polar::ONE.powi(n), Polar::ONE);
234        }
235    }
236
237    #[test]
238    fn powf() {
239        for z in random_samples::<Polar>() {
240            assert_eq!(z.powf(0.0), Polar::ONE);
241            assert_eq!(z.powf(1.0), z);
242            for n in random_samples::<i32>() {
243                let x = n as FT * 0.01;
244                assert_eq!(z.powf(x).abs, z.abs.powf(x));
245                assert_eq!(z.powf(x).arg, z.arg * x);
246            }
247        }
248        for n in random_samples::<i32>() {
249            let x = n as FT * 0.01;
250            assert_eq!(Polar::ZERO.powf(x.abs()), Polar::ZERO);
251            assert_eq!(Polar::ONE.powf(x), Polar::ONE);
252        }
253    }
254
255    #[test]
256    fn normalize() {
257        for z in random_samples::<Polar>() {
258            for n in uniform_samples::<i32>(-99, 99) {
259                let w = Polar::new(z.abs, z.arg + n as FT * TAU);
260                assert_ulps_eq!(z, w.normalize(), epsilon = 2000.0 * FT::EPSILON);
261
262                assert_ulps_eq!(
263                    Polar::new(-z.abs, z.arg).normalize(),
264                    Polar::new(z.abs, z.arg + PI).normalize(),
265                    epsilon = 2000.0 * FT::EPSILON
266                );
267            }
268        }
269    }
270}
abels_complex/complex/f64/polar.rs

abels_complex/complex/f64/
polar.rs
