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    pub fn inv(self) -> C64 {
61        let d = self.norm_sqr();
62        C64::new(self.re / d, -self.im / d)
63    }
64
65    pub fn exp(self) -> C64 {
66        let r = self.re.exp();
67        C64::new(r * self.im.cos(), r * self.im.sin())
68    }
69
70    pub fn powi(self, n: u32) -> C64 {
71        let mut out = C64::ONE;
72        for _ in 0..n {
73            out = out * self;
74        }
75        out
76    }
77
78    pub fn scale(self, k: f64) -> C64 {
79        C64::new(self.re * k, self.im * k)
80    }
81
82    pub fn over(self, k: f64) -> C64 {
83        C64::new(self.re / k, self.im / k)
84    }
85}
86
87pub fn canonical(x: f64) -> u64 {
88    if x == 0.0 {
89        0f64.to_bits()
90    } else {
91        x.to_bits()
92    }
93}
94
95impl Add for C64 {
96    type Output = C64;
97    fn add(self, o: C64) -> C64 {
98        C64::new(self.re + o.re, self.im + o.im)
99    }
100}
101impl Sub for C64 {
102    type Output = C64;
103    fn sub(self, o: C64) -> C64 {
104        C64::new(self.re - o.re, self.im - o.im)
105    }
106}
107impl Mul for C64 {
108    type Output = C64;
109    fn mul(self, o: C64) -> C64 {
110        C64::new(
111            self.re * o.re - self.im * o.im,
112            self.re * o.im + self.im * o.re,
113        )
114    }
115}
116impl Div for C64 {
117    type Output = C64;
118    fn div(self, o: C64) -> C64 {
119        let d = o.norm_sqr();
120        C64::new(
121            (self.re * o.re + self.im * o.im) / d,
122            (self.im * o.re - self.re * o.im) / d,
123        )
124    }
125}
126impl Neg for C64 {
127    type Output = C64;
128    fn neg(self) -> C64 {
129        C64::new(-self.re, -self.im)
130    }
131}
sva_formula/complex.rs

sva_formula/
complex.rs
