1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
//! `num_complex::Complex` transcendentals evaluated through `SemiflowFloat::libm_*`.
//!
//! `num-complex` routes its `Float`-based methods through `num-traits`, whose
//! backend is the platform math library in `std` builds and `libm` in `no_std`
//! builds. These ports keep num-complex 0.4.6's formulas and special cases
//! verbatim but always use `libm`, so `std` and `no_std` produce identical bits
//! (ADR-0200). `clippy::disallowed_methods` forbids the num-complex originals.
use num_complex::Complex;
use crate::float::SemiflowFloat;
/// `|z|` (num-complex `Complex::norm`).
#[inline]
pub(crate) fn norm<T: SemiflowFloat>(z: Complex<T>) -> T {
z.re.libm_hypot(z.im)
}
/// Principal argument `arg z ∈ (−π, π]` (num-complex `Complex::arg`).
#[inline]
pub(crate) fn arg<T: SemiflowFloat>(z: Complex<T>) -> T {
z.im.libm_atan2(z.re)
}
/// `(|z|, arg z)` (num-complex `Complex::to_polar`).
#[inline]
pub(crate) fn to_polar<T: SemiflowFloat>(z: Complex<T>) -> (T, T) {
(norm(z), arg(z))
}
/// `r·e^{iθ}` (num-complex `Complex::from_polar`).
#[inline]
pub(crate) fn from_polar<T: SemiflowFloat>(r: T, theta: T) -> Complex<T> {
Complex::new(r * theta.libm_cos(), r * theta.libm_sin())
}
/// `e^z` (num-complex `Complex::exp`, including its ±∞/NaN corner cases).
#[inline]
pub(crate) fn exp<T: SemiflowFloat>(z: Complex<T>) -> Complex<T> {
let Complex { re, mut im } = z;
if re.is_infinite() {
if re < T::zero() {
if !im.is_finite() {
return Complex::new(T::zero(), T::zero());
}
} else if im == T::zero() || !im.is_finite() {
if im.is_infinite() {
im = T::nan();
}
return Complex::new(re, im);
}
} else if re.is_nan() && im == T::zero() {
return z;
}
from_polar(re.libm_exp(), im)
}
/// Principal square root, branch cut `(−∞, 0)` continuous from above
/// (num-complex `Complex::sqrt`).
#[inline]
pub(crate) fn sqrt<T: SemiflowFloat>(z: Complex<T>) -> Complex<T> {
let two = T::one() + T::one();
if z.im.is_zero() {
if z.re.is_sign_positive() {
Complex::new(z.re.sqrt(), z.im)
} else {
let re = T::zero();
let im = (-z.re).sqrt();
if z.im.is_sign_positive() {
Complex::new(re, im)
} else {
Complex::new(re, -im)
}
}
} else if z.re.is_zero() {
let x = (z.im.abs() / two).sqrt();
if z.im.is_sign_positive() {
Complex::new(x, x)
} else {
Complex::new(x, -x)
}
} else {
let (r, theta) = to_polar(z);
from_polar(r.sqrt(), theta / two)
}
}