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

sva_formula/
complex.rs
