abels_complex/complex/
rectangular.rs

1use core::ops::*;
2pub type Complex32 = Complex<f32>;
3pub type Complex64 = Complex<f64>;
4use crate::traits::Number;
5
6use super::ComplexPolar as Polar;
7
8/// Creates a complex number in rectangular form.
9#[inline(always)]
10#[must_use]
11pub const fn complex<FT>(re: FT, im: FT) -> Complex<FT> {
12    Complex::new(re, im)
13}
14
15/// A complex number in rectangular form.
16#[derive(Clone, Copy, PartialEq, Debug, Default)]
17#[repr(C)]
18pub struct Complex<FT> {
19    pub re: FT,
20    pub im: FT,
21}
22impl<FT> Complex<FT> {
23    /// Creates a complex number.
24    pub const fn new(re: FT, im: FT) -> Self {
25        Self { re, im }
26    }
27}
28
29impl<FT: Number> Complex<FT> {
30    pub const ZERO: Self = Self::new(FT::ZERO, FT::ZERO);
31    pub const ONE: Self = Self::new(FT::ONE, FT::ZERO);
32    pub const I: Self = Self::new(FT::ZERO, FT::ONE);
33
34    /// Computes the conjugate.
35    pub fn conjugate(self) -> Self {
36        Self::new(self.re, -self.im)
37    }
38
39    /// Computes the absolute value.
40    pub fn abs(self) -> FT {
41        self.abs_sq().sqrt()
42    }
43
44    pub fn square(mut self) -> Self {
45        let two = FT::ONE + FT::ONE;
46        let re = self.re * self.re - self.im * self.im;
47        self.im = self.re * self.im * two;
48        self.re = re;
49        self
50    }
51
52    /// Computes the squared absolute value.
53    ///
54    /// This is faster than `abs()` as it avoids a square root operation.
55    pub fn abs_sq(self) -> FT {
56        self.re * self.re + self.im * self.im
57    }
58
59    /// Computes the argument in the range `(-π, +π]`.
60    pub fn arg(self) -> FT {
61        self.im.atan2(self.re)
62    }
63
64    /// Computes the reciprocal.
65    pub fn recip(self) -> Self {
66        self.conjugate() / self.abs_sq()
67    }
68
69    /// Convert to polar form.
70    pub fn to_polar(self) -> Polar<FT> {
71        Polar::new(self.abs(), self.arg())
72    }
73
74    /// Computes `e^self` where `e` is the base of the natural logarithm.
75    pub fn exp(self) -> Polar<FT> {
76        Polar::new(self.re.exp(), self.im)
77    }
78
79    /// Computes `2^self`.
80    pub fn exp2(self) -> Polar<FT> {
81        Polar::new(self.re.exp2(), self.im * FT::LN_2())
82    }
83
84    /// Computes the principal natural logarithm.
85    pub fn ln(self) -> Self {
86        self.to_polar().ln()
87    }
88
89    /// Computes the principal natural logarithm of `1 + self`.
90    ///
91    /// More numerically stable than `(self + 1).ln()` when `self ≈ 0`.
92    pub fn ln_1p(self) -> Self {
93        let two = FT::ONE + FT::ONE;
94        let re = (two * self.re + self.abs_sq()).ln_1p() / two;
95        let im = (self + FT::ONE).arg();
96        Self::new(re, im)
97    }
98
99    /// Computes the principal logarithm in base 2.
100    pub fn log2(self) -> Self {
101        self.ln() / FT::LN_2()
102    }
103
104    /// Computes the principal logarithm in base 10.
105    pub fn log10(self) -> Self {
106        self.ln() / FT::LN_10()
107    }
108
109    /// Computes the k-th branch of the natural logarithm.
110    ///
111    /// The principal value is `k = 0`. Each increment of `k` adds `2πi`.
112    pub fn ln_branch(self, k: i32) -> Self {
113        self.to_polar().ln_branch(k)
114    }
115
116    /// Computes the k-th branch of the natural logarithm of `1 + self`.
117    ///
118    /// The principal value is `k = 0`. Each increment of `k` adds `2πi`.
119    pub fn ln_1p_branch(self, k: i32) -> Self {
120        let p = self.ln_1p();
121        Self::new(p.re, p.im + FT::TAU() * FT::from_i32(k))
122    }
123
124    /// Computes the k-th branch of the base-2 logarithm.
125    pub fn log2_branch(self, k: i32) -> Self {
126        self.to_polar().log2_branch(k)
127    }
128
129    /// Computes the k-th branch of the base-10 logarithm.
130    pub fn log10_branch(self, k: i32) -> Self {
131        self.to_polar().log10_branch(k)
132    }
133
134    /// Computes the k-th square root.
135    ///
136    /// The principal value is `k = 0`. Only `k = 0` and `k = 1` give distinct values.
137    pub fn sqrt_branch(self, k: i32) -> Polar<FT> {
138        self.to_polar().sqrt_branch(k)
139    }
140
141    /// Computes the k-th value of the n-th root.
142    ///
143    /// The `n` distinct values correspond to `k = 0..n-1`.
144    pub fn nth_root(self, n: i32, k: i32) -> Polar<FT> {
145        self.to_polar().nth_root(n, k)
146    }
147
148    /// Raises `self` to the rational power `p/q`, selecting the k-th branch.
149    ///
150    /// There are `q` distinct values corresponding to `k = 0..q-1`, provided `p/q` is in lowest
151    /// terms. If `p/q` is not reduced, first reduce it to find the true number of distinct values.
152    pub fn pow_rational(self, p: i32, q: i32, k: i32) -> Polar<FT> {
153        self.to_polar().pow_rational(p, q, k)
154    }
155
156    /// Raises `self` to an integer power.
157    pub fn powi(self, n: i32) -> Polar<FT> {
158        self.to_polar().powi(n)
159    }
160
161    /// Raises `self` to a floating point power.
162    pub fn powf(self, x: FT) -> Polar<FT> {
163        self.to_polar().powf(x)
164    }
165
166    /// Computes the principal square root.
167    pub fn sqrt(self) -> Self {
168        let two = FT::ONE + FT::ONE;
169        let abs = self.abs();
170        Self::new(
171            ((abs + self.re) / two).sqrt(),
172            ((abs - self.re) / two).sqrt().copysign(self.im),
173        )
174    }
175
176    /// Computes the euclidian distance between two points.
177    pub fn distance(self, other: Self) -> FT {
178        (self - other).abs()
179    }
180
181    /// Computes the squared euclidian distance between two points.
182    pub fn distance_squared(self, other: Self) -> FT {
183        (self - other).abs_sq()
184    }
185
186    /// Computes the linear interpolation between two points based on the value `t`.
187    pub fn lerp(self, other: Self, t: FT) -> Self {
188        self + (other - self) * t
189    }
190}
191
192impl<FT: Number> Add for Complex<FT> {
193    type Output = Self;
194    fn add(self, other: Self) -> Self::Output {
195        Complex::new(self.re + other.re, self.im + other.im)
196    }
197}
198
199impl<FT: Number> Add<FT> for Complex<FT> {
200    type Output = Self;
201    fn add(self, re: FT) -> Self::Output {
202        Complex::new(self.re + re, self.im)
203    }
204}
205
206impl<FT: Number> AddAssign for Complex<FT> {
207    fn add_assign(&mut self, other: Self) {
208        self.re += other.re;
209        self.im += other.im;
210    }
211}
212
213impl<FT: Number> AddAssign<FT> for Complex<FT> {
214    fn add_assign(&mut self, re: FT) {
215        self.re += re;
216    }
217}
218
219impl<FT: Number> Sub for Complex<FT> {
220    type Output = Self;
221    fn sub(self, other: Self) -> Self::Output {
222        Complex::new(self.re - other.re, self.im - other.im)
223    }
224}
225
226impl<FT: Number> Sub<FT> for Complex<FT> {
227    type Output = Self;
228    fn sub(self, re: FT) -> Self::Output {
229        Complex::new(self.re - re, self.im)
230    }
231}
232
233impl<FT: Number> SubAssign for Complex<FT> {
234    fn sub_assign(&mut self, other: Self) {
235        self.re -= other.re;
236        self.im -= other.im;
237    }
238}
239
240impl<FT: Number> SubAssign<FT> for Complex<FT> {
241    fn sub_assign(&mut self, re: FT) {
242        self.re -= re;
243    }
244}
245
246impl<FT: Number> Mul for Complex<FT> {
247    type Output = Self;
248    fn mul(mut self, other: Self) -> Self {
249        self *= other;
250        self
251    }
252}
253
254impl<FT: Number> Mul<FT> for Complex<FT> {
255    type Output = Self;
256    fn mul(self, re: FT) -> Self {
257        Complex::new(self.re * re, self.im * re)
258    }
259}
260
261impl<FT: Number> MulAssign for Complex<FT> {
262    fn mul_assign(&mut self, other: Self) {
263        let re = self.re * other.re - self.im * other.im;
264        self.im = self.re * other.im + self.im * other.re;
265        self.re = re;
266    }
267}
268
269impl<FT: Number> MulAssign<FT> for Complex<FT> {
270    fn mul_assign(&mut self, re: FT) {
271        self.re *= re;
272        self.im *= re;
273    }
274}
275
276impl<FT: Number> Div for Complex<FT> {
277    type Output = Self;
278    fn div(self, other: Self) -> Self::Output {
279        self * other.recip()
280    }
281}
282
283impl<FT: Number> Div<FT> for Complex<FT> {
284    type Output = Self;
285    fn div(self, re: FT) -> Self::Output {
286        Complex::new(self.re / re, self.im / re)
287    }
288}
289
290impl<FT: Number> DivAssign for Complex<FT> {
291    fn div_assign(&mut self, other: Self) {
292        *self = *self / other;
293    }
294}
295
296impl<FT: Number> DivAssign<FT> for Complex<FT> {
297    fn div_assign(&mut self, re: FT) {
298        self.re /= re;
299        self.im /= re;
300    }
301}
302
303impl<FT: Number> Neg for Complex<FT> {
304    type Output = Self;
305    fn neg(self) -> Self::Output {
306        Self::new(-self.re, -self.im)
307    }
308}
309
310impl<FT: Number> From<FT> for Complex<FT> {
311    fn from(value: FT) -> Self {
312        Self::new(value, FT::ZERO)
313    }
314}
315
316#[cfg(feature = "approx")]
317use approx::{AbsDiffEq, RelativeEq, UlpsEq};
318
319#[cfg(feature = "approx")]
320impl<FT: AbsDiffEq + Copy> AbsDiffEq for Complex<FT>
321where
322    <FT as AbsDiffEq>::Epsilon: Copy,
323{
324    type Epsilon = <FT as AbsDiffEq>::Epsilon;
325    fn default_epsilon() -> Self::Epsilon {
326        FT::default_epsilon()
327    }
328    fn abs_diff_eq(&self, other: &Self, epsilon: Self::Epsilon) -> bool {
329        FT::abs_diff_eq(&self.re, &other.re, epsilon)
330            && FT::abs_diff_eq(&self.im, &other.im, epsilon)
331    }
332}
333
334#[cfg(feature = "approx")]
335impl<FT: RelativeEq + Copy> RelativeEq for Complex<FT>
336where
337    <FT as AbsDiffEq>::Epsilon: Copy,
338{
339    fn default_max_relative() -> Self::Epsilon {
340        FT::default_max_relative()
341    }
342    fn relative_eq(
343        &self,
344        other: &Self,
345        epsilon: Self::Epsilon,
346        max_relative: Self::Epsilon,
347    ) -> bool {
348        FT::relative_eq(&self.re, &other.re, epsilon, max_relative)
349            && FT::relative_eq(&self.im, &other.im, epsilon, max_relative)
350    }
351}
352
353#[cfg(feature = "approx")]
354impl<FT: UlpsEq + Copy> UlpsEq for Complex<FT>
355where
356    <FT as AbsDiffEq>::Epsilon: Copy,
357{
358    fn default_max_ulps() -> u32 {
359        FT::default_max_ulps()
360    }
361    fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool {
362        FT::ulps_eq(&self.re, &other.re, epsilon, max_ulps)
363            && FT::ulps_eq(&self.im, &other.im, epsilon, max_ulps)
364    }
365}
abels_complex/complex/rectangular.rs

abels_complex/complex/
rectangular.rs
