abels_complex/complex/f32/
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 = f32;
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::f32::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 exp2() {
340        for z in random_samples::<Rectangular>() {
341            assert_eq!(z.exp2().abs, z.re.exp2());
342            assert_eq!(z.exp2().arg, z.im * LN_2);
343            assert_ulps_eq!(z.exp2().log2(), z);
344        }
345        assert_eq!(Rectangular::ONE.exp2(), Polar::new(2.0, 0.0));
346        assert_eq!(Rectangular::I.exp2(), Polar::new(1.0, LN_2));
347        assert_eq!(Rectangular::NEG_ONE.exp2(), Polar::new(0.5, 0.0));
348        assert_eq!(Rectangular::NEG_I.exp2(), Polar::new(1.0, -LN_2));
349    }
350
351    #[test]
352    fn expm1() {
353        for z in random_samples::<Rectangular>() {
354            assert_ulps_eq!(z.expm1(), z.exp().to_rectangular() - 1.0, max_ulps = 5);
355        }
356        assert_eq!(Rectangular::ZERO.expm1(), Rectangular::ZERO);
357        assert_ulps_eq!(Rectangular::ONE.expm1(), Rectangular::new(E - 1.0, 0.0));
358        assert_ulps_eq!(
359            Rectangular::I.expm1(),
360            Rectangular::I.exp().to_rectangular() - 1.0
361        );
362
363        let tiny = Rectangular::new(1e-8, 0.0);
364        assert_eq!(tiny.exp().to_rectangular() - 1.0, Rectangular::ZERO);
365        assert_ne!(tiny.expm1(), Rectangular::ZERO);
366    }
367
368    #[test]
369    fn log() {
370        for z in random_samples::<Rectangular>() {
371            assert_eq!(z.ln().re, z.abs().ln());
372            assert_eq!(z.ln().im, z.arg());
373            assert_ulps_eq!(z.ln().exp(), z.to_polar());
374        }
375        assert_eq!(Rectangular::ONE.ln(), Rectangular::ZERO);
376        assert_eq!(Rectangular::I.ln(), Rectangular::I * FRAC_PI_2);
377        assert_eq!(Rectangular::NEG_ONE.ln(), Rectangular::I * PI);
378        assert_eq!(Rectangular::NEG_I.ln(), Rectangular::I * -FRAC_PI_2);
379
380        assert_ulps_eq!(Rectangular::new(E, 0.0).ln(), Rectangular::ONE);
381        assert_ulps_eq!(Rectangular::new(2.0, 0.0).log2(), Rectangular::ONE);
382        assert_ulps_eq!(Rectangular::new(10.0, 0.0).log10(), Rectangular::ONE);
383    }
384
385    #[test]
386    fn ln_1p() {
387        for z in random_samples::<Rectangular>() {
388            // (z + 1).ln() is unreliable near the singularity at z = -1
389            let z = Rectangular::new(z.re.max(-0.75), z.im);
390            assert_ulps_eq!(z.ln_1p(), (z + 1.0).ln(), max_ulps = 5);
391        }
392        assert_eq!(Rectangular::ZERO.ln_1p(), Rectangular::ZERO);
393        assert_ulps_eq!(Rectangular::ONE.ln_1p(), Rectangular::new(LN_2, 0.0));
394        assert_ulps_eq!(
395            Rectangular::I.ln_1p(),
396            Rectangular::new(LN_2 / 2.0, FRAC_PI_4)
397        );
398        // Near zero: (z + 1).ln() loses all precision, ln_1p preserves it
399        let tiny = Rectangular::new(1e-8, 0.0);
400        assert_eq!((tiny + 1.0).ln(), Rectangular::ZERO); // 1.0 + 1e-8 == 1.0 in f32
401        assert_ne!(tiny.ln_1p(), Rectangular::ZERO);
402    }
403
404    #[test]
405    fn powi() {
406        for z in random_samples::<Rectangular>() {
407            assert_eq!(z.powi(0), Polar::ONE);
408            assert_eq!(z.powi(1), z.to_polar());
409            for n in random_samples::<i32>() {
410                assert_eq!(z.powi(n).abs, z.abs().powi(n));
411                assert_eq!(z.powi(n).arg, z.arg() * n as FT);
412            }
413        }
414        for n in random_samples::<i32>() {
415            assert_eq!(Rectangular::ZERO.powi(n.abs()), Polar::ZERO);
416            assert_eq!(Rectangular::ONE.powi(n), Polar::ONE);
417        }
418    }
419
420    #[test]
421    fn powf() {
422        for z in random_samples::<Rectangular>() {
423            assert_eq!(z.powf(0.0), Polar::ONE);
424            assert_eq!(z.powf(1.0), z.to_polar());
425            for n in random_samples::<i32>() {
426                let x = n as FT * 0.01;
427                assert_eq!(z.powf(x).abs, z.abs().powf(x));
428                assert_eq!(z.powf(x).arg, z.arg() * x);
429            }
430        }
431        for n in random_samples::<i32>() {
432            let x = n as FT * 0.01;
433            assert_eq!(Rectangular::ZERO.powf(x.abs()), Polar::ZERO);
434            assert_eq!(Rectangular::ONE.powf(x), Polar::ONE);
435        }
436    }
437}
abels_complex/complex/f32/rectangular.rs

abels_complex/complex/f32/
rectangular.rs
