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

abels_complex/complex/f64/
rectangular.rs
