sva_formula/
complex.rs

1// Concern: the complex amplitude every atom factor is measured in | Non-concern: what an amplitude multiplies (spectral_sum/atom.rs) | IO: (C64, C64) -> C64
2
3use std::ops::{Add, Div, Mul, Neg, Sub};
4
5#[derive(Clone, Copy, Debug)]
6pub struct C64 {
7    pub re: f64,
8    pub im: f64,
9}
10
11/// Bitwise after canonicalizing `-0.0`, so equality is transitive.
12impl PartialEq for C64 {
13    fn eq(&self, other: &C64) -> bool {
14        self.bits() == other.bits()
15    }
16}
17impl Eq for C64 {}
18
19impl C64 {
20    pub const ZERO: C64 = C64 { re: 0.0, im: 0.0 };
21    pub const ONE: C64 = C64 { re: 1.0, im: 0.0 };
22    pub const I: C64 = C64 { re: 0.0, im: 1.0 };
23
24    pub fn new(re: f64, im: f64) -> C64 {
25        C64 { re, im }
26    }
27
28    pub fn real(re: f64) -> C64 {
29        C64 { re, im: 0.0 }
30    }
31
32    pub fn bits(self) -> (u64, u64) {
33        (canonical(self.re), canonical(self.im))
34    }
35
36    pub fn is_finite(self) -> bool {
37        self.re.is_finite() && self.im.is_finite()
38    }
39
40    pub fn is_zero(self) -> bool {
41        self.bits() == C64::ZERO.bits()
42    }
43
44    pub fn is_real(self) -> bool {
45        canonical(self.im) == canonical(0.0)
46    }
47
48    pub fn conj(self) -> C64 {
49        C64::new(self.re, -self.im)
50    }
51
52    pub fn norm_sqr(self) -> f64 {
53        self.re * self.re + self.im * self.im
54    }
55
56    pub fn abs(self) -> f64 {
57        libm::hypot(self.re, self.im)
58    }
59
60    /// A component alone on its axis inverts by one division, correctly rounded; only a
61    /// value off both axes takes the conjugate over the squared norm.
62    pub fn inv(self) -> C64 {
63        let d = self.norm_sqr();
64        match (self.re == 0.0, self.im == 0.0) {
65            (false, true) => C64::new(1.0 / self.re, -self.im / d),
66            (true, false) => C64::new(self.re / d, -1.0 / self.im),
67            _ => C64::new(self.re / d, -self.im / d),
68        }
69    }
70
71    pub fn exp(self) -> C64 {
72        let r = crate::exp(self.re);
73        let (sin, cos) = libm::sincos(self.im);
74        C64::new(r * cos, r * sin)
75    }
76
77    /// `e^z - 1` without the cancellation `exp() - ONE` meets near zero.
78    pub fn exp_m1(self) -> C64 {
79        let half = libm::sin(self.im * 0.5);
80        C64::new(
81            libm::expm1(self.re) * libm::cos(self.im) - 2.0 * half * half,
82            crate::exp(self.re) * libm::sin(self.im),
83        )
84    }
85
86    pub fn powi(self, n: u32) -> C64 {
87        let mut out = C64::ONE;
88        for _ in 0..n {
89            out = out * self;
90        }
91        out
92    }
93
94    pub fn scale(self, k: f64) -> C64 {
95        C64::new(self.re * k, self.im * k)
96    }
97
98    pub fn over(self, k: f64) -> C64 {
99        C64::new(self.re / k, self.im / k)
100    }
101}
102
103pub fn canonical(x: f64) -> u64 {
104    if x == 0.0 {
105        0f64.to_bits()
106    } else {
107        x.to_bits()
108    }
109}
110
111impl Add for C64 {
112    type Output = C64;
113    fn add(self, o: C64) -> C64 {
114        C64::new(self.re + o.re, self.im + o.im)
115    }
116}
117impl Sub for C64 {
118    type Output = C64;
119    fn sub(self, o: C64) -> C64 {
120        C64::new(self.re - o.re, self.im - o.im)
121    }
122}
123impl Mul for C64 {
124    type Output = C64;
125    fn mul(self, o: C64) -> C64 {
126        C64::new(
127            self.re * o.re - self.im * o.im,
128            self.re * o.im + self.im * o.re,
129        )
130    }
131}
132impl Div for C64 {
133    type Output = C64;
134    fn div(self, o: C64) -> C64 {
135        if o.im == 0.0 && o.re != 0.0 {
136            return C64::new(self.re / o.re, self.im / o.re);
137        }
138        if o.re == 0.0 && o.im != 0.0 {
139            return C64::new(self.im / o.im, -self.re / o.im);
140        }
141        let d = o.norm_sqr();
142        C64::new(
143            (self.re * o.re + self.im * o.im) / d,
144            (self.im * o.re - self.re * o.im) / d,
145        )
146    }
147}
148impl Neg for C64 {
149    type Output = C64;
150    fn neg(self) -> C64 {
151        C64::new(-self.re, -self.im)
152    }
153}
sva_formula/complex.rs

sva_formula/
complex.rs
