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 `e^self - 1`.
85    ///
86    /// More numerically stable than `self.exp().to_rectangular() - 1` when `self ≈ 0`.
87    pub fn expm1(self) -> Self {
88        let two = FT::ONE + FT::ONE;
89        let (sin, cos) = self.im.sin_cos();
90        Self::new(
91            self.re.exp_m1() * cos - two * (self.im / two).sin().powi(2),
92            self.re.exp() * sin,
93        )
94    }
95
96    /// Computes the principal natural logarithm.
97    pub fn ln(self) -> Self {
98        self.to_polar().ln()
99    }
100
101    /// Computes the principal natural logarithm of `1 + self`.
102    ///
103    /// More numerically stable than `(self + 1).ln()` when `self ≈ 0`.
104    pub fn ln_1p(self) -> Self {
105        let two = FT::ONE + FT::ONE;
106        let re = (two * self.re + self.abs_sq()).ln_1p() / two;
107        let im = (self + FT::ONE).arg();
108        Self::new(re, im)
109    }
110
111    /// Computes the principal logarithm in base 2.
112    pub fn log2(self) -> Self {
113        self.ln() / FT::LN_2()
114    }
115
116    /// Computes the principal logarithm in base 10.
117    pub fn log10(self) -> Self {
118        self.ln() / FT::LN_10()
119    }
120
121    /// Computes the k-th branch of the natural logarithm.
122    ///
123    /// The principal value is `k = 0`. Each increment of `k` adds `2πi`.
124    pub fn ln_branch(self, k: i32) -> Self {
125        self.to_polar().ln_branch(k)
126    }
127
128    /// Computes the k-th branch of the natural logarithm of `1 + self`.
129    ///
130    /// The principal value is `k = 0`. Each increment of `k` adds `2πi`.
131    pub fn ln_1p_branch(self, k: i32) -> Self {
132        let p = self.ln_1p();
133        Self::new(p.re, p.im + FT::TAU() * FT::from_i32(k))
134    }
135
136    /// Computes the k-th branch of the base-2 logarithm.
137    pub fn log2_branch(self, k: i32) -> Self {
138        self.to_polar().log2_branch(k)
139    }
140
141    /// Computes the k-th branch of the base-10 logarithm.
142    pub fn log10_branch(self, k: i32) -> Self {
143        self.to_polar().log10_branch(k)
144    }
145
146    /// Computes the k-th square root.
147    ///
148    /// The principal value is `k = 0`. Only `k = 0` and `k = 1` give distinct values.
149    pub fn sqrt_branch(self, k: i32) -> Polar<FT> {
150        self.to_polar().sqrt_branch(k)
151    }
152
153    /// Computes the k-th value of the n-th root.
154    ///
155    /// The `n` distinct values correspond to `k = 0..n-1`.
156    pub fn nth_root(self, n: i32, k: i32) -> Polar<FT> {
157        self.to_polar().nth_root(n, k)
158    }
159
160    /// Raises `self` to the rational power `p/q`, selecting the k-th branch.
161    ///
162    /// There are `q` distinct values corresponding to `k = 0..q-1`, provided `p/q` is in lowest
163    /// terms. If `p/q` is not reduced, first reduce it to find the true number of distinct values.
164    pub fn pow_rational(self, p: i32, q: i32, k: i32) -> Polar<FT> {
165        self.to_polar().pow_rational(p, q, k)
166    }
167
168    /// Raises `self` to an integer power.
169    pub fn powi(self, n: i32) -> Polar<FT> {
170        self.to_polar().powi(n)
171    }
172
173    /// Raises `self` to a floating point power.
174    pub fn powf(self, x: FT) -> Polar<FT> {
175        self.to_polar().powf(x)
176    }
177
178    /// Computes the principal square root.
179    pub fn sqrt(self) -> Self {
180        let two = FT::ONE + FT::ONE;
181        let abs = self.abs();
182        Self::new(
183            ((abs + self.re) / two).sqrt(),
184            ((abs - self.re) / two).sqrt().copysign(self.im),
185        )
186    }
187
188    /// Computes the euclidian distance between two points.
189    pub fn distance(self, other: Self) -> FT {
190        (self - other).abs()
191    }
192
193    /// Computes the squared euclidian distance between two points.
194    pub fn distance_squared(self, other: Self) -> FT {
195        (self - other).abs_sq()
196    }
197
198    /// Computes the linear interpolation between two points based on the value `t`.
199    pub fn lerp(self, other: Self, t: FT) -> Self {
200        self + (other - self) * t
201    }
202
203    /// Returns whether both the real and imaginary component are finite.
204    pub fn is_finite(self) -> bool {
205        self.re.is_finite() && self.im.is_finite()
206    }
207}
208
209impl<FT: Number> Add for Complex<FT> {
210    type Output = Self;
211    fn add(self, other: Self) -> Self::Output {
212        Complex::new(self.re + other.re, self.im + other.im)
213    }
214}
215
216impl<FT: Number> Add<FT> for Complex<FT> {
217    type Output = Self;
218    fn add(self, re: FT) -> Self::Output {
219        Complex::new(self.re + re, self.im)
220    }
221}
222
223impl<FT: Number> AddAssign for Complex<FT> {
224    fn add_assign(&mut self, other: Self) {
225        self.re += other.re;
226        self.im += other.im;
227    }
228}
229
230impl<FT: Number> AddAssign<FT> for Complex<FT> {
231    fn add_assign(&mut self, re: FT) {
232        self.re += re;
233    }
234}
235
236impl<FT: Number> Sub for Complex<FT> {
237    type Output = Self;
238    fn sub(self, other: Self) -> Self::Output {
239        Complex::new(self.re - other.re, self.im - other.im)
240    }
241}
242
243impl<FT: Number> Sub<FT> for Complex<FT> {
244    type Output = Self;
245    fn sub(self, re: FT) -> Self::Output {
246        Complex::new(self.re - re, self.im)
247    }
248}
249
250impl<FT: Number> SubAssign for Complex<FT> {
251    fn sub_assign(&mut self, other: Self) {
252        self.re -= other.re;
253        self.im -= other.im;
254    }
255}
256
257impl<FT: Number> SubAssign<FT> for Complex<FT> {
258    fn sub_assign(&mut self, re: FT) {
259        self.re -= re;
260    }
261}
262
263impl<FT: Number> Mul for Complex<FT> {
264    type Output = Self;
265    fn mul(mut self, other: Self) -> Self {
266        self *= other;
267        self
268    }
269}
270
271impl<FT: Number> Mul<FT> for Complex<FT> {
272    type Output = Self;
273    fn mul(self, re: FT) -> Self {
274        Complex::new(self.re * re, self.im * re)
275    }
276}
277
278impl<FT: Number> MulAssign for Complex<FT> {
279    fn mul_assign(&mut self, other: Self) {
280        let re = self.re * other.re - self.im * other.im;
281        self.im = self.re * other.im + self.im * other.re;
282        self.re = re;
283    }
284}
285
286impl<FT: Number> MulAssign<FT> for Complex<FT> {
287    fn mul_assign(&mut self, re: FT) {
288        self.re *= re;
289        self.im *= re;
290    }
291}
292
293impl<FT: Number> Div for Complex<FT> {
294    type Output = Self;
295    fn div(self, other: Self) -> Self::Output {
296        self * other.recip()
297    }
298}
299
300impl<FT: Number> Div<FT> for Complex<FT> {
301    type Output = Self;
302    fn div(self, re: FT) -> Self::Output {
303        Complex::new(self.re / re, self.im / re)
304    }
305}
306
307impl<FT: Number> DivAssign for Complex<FT> {
308    fn div_assign(&mut self, other: Self) {
309        *self = *self / other;
310    }
311}
312
313impl<FT: Number> DivAssign<FT> for Complex<FT> {
314    fn div_assign(&mut self, re: FT) {
315        self.re /= re;
316        self.im /= re;
317    }
318}
319
320impl<FT: Number> Neg for Complex<FT> {
321    type Output = Self;
322    fn neg(self) -> Self::Output {
323        Self::new(-self.re, -self.im)
324    }
325}
326
327impl<FT: Number> From<FT> for Complex<FT> {
328    fn from(value: FT) -> Self {
329        Self::new(value, FT::ZERO)
330    }
331}
332
333#[cfg(feature = "approx")]
334use approx::{AbsDiffEq, RelativeEq, UlpsEq};
335
336#[cfg(feature = "approx")]
337impl<FT: AbsDiffEq + Copy> AbsDiffEq for Complex<FT>
338where
339    <FT as AbsDiffEq>::Epsilon: Copy,
340{
341    type Epsilon = <FT as AbsDiffEq>::Epsilon;
342    fn default_epsilon() -> Self::Epsilon {
343        FT::default_epsilon()
344    }
345    fn abs_diff_eq(&self, other: &Self, epsilon: Self::Epsilon) -> bool {
346        FT::abs_diff_eq(&self.re, &other.re, epsilon)
347            && FT::abs_diff_eq(&self.im, &other.im, epsilon)
348    }
349}
350
351#[cfg(feature = "approx")]
352impl<FT: RelativeEq + Copy> RelativeEq for Complex<FT>
353where
354    <FT as AbsDiffEq>::Epsilon: Copy,
355{
356    fn default_max_relative() -> Self::Epsilon {
357        FT::default_max_relative()
358    }
359    fn relative_eq(
360        &self,
361        other: &Self,
362        epsilon: Self::Epsilon,
363        max_relative: Self::Epsilon,
364    ) -> bool {
365        FT::relative_eq(&self.re, &other.re, epsilon, max_relative)
366            && FT::relative_eq(&self.im, &other.im, epsilon, max_relative)
367    }
368}
369
370#[cfg(feature = "approx")]
371impl<FT: UlpsEq + Copy> UlpsEq for Complex<FT>
372where
373    <FT as AbsDiffEq>::Epsilon: Copy,
374{
375    fn default_max_ulps() -> u32 {
376        FT::default_max_ulps()
377    }
378    fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool {
379        FT::ulps_eq(&self.re, &other.re, epsilon, max_ulps)
380            && FT::ulps_eq(&self.im, &other.im, epsilon, max_ulps)
381    }
382}
383
384impl<FT: num_traits::AsPrimitive<f32>> Complex<FT> {
385    /// Casts to a glam::Vec2.
386    #[cfg(feature = "glam")]
387    pub fn as_vec2(self) -> glam::Vec2 {
388        glam::vec2(self.re.as_(), self.im.as_())
389    }
390}
391
392impl<FT: num_traits::AsPrimitive<f64>> Complex<FT> {
393    /// Casts to a glam::DVec2.
394    #[cfg(feature = "glam")]
395    pub fn as_dvec2(self) -> glam::DVec2 {
396        glam::dvec2(self.re.as_(), self.im.as_())
397    }
398}
399
400#[cfg(feature = "glam")]
401impl<T: From<f32>> From<glam::Vec2> for Complex<T> {
402    fn from(v: glam::Vec2) -> Self {
403        Self::new(v.x.into(), v.y.into())
404    }
405}
406
407#[cfg(feature = "glam")]
408impl<T: Into<f32>> From<Complex<T>> for glam::Vec2 {
409    fn from(z: Complex<T>) -> Self {
410        glam::vec2(z.re.into(), z.im.into())
411    }
412}
413
414#[cfg(feature = "glam")]
415impl<T: From<f64>> From<glam::DVec2> for Complex<T> {
416    fn from(v: glam::DVec2) -> Self {
417        Self::new(v.x.into(), v.y.into())
418    }
419}
420
421#[cfg(feature = "glam")]
422impl<T: Into<f64>> From<Complex<T>> for glam::DVec2 {
423    fn from(z: Complex<T>) -> Self {
424        glam::dvec2(z.re.into(), z.im.into())
425    }
426}
abels_complex/complex/rectangular.rs

abels_complex/complex/
rectangular.rs
