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        self.re.hypot(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 = self.re.exp();
73        C64::new(r * self.im.cos(), r * self.im.sin())
74    }
75
76    /// `e^z - 1` without the cancellation `exp() - ONE` meets near zero.
77    pub fn exp_m1(self) -> C64 {
78        let half = (self.im * 0.5).sin();
79        C64::new(
80            self.re.exp_m1() * self.im.cos() - 2.0 * half * half,
81            self.re.exp() * self.im.sin(),
82        )
83    }
84
85    pub fn powi(self, n: u32) -> C64 {
86        let mut out = C64::ONE;
87        for _ in 0..n {
88            out = out * self;
89        }
90        out
91    }
92
93    pub fn scale(self, k: f64) -> C64 {
94        C64::new(self.re * k, self.im * k)
95    }
96
97    pub fn over(self, k: f64) -> C64 {
98        C64::new(self.re / k, self.im / k)
99    }
100}
101
102pub fn canonical(x: f64) -> u64 {
103    if x == 0.0 {
104        0f64.to_bits()
105    } else {
106        x.to_bits()
107    }
108}
109
110impl Add for C64 {
111    type Output = C64;
112    fn add(self, o: C64) -> C64 {
113        C64::new(self.re + o.re, self.im + o.im)
114    }
115}
116impl Sub for C64 {
117    type Output = C64;
118    fn sub(self, o: C64) -> C64 {
119        C64::new(self.re - o.re, self.im - o.im)
120    }
121}
122impl Mul for C64 {
123    type Output = C64;
124    fn mul(self, o: C64) -> C64 {
125        C64::new(
126            self.re * o.re - self.im * o.im,
127            self.re * o.im + self.im * o.re,
128        )
129    }
130}
131impl Div for C64 {
132    type Output = C64;
133    fn div(self, o: C64) -> C64 {
134        if o.im == 0.0 && o.re != 0.0 {
135            return C64::new(self.re / o.re, self.im / o.re);
136        }
137        if o.re == 0.0 && o.im != 0.0 {
138            return C64::new(self.im / o.im, -self.re / o.im);
139        }
140        let d = o.norm_sqr();
141        C64::new(
142            (self.re * o.re + self.im * o.im) / d,
143            (self.im * o.re - self.re * o.im) / d,
144        )
145    }
146}
147impl Neg for C64 {
148    type Output = C64;
149    fn neg(self) -> C64 {
150        C64::new(-self.re, -self.im)
151    }
152}
sva_formula/complex.rs

sva_formula/
complex.rs
