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    /// Casts to a glam::Vec2.
43    #[cfg(feature = "glam")]
44    pub fn as_vec2(self) -> glam::Vec2 {
45        glam::vec2(self.re, self.im)
46    }
47}
48
49impl fmt::Display for Rectangular {
50    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
51        fn fmt_x(f: &mut fmt::Formatter, x: FT, sign_plus: bool) -> fmt::Result {
52            match (f.precision(), sign_plus) {
53                (None, false) => write!(f, "{}", x),
54                (None, true) => write!(f, "{:+}", x),
55                (Some(p), false) => write!(f, "{:.*}", p, x),
56                (Some(p), true) => write!(f, "{:+.*}", p, x),
57            }
58        }
59        match (self.re, self.im, f.sign_plus()) {
60            (re, 0.0, sp) => fmt_x(f, re, sp),
61            (0.0, 1.0, false) => write!(f, "i"),
62            (0.0, 1.0, true) => write!(f, "+i"),
63            (0.0, -1.0, _) => write!(f, "-i"),
64            (0.0, im, sp) => {
65                fmt_x(f, im, sp)?;
66                write!(f, "i")
67            }
68            (re, 1.0, sp) => {
69                fmt_x(f, re, sp)?;
70                write!(f, "+i")
71            }
72            (re, -1.0, sp) => {
73                fmt_x(f, re, sp)?;
74                write!(f, "-i")
75            }
76            (re, im, sp) => {
77                fmt_x(f, re, sp)?;
78                fmt_x(f, im, true)?;
79                write!(f, "i")
80            }
81        }
82    }
83}
84
85#[cfg(feature = "rand")]
86impl rand::distr::Distribution<Rectangular> for rand::distr::StandardUniform {
87    fn sample<R: rand::Rng + ?Sized>(&self, rng: &mut R) -> Rectangular {
88        use rand::RngExt;
89        rng.sample::<Polar, _>(self).to_rectangular()
90    }
91}
92
93#[cfg(feature = "glam")]
94impl From<glam::Vec2> for Rectangular {
95    fn from(v: glam::Vec2) -> Self {
96        Self::new(v.x, v.y)
97    }
98}
99
100#[cfg(feature = "glam")]
101impl From<Rectangular> for glam::Vec2 {
102    fn from(z: Rectangular) -> Self {
103        z.as_vec2()
104    }
105}
106
107#[cfg(test)]
108mod tests {
109    use super::*;
110    use approx::*;
111    use core::f32::consts::{E, FRAC_PI_2, FRAC_PI_4, LN_2, PI, SQRT_2};
112    use core::iter::Iterator;
113    use rand::{
114        RngExt, SeedableRng,
115        distr::{Distribution, StandardUniform},
116        rngs::StdRng,
117    };
118
119    const NUM_SAMPLES: usize = 100;
120    const SEED: u64 = 21;
121
122    fn random_samples<T>() -> impl Iterator<Item = T>
123    where
124        StandardUniform: Distribution<T>,
125    {
126        StdRng::seed_from_u64(SEED)
127            .sample_iter(StandardUniform)
128            .take(NUM_SAMPLES)
129    }
130
131    #[test]
132    fn addition() {
133        for z0 in random_samples::<Rectangular>() {
134            for z1 in random_samples::<Rectangular>() {
135                let z = z0 + z1;
136                assert_eq!(z.re, z0.re + z1.re);
137                assert_eq!(z.im, z0.im + z1.im);
138
139                let z = z0 + z1.re;
140                assert_eq!(z.re, z0.re + z1.re);
141                assert_eq!(z.im, z0.im);
142
143                let z = z0.re + z1;
144                assert_eq!(z.re, z0.re + z1.re);
145                assert_eq!(z.im, z1.im);
146
147                let mut z = z0;
148                z += z1;
149                assert_eq!(z, z0 + z1);
150
151                let mut z = z0;
152                z += z1.re;
153                assert_eq!(z, z0 + z1.re);
154            }
155            assert_eq!(z0 + Rectangular::ZERO, z0);
156        }
157    }
158
159    #[test]
160    fn subtraction() {
161        for z0 in random_samples::<Rectangular>() {
162            for z1 in random_samples::<Rectangular>() {
163                let z = z0 - z1;
164                assert_eq!(z.re, z0.re - z1.re);
165                assert_eq!(z.im, z0.im - z1.im);
166
167                let z = z0 - z1.re;
168                assert_eq!(z.re, z0.re - z1.re);
169                assert_eq!(z.im, z0.im);
170
171                let z = z0.re - z1;
172                assert_eq!(z.re, z0.re - z1.re);
173                assert_eq!(z.im, -z1.im);
174
175                let mut z = z0;
176                z -= z1;
177                assert_eq!(z, z0 - z1);
178
179                let mut z = z0;
180                z -= z1.re;
181                assert_eq!(z, z0 - z1.re);
182            }
183            assert_eq!(z0 - z0, Rectangular::ZERO);
184            assert_eq!(z0 - Rectangular::ZERO, z0);
185        }
186    }
187
188    #[test]
189    fn multiplication() {
190        for z0 in random_samples::<Rectangular>() {
191            for z1 in random_samples::<Rectangular>() {
192                let z = z0 * z1;
193                assert_ulps_eq!(z.abs(), z0.abs() * z1.abs());
194                assert_ulps_eq!(
195                    z.arg().sin(),
196                    (z0.arg() + z1.arg()).sin(),
197                    epsilon = 4.0 * FT::EPSILON
198                );
199
200                let z = z0 * z1.re;
201                assert_eq!(z.re, z0.re * z1.re);
202                assert_eq!(z.im, z0.im * z1.re);
203
204                let z = z0.re * z1;
205                assert_eq!(z.re, z0.re * z1.re);
206                assert_eq!(z.im, z0.re * z1.im);
207
208                let mut z = z0;
209                z *= z1;
210                assert_eq!(z, z0 * z1);
211
212                let mut z = z0;
213                z *= z1.re;
214                assert_eq!(z, z0 * z1.re);
215            }
216            assert_eq!(z0 * Rectangular::ONE, z0);
217            assert_eq!(z0 * Rectangular::ZERO, Rectangular::ZERO);
218            assert_eq!(z0 * 0.0, Rectangular::ZERO);
219        }
220    }
221
222    #[test]
223    fn division() {
224        for z0 in random_samples::<Rectangular>() {
225            for z1 in random_samples::<Rectangular>() {
226                let z = z0 / z1;
227                assert_relative_eq!(
228                    z.abs(),
229                    z0.abs() / z1.abs(),
230                    max_relative = 3.0 * FT::EPSILON
231                );
232                assert_ulps_eq!(
233                    z.arg().sin(),
234                    (z0.arg() - z1.arg()).sin(),
235                    epsilon = 4.0 * FT::EPSILON
236                );
237
238                let z = z0 / z1.re;
239                assert_eq!(z.re, z0.re / z1.re);
240                assert_eq!(z.im, z0.im / z1.re);
241
242                let z = z0.re / z1;
243                assert_ulps_eq!(z.abs(), z0.re.abs() / z1.abs());
244                assert_ulps_eq!(
245                    z.arg().sin(),
246                    (-z0.re.signum() * z1.arg()).sin(),
247                    epsilon = 2.0 * FT::EPSILON
248                );
249
250                let mut z = z0;
251                z /= z1;
252                assert_eq!(z, z0 / z1);
253
254                let mut z = z0;
255                z /= z1.re;
256                assert_eq!(z, z0 / z1.re);
257            }
258            assert_ulps_eq!(z0 / z0, Rectangular::ONE);
259            assert_eq!(Rectangular::ZERO / z0, Rectangular::ZERO);
260        }
261    }
262
263    #[test]
264    fn negation() {
265        for z in random_samples::<Rectangular>() {
266            assert_eq!(-z, Rectangular::new(-z.re, -z.im));
267        }
268        assert_eq!(-Rectangular::ONE, Rectangular::NEG_ONE);
269        assert_eq!(-Rectangular::I, Rectangular::NEG_I);
270        assert_eq!(-Rectangular::NEG_ONE, Rectangular::ONE);
271        assert_eq!(-Rectangular::NEG_I, Rectangular::I);
272    }
273
274    #[test]
275    fn reciprocal() {
276        for z in random_samples::<Rectangular>() {
277            assert_eq!(z.recip(), 1.0 / z);
278            assert_ulps_eq!(z * z.recip(), Rectangular::ONE);
279        }
280        assert_eq!(Rectangular::ONE.recip(), Rectangular::ONE);
281        assert_eq!(Rectangular::I.recip(), Rectangular::NEG_I);
282        assert_eq!(Rectangular::NEG_ONE.recip(), Rectangular::NEG_ONE);
283        assert_eq!(Rectangular::NEG_I.recip(), Rectangular::I);
284    }
285
286    #[test]
287    fn sqrt() {
288        for z in random_samples::<Rectangular>() {
289            assert_ulps_eq!(z.sqrt().abs(), z.abs().sqrt());
290            assert_ulps_eq!(
291                z.sqrt().arg(),
292                z.arg() / 2.0,
293                epsilon = 1400.0 * FT::EPSILON
294            );
295        }
296        assert_eq!(Rectangular::ONE.sqrt(), Rectangular::ONE);
297        assert_eq!(Rectangular::NEG_ONE.sqrt(), Rectangular::I);
298        assert_eq!(
299            Rectangular::new(0.0, 2.0).sqrt(),
300            Rectangular::new(1.0, 1.0)
301        );
302        assert_eq!(
303            Rectangular::new(0.0, -2.0).sqrt(),
304            Rectangular::new(1.0, -1.0)
305        );
306    }
307
308    #[test]
309    fn abs() {
310        for z in random_samples::<Rectangular>() {
311            assert_ulps_eq!(z.abs_sq(), z.abs() * z.abs());
312            assert_eq!(z.abs_sq(), z.re * z.re + z.im * z.im);
313        }
314        assert_eq!(Rectangular::ONE.abs(), 1.0);
315        assert_eq!(Rectangular::I.abs(), 1.0);
316        assert_eq!(Rectangular::NEG_ONE.abs(), 1.0);
317        assert_eq!(Rectangular::NEG_I.abs(), 1.0);
318        assert_eq!(Rectangular::new(1.0, 1.0).abs(), SQRT_2);
319        assert_eq!(Rectangular::new(-1.0, 1.0).abs(), SQRT_2);
320        assert_eq!(Rectangular::new(-1.0, -1.0).abs(), SQRT_2);
321        assert_eq!(Rectangular::new(1.0, -1.0).abs(), SQRT_2);
322    }
323
324    #[test]
325    fn conjugate() {
326        for z in random_samples::<Rectangular>() {
327            assert_eq!(z.conjugate().re, z.re);
328            assert_eq!(z.conjugate().im, -z.im);
329            assert_eq!(z.conjugate().conjugate(), z);
330        }
331        assert_eq!(Rectangular::ONE.conjugate(), Rectangular::ONE);
332        assert_eq!(Rectangular::I.conjugate(), Rectangular::NEG_I);
333        assert_eq!(Rectangular::NEG_ONE.conjugate(), Rectangular::NEG_ONE);
334        assert_eq!(Rectangular::NEG_I.conjugate(), Rectangular::I);
335    }
336
337    #[test]
338    fn arg() {
339        assert_eq!(Rectangular::ONE.arg(), 0.0);
340        assert_eq!(Rectangular::I.arg(), FRAC_PI_2);
341        assert_eq!(Rectangular::NEG_ONE.arg(), PI);
342        assert_eq!(Rectangular::NEG_I.arg(), -FRAC_PI_2);
343    }
344
345    #[test]
346    fn exp() {
347        for z in random_samples::<Rectangular>() {
348            assert_eq!(z.exp().abs, z.re.exp());
349            assert_eq!(z.exp().arg, z.im);
350            assert_ulps_eq!(z.exp().ln(), z);
351        }
352        assert_eq!(Rectangular::ONE.exp(), Polar::new(E, 0.0));
353        assert_eq!(Rectangular::I.exp(), Polar::new(1.0, 1.0));
354        assert_eq!(Rectangular::NEG_ONE.exp(), Polar::new(E.recip(), 0.0));
355        assert_eq!(Rectangular::NEG_I.exp(), Polar::new(1.0, -1.0));
356    }
357
358    #[test]
359    fn exp2() {
360        for z in random_samples::<Rectangular>() {
361            assert_eq!(z.exp2().abs, z.re.exp2());
362            assert_eq!(z.exp2().arg, z.im * LN_2);
363            assert_ulps_eq!(z.exp2().log2(), z);
364        }
365        assert_eq!(Rectangular::ONE.exp2(), Polar::new(2.0, 0.0));
366        assert_eq!(Rectangular::I.exp2(), Polar::new(1.0, LN_2));
367        assert_eq!(Rectangular::NEG_ONE.exp2(), Polar::new(0.5, 0.0));
368        assert_eq!(Rectangular::NEG_I.exp2(), Polar::new(1.0, -LN_2));
369    }
370
371    #[test]
372    fn expm1() {
373        for z in random_samples::<Rectangular>() {
374            assert_ulps_eq!(z.expm1(), z.exp().to_rectangular() - 1.0, max_ulps = 5);
375        }
376        assert_eq!(Rectangular::ZERO.expm1(), Rectangular::ZERO);
377        assert_ulps_eq!(Rectangular::ONE.expm1(), Rectangular::new(E - 1.0, 0.0));
378        assert_ulps_eq!(
379            Rectangular::I.expm1(),
380            Rectangular::I.exp().to_rectangular() - 1.0
381        );
382
383        let tiny = Rectangular::new(1e-8, 0.0);
384        assert_eq!(tiny.exp().to_rectangular() - 1.0, Rectangular::ZERO);
385        assert_ne!(tiny.expm1(), Rectangular::ZERO);
386    }
387
388    #[test]
389    fn log() {
390        for z in random_samples::<Rectangular>() {
391            assert_eq!(z.ln().re, z.abs().ln());
392            assert_eq!(z.ln().im, z.arg());
393            assert_ulps_eq!(z.ln().exp(), z.to_polar());
394        }
395        assert_eq!(Rectangular::ONE.ln(), Rectangular::ZERO);
396        assert_eq!(Rectangular::I.ln(), Rectangular::I * FRAC_PI_2);
397        assert_eq!(Rectangular::NEG_ONE.ln(), Rectangular::I * PI);
398        assert_eq!(Rectangular::NEG_I.ln(), Rectangular::I * -FRAC_PI_2);
399
400        assert_ulps_eq!(Rectangular::new(E, 0.0).ln(), Rectangular::ONE);
401        assert_ulps_eq!(Rectangular::new(2.0, 0.0).log2(), Rectangular::ONE);
402        assert_ulps_eq!(Rectangular::new(10.0, 0.0).log10(), Rectangular::ONE);
403    }
404
405    #[test]
406    fn ln_1p() {
407        for z in random_samples::<Rectangular>() {
408            // (z + 1).ln() is unreliable near the singularity at z = -1
409            let z = Rectangular::new(z.re.max(-0.75), z.im);
410            assert_ulps_eq!(z.ln_1p(), (z + 1.0).ln(), max_ulps = 5);
411        }
412        assert_eq!(Rectangular::ZERO.ln_1p(), Rectangular::ZERO);
413        assert_ulps_eq!(Rectangular::ONE.ln_1p(), Rectangular::new(LN_2, 0.0));
414        assert_ulps_eq!(
415            Rectangular::I.ln_1p(),
416            Rectangular::new(LN_2 / 2.0, FRAC_PI_4)
417        );
418        // Near zero: (z + 1).ln() loses all precision, ln_1p preserves it
419        let tiny = Rectangular::new(1e-8, 0.0);
420        assert_eq!((tiny + 1.0).ln(), Rectangular::ZERO); // 1.0 + 1e-8 == 1.0 in f32
421        assert_ne!(tiny.ln_1p(), Rectangular::ZERO);
422    }
423
424    #[test]
425    fn powi() {
426        for z in random_samples::<Rectangular>() {
427            assert_eq!(z.powi(0), Polar::ONE);
428            assert_eq!(z.powi(1), z.to_polar());
429            for n in random_samples::<i32>() {
430                assert_eq!(z.powi(n).abs, z.abs().powi(n));
431                assert_eq!(z.powi(n).arg, z.arg() * n as FT);
432            }
433        }
434        for n in random_samples::<i32>() {
435            assert_eq!(Rectangular::ZERO.powi(n.abs()), Polar::ZERO);
436            assert_eq!(Rectangular::ONE.powi(n), Polar::ONE);
437        }
438    }
439
440    #[test]
441    fn powf() {
442        for z in random_samples::<Rectangular>() {
443            assert_eq!(z.powf(0.0), Polar::ONE);
444            assert_eq!(z.powf(1.0), z.to_polar());
445            for n in random_samples::<i32>() {
446                let x = n as FT * 0.01;
447                assert_eq!(z.powf(x).abs, z.abs().powf(x));
448                assert_eq!(z.powf(x).arg, z.arg() * x);
449            }
450        }
451        for n in random_samples::<i32>() {
452            let x = n as FT * 0.01;
453            assert_eq!(Rectangular::ZERO.powf(x.abs()), Polar::ZERO);
454            assert_eq!(Rectangular::ONE.powf(x), Polar::ONE);
455        }
456    }
457}
abels_complex/complex/f32/rectangular.rs

abels_complex/complex/f32/
rectangular.rs
