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        use rand::RngExt;
32        Polar::new(self.sample(rng), rng.random_range((-PI).next_up()..=PI))
33    }
34}
35
36#[cfg(test)]
37mod tests {
38    use super::*;
39    use approx::*;
40    use core::f64::consts::{E, FRAC_PI_4, LN_2, TAU};
41    use rand::{
42        RngExt, SeedableRng,
43        distr::{Distribution, StandardUniform, Uniform, uniform::*},
44        rngs::StdRng,
45    };
46
47    type Rectangular = crate::complex::rectangular::Complex<FT>;
48
49    const NUM_SAMPLES: usize = 100;
50    const SEED: u64 = 21;
51
52    fn random_samples<T>() -> impl core::iter::Iterator<Item = T>
53    where
54        StandardUniform: Distribution<T>,
55    {
56        StdRng::seed_from_u64(SEED)
57            .sample_iter(StandardUniform)
58            .take(NUM_SAMPLES)
59    }
60
61    fn uniform_samples<T>(low: T, high: T) -> impl core::iter::Iterator<Item = T>
62    where
63        T: SampleUniform,
64    {
65        StdRng::seed_from_u64(SEED)
66            .sample_iter(Uniform::new(low, high).unwrap())
67            .take(NUM_SAMPLES)
68    }
69
70    #[test]
71    fn multiplication() {
72        for z0 in random_samples::<Polar>() {
73            for z1 in random_samples::<Polar>() {
74                let z = z0 * z1;
75                assert_eq!(z.abs, z0.abs * z1.abs);
76                assert_eq!(z.arg, z0.arg + z1.arg);
77
78                let z = z0 * z1.re();
79                assert_eq!(z.normalize().abs, z0.abs * z1.re().abs());
80                assert_eq!(z.arg, z0.arg);
81
82                let z = z0.re() * z1;
83                assert_eq!(z.normalize().abs, z0.re().abs() * z1.abs);
84                assert_eq!(z.arg, z1.arg);
85
86                let mut z = z0;
87                z *= z1;
88                assert_eq!(z, z0 * z1);
89
90                let mut z = z0;
91                z *= z1.re();
92                assert_eq!(z, z0 * z1.re());
93            }
94            assert_eq!(z0 * Polar::ONE, z0);
95            assert_eq!((z0 * Polar::ZERO).abs, 0.0);
96            assert_eq!((z0 * 0.0).abs, 0.0);
97        }
98    }
99
100    #[test]
101    fn division() {
102        for z0 in random_samples::<Polar>() {
103            for z1 in random_samples::<Polar>() {
104                let z = z0 / z1;
105                assert_ulps_eq!(z.abs, z0.abs / z1.abs);
106                assert_ulps_eq!(z.arg, z0.arg - z1.arg);
107
108                let z = z0 / z1.re();
109                assert_ulps_eq!(z.normalize().abs, z0.abs / z1.re().abs());
110                assert_ulps_eq!(z.arg, z0.arg);
111
112                let z = z0.re() / z1;
113                assert_ulps_eq!(z.normalize().abs, z0.re().abs() / z1.abs);
114                assert_ulps_eq!(z.arg, -z1.arg);
115
116                let mut z = z0;
117                z /= z1;
118                assert_eq!(z, z0 / z1);
119
120                let mut z = z0;
121                z /= z1.re();
122                assert_eq!(z, z0 / z1.re());
123            }
124            assert_ulps_eq!(z0 / z0, Polar::ONE);
125            assert_eq!((Polar::ZERO / z0).abs, 0.0);
126        }
127    }
128
129    #[test]
130    fn negation() {
131        assert_eq!((-Polar::ONE).normalize(), Polar::NEG_ONE);
132        assert_eq!((-Polar::I).normalize(), Polar::NEG_I);
133        assert_eq!((-Polar::NEG_ONE).normalize(), Polar::ONE);
134        assert_ulps_eq!((-Polar::NEG_I).normalize(), Polar::I);
135    }
136
137    #[test]
138    fn reciprocal() {
139        for z in random_samples::<Polar>() {
140            assert_ulps_eq!(z.recip(), 1.0 / z);
141            assert_ulps_eq!(z * z.recip(), Polar::ONE);
142        }
143        assert_eq!(Polar::ONE.recip(), Polar::ONE);
144        assert_eq!(Polar::I.recip(), Polar::NEG_I);
145        assert_eq!(Polar::NEG_ONE.recip().normalize(), Polar::NEG_ONE);
146        assert_eq!(Polar::NEG_I.recip(), Polar::I);
147    }
148
149    #[test]
150    fn sqrt() {
151        for z in random_samples::<Polar>() {
152            assert_eq!(z.sqrt().abs, z.abs.sqrt());
153            assert_eq!(z.sqrt().arg, z.arg / 2.0);
154        }
155        assert_eq!(Polar::ONE.sqrt(), Polar::ONE);
156        assert_eq!(Polar::NEG_ONE.sqrt(), Polar::I);
157        assert_eq!(Polar::ONE.sqrt(), Polar::ONE);
158    }
159
160    #[test]
161    fn abs() {
162        for z in random_samples::<Polar>() {
163            assert_eq!(z.abs_sq(), z.abs * z.abs);
164        }
165        assert_eq!(Polar::ONE.abs, 1.0);
166        assert_eq!(Polar::I.abs, 1.0);
167        assert_eq!(Polar::NEG_ONE.abs, 1.0);
168        assert_eq!(Polar::NEG_I.abs, 1.0);
169    }
170
171    #[test]
172    fn conjugate() {
173        for z in random_samples::<Polar>() {
174            assert_eq!(z.conjugate().re(), z.re());
175            assert_eq!(z.conjugate().im(), -z.im());
176            assert_eq!(z.conjugate().conjugate(), z);
177        }
178        assert_ulps_eq!(Polar::ONE.conjugate().normalize(), Polar::ONE);
179        assert_ulps_eq!(Polar::I.conjugate().normalize(), Polar::NEG_I);
180        assert_ulps_eq!(Polar::NEG_ONE.conjugate().normalize(), Polar::NEG_ONE);
181        assert_ulps_eq!(Polar::NEG_I.conjugate().normalize(), Polar::I);
182    }
183
184    #[test]
185    fn arg() {
186        assert_eq!(Polar::ONE.arg, 0.0);
187        assert_eq!(Polar::I.arg, FRAC_PI_2);
188        assert_eq!(Polar::NEG_ONE.arg, PI);
189        assert_eq!(Polar::NEG_I.arg, -FRAC_PI_2);
190    }
191
192    #[test]
193    fn exp() {
194        for z in random_samples::<Polar>() {
195            assert_eq!(z.exp().abs, z.re().exp());
196            assert_eq!(z.exp().arg, z.im());
197            assert_ulps_eq!(z.exp().ln(), z.to_rectangular());
198        }
199        assert_ulps_eq!(Polar::ONE.exp(), Polar::new(E, 0.0));
200        assert_ulps_eq!(Polar::I.exp(), Polar::new(1.0, 1.0));
201        assert_ulps_eq!(Polar::NEG_ONE.exp(), Polar::new(E.recip(), 0.0));
202        assert_ulps_eq!(Polar::NEG_I.exp(), Polar::new(1.0, -1.0));
203    }
204
205    #[test]
206    fn expm1() {
207        for z in random_samples::<Polar>() {
208            assert_ulps_eq!(z.expm1(), z.to_rectangular().expm1(), max_ulps = 2);
209        }
210        assert_eq!(Polar::ZERO.expm1(), Rectangular::ZERO);
211        assert_ulps_eq!(Polar::ONE.expm1(), Rectangular::new(E - 1.0, 0.0));
212        assert_ulps_eq!(Polar::I.expm1(), Polar::I.exp().to_rectangular() - 1.0);
213    }
214
215    #[test]
216    fn log() {
217        for z in random_samples::<Polar>() {
218            assert_eq!(z.ln().re, z.abs.ln());
219            assert_eq!(z.ln().im, z.arg);
220            assert_ulps_eq!(z.ln().exp(), z);
221        }
222        assert_eq!(Polar::ONE.ln(), Rectangular::ZERO);
223        assert_eq!(Polar::I.ln(), Rectangular::I * FRAC_PI_2);
224        assert_eq!(Polar::NEG_ONE.ln(), Rectangular::I * PI);
225        assert_eq!(Polar::NEG_I.ln(), Rectangular::I * -FRAC_PI_2);
226
227        assert_ulps_eq!(Polar::new(E, 0.0).ln(), Rectangular::ONE);
228        assert_ulps_eq!(Polar::new(2.0, 0.0).log2(), Rectangular::ONE);
229        assert_ulps_eq!(Polar::new(10.0, 0.0).log10(), Rectangular::ONE);
230    }
231
232    #[test]
233    fn ln_1p() {
234        for z in random_samples::<Polar>() {
235            assert_ulps_eq!(
236                z.ln_1p(),
237                z.to_rectangular().ln_1p(),
238                epsilon = 2.0 * FT::EPSILON
239            );
240        }
241        assert_eq!(Polar::ZERO.ln_1p(), Rectangular::ZERO);
242        assert_ulps_eq!(Polar::ONE.ln_1p(), Rectangular::new(LN_2, 0.0));
243        assert_ulps_eq!(Polar::I.ln_1p(), Rectangular::new(LN_2 / 2.0, FRAC_PI_4));
244    }
245
246    #[test]
247    fn powi() {
248        for z in random_samples::<Polar>() {
249            assert_eq!(z.powi(0), Polar::ONE);
250            assert_eq!(z.powi(1), z);
251            for n in random_samples::<i32>() {
252                assert_eq!(z.powi(n).abs, z.abs.powi(n));
253                assert_eq!(z.powi(n).arg, z.arg * n as FT);
254            }
255        }
256        for n in random_samples::<i32>() {
257            assert_eq!(Polar::ZERO.powi(n.abs()), Polar::ZERO);
258            assert_eq!(Polar::ONE.powi(n), Polar::ONE);
259        }
260    }
261
262    #[test]
263    fn powf() {
264        for z in random_samples::<Polar>() {
265            assert_eq!(z.powf(0.0), Polar::ONE);
266            assert_eq!(z.powf(1.0), z);
267            for n in random_samples::<i32>() {
268                let x = n as FT * 0.01;
269                assert_eq!(z.powf(x).abs, z.abs.powf(x));
270                assert_eq!(z.powf(x).arg, z.arg * x);
271            }
272        }
273        for n in random_samples::<i32>() {
274            let x = n as FT * 0.01;
275            assert_eq!(Polar::ZERO.powf(x.abs()), Polar::ZERO);
276            assert_eq!(Polar::ONE.powf(x), Polar::ONE);
277        }
278    }
279
280    #[test]
281    fn ln_branch() {
282        // k=0 matches principal value
283        for z in random_samples::<Polar>() {
284            assert_eq!(z.ln_branch(0), z.ln());
285        }
286        // ln_k(z) = ln(z) + 2πki
287        for z in random_samples::<Polar>() {
288            for k in -3_i32..=3 {
289                assert_ulps_eq!(z.ln_branch(k), z.ln() + Rectangular::I * (k as FT * TAU));
290            }
291        }
292        // Known values: ln_1(1) = 2πi
293        assert_ulps_eq!(Polar::ONE.ln_branch(1), Rectangular::new(0.0, TAU));
294        // ln_(-1)(-1) = ln|-1| + i(π - 2π) = -πi
295        assert_ulps_eq!(Polar::NEG_ONE.ln_branch(-1), Rectangular::new(0.0, -PI));
296        // log2_branch and log10_branch consistency
297        for z in random_samples::<Polar>() {
298            assert_ulps_eq!(z.log2_branch(1), z.ln_branch(1) / LN_2);
299            assert_ulps_eq!(z.log10_branch(1), z.ln_branch(1) / core::f64::consts::LN_10);
300        }
301    }
302
303    #[test]
304    fn ln_1p_branch() {
305        // k=0 matches principal value
306        for z in random_samples::<Polar>() {
307            assert_eq!(z.ln_1p_branch(0), z.ln_1p());
308        }
309        // ln_1p_branch(k) = ln_1p + 2πki
310        for z in random_samples::<Polar>() {
311            for k in -3_i32..=3 {
312                assert_ulps_eq!(
313                    z.ln_1p_branch(k),
314                    z.ln_1p() + Rectangular::I * (k as FT * TAU)
315                );
316            }
317        }
318    }
319
320    #[test]
321    fn sqrt_branch() {
322        // k=0 matches principal value
323        for z in random_samples::<Polar>() {
324            assert_eq!(z.sqrt_branch(0), z.sqrt());
325        }
326        // k=1 gives the other square root: arg offset by π
327        for z in random_samples::<Polar>() {
328            assert_ulps_eq!(z.sqrt_branch(1).abs, z.sqrt().abs);
329            assert_ulps_eq!(z.sqrt_branch(1).arg, z.sqrt().arg + PI);
330        }
331        // The two roots multiply to give z (after normalizing abs)
332        assert_ulps_eq!(Polar::ONE.sqrt_branch(1), Polar::new(1.0, PI));
333    }
334
335    #[test]
336    fn nth_root() {
337        // k=0 matches powf(1/n) in terms of abs; arg = arg/n
338        for z in random_samples::<Polar>() {
339            let r = z.nth_root(3, 0);
340            assert_ulps_eq!(r.abs, z.abs.powf(1.0 / 3.0));
341            assert_ulps_eq!(r.arg, z.arg / 3.0);
342        }
343        // Cube roots of unity: nth_root(1+0i, 3, k) for k=0,1,2
344        let one = Polar::ONE;
345        let r0 = one.nth_root(3, 0);
346        let r1 = one.nth_root(3, 1);
347        let r2 = one.nth_root(3, 2);
348        assert_ulps_eq!(r0, Polar::ONE);
349        assert_ulps_eq!(r1.abs, 1.0);
350        assert_ulps_eq!(r1.arg, TAU / 3.0);
351        assert_ulps_eq!(r2.abs, 1.0);
352        assert_ulps_eq!(r2.arg, TAU * 2.0 / 3.0);
353        // Each root cubed gives the original
354        for z in random_samples::<Polar>() {
355            for k in 0..3_i32 {
356                assert_ulps_eq!(z.nth_root(3, k).powi(3).abs, z.abs, epsilon = 1e-13);
357            }
358        }
359    }
360
361    #[test]
362    fn pow_rational() {
363        // p/q = 1/1: same as powi(1)
364        for z in random_samples::<Polar>() {
365            assert_ulps_eq!(z.pow_rational(1, 1, 0), z.powi(1));
366        }
367        // p/q = 1/2, k=0: same as sqrt
368        for z in random_samples::<Polar>() {
369            assert_eq!(z.pow_rational(1, 2, 0), z.sqrt());
370        }
371        // p/q = 1/3: cube root, k=0..2 give distinct values
372        for z in random_samples::<Polar>() {
373            for k in 0..3_i32 {
374                let r = z.pow_rational(1, 3, k);
375                assert_ulps_eq!(r.powi(3).abs, z.abs, epsilon = 1e-13);
376                assert_ulps_eq!(r.abs, z.abs.powf(1.0 / 3.0), epsilon = 1e-13);
377            }
378        }
379        // p/q = 2/3: two-thirds power — 3 distinct values
380        let z = Polar::new(8.0, 0.0); // 8^(2/3) = 4
381        assert_ulps_eq!(z.pow_rational(2, 3, 0).abs, 4.0);
382        // 2/4 has only 2 distinct values (same as 1/2 after normalization)
383        for z in random_samples::<Polar>() {
384            assert_relative_eq!(
385                z.pow_rational(2, 4, 0).normalize(),
386                z.pow_rational(1, 2, 0).normalize(),
387                max_relative = 1e-14
388            );
389            assert_relative_eq!(
390                z.pow_rational(2, 4, 2).normalize(),
391                z.pow_rational(1, 2, 0).normalize(),
392                max_relative = 1e-14
393            );
394        }
395    }
396
397    #[test]
398    fn normalize() {
399        for z in random_samples::<Polar>() {
400            for n in uniform_samples::<i32>(-99, 99) {
401                let w = Polar::new(z.abs, z.arg + n as FT * TAU);
402                assert_ulps_eq!(z, w.normalize(), epsilon = 2000.0 * FT::EPSILON);
403
404                assert_ulps_eq!(
405                    Polar::new(-z.abs, z.arg).normalize(),
406                    Polar::new(z.abs, z.arg + PI).normalize(),
407                    epsilon = 2000.0 * FT::EPSILON
408                );
409            }
410        }
411    }
412}
abels_complex/complex/f64/polar.rs

abels_complex/complex/f64/
polar.rs
