abels_complex/complex/f64/
rectangular.rs

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

abels_complex/complex/f64/
rectangular.rs
