abels_complex/complex/
polar.rs

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

abels_complex/complex/
polar.rs
