complex_bessel/
machine.rs

1//! Machine constants and the `BesselFloat` trait.
2//!
3//! Constants follow the Fortran I1MACH/D1MACH conventions from TOMS 644.
4
5use num_traits::Float;
6
7/// Floating-point trait for Bessel and Airy function computation.
8///
9/// Implemented for `f64` and `f32`. Provides machine constants and
10/// derived thresholds used by the Amos algorithm.
11pub trait BesselFloat: Float + core::fmt::Debug + 'static {
12    /// Machine epsilon (D1MACH(3)).
13    const MACH_EPSILON: Self;
14    /// Smallest positive normal number (D1MACH(1)).
15    const MACH_TINY: Self;
16    /// Largest representable number (D1MACH(2)).
17    const MACH_HUGE: Self;
18    /// Number of binary digits in the mantissa (I1MACH(14)).
19    const MACH_DIGITS: i32;
20    /// Minimum binary exponent (I1MACH(12)).
21    const MACH_MIN_EXP: i32;
22    /// Maximum binary exponent (I1MACH(11)).
23    const MACH_MAX_EXP: i32;
24
25    /// Infallible conversion from f64.
26    ///
27    /// For f64 this is the identity; for f32 it truncates via `as f32`.
28    /// All Amos algorithm constants originate as f64 literals, so this
29    /// conversion always succeeds for the supported types.
30    fn from_f64(x: f64) -> Self;
31
32    /// Tolerance: max(MACH_EPSILON, 1e-18).
33    fn tol() -> Self;
34    /// Large order threshold: 10 + 6*(DIG - 3), where DIG = log10(2) * (DIGITS - 1).
35    fn fnul() -> Self;
36    /// Asymptotic region boundary: 1.2*DIG + 3, where DIG = log10(2) * (DIGITS - 1).
37    fn rl() -> Self;
38    /// Underflow elimination threshold: 2.303*(K*R1M5 - 3), K = min(|MIN_EXP|, MAX_EXP).
39    fn elim() -> Self;
40    /// Overflow elimination threshold: ELIM + max(-2.303*R1M5*(DIGITS-1), -41.45).
41    fn alim() -> Self;
42
43    /// Fused multiply-add: `self * a + b`.
44    ///
45    /// With `std` enabled, uses hardware FMA via the C library `fma()`.
46    /// Without `std`, falls back to plain `self * a + b` to avoid the
47    /// slow software FMA in libm.
48    ///
49    /// Named `fma` to avoid ambiguity with [`Float::mul_add`].
50    fn fma(self, a: Self, b: Self) -> Self;
51}
52
53impl BesselFloat for f64 {
54    const MACH_EPSILON: f64 = 2.220446049250313e-16;
55    const MACH_TINY: f64 = 2.2250738585072014e-308;
56    const MACH_HUGE: f64 = 1.7976931348623157e+308;
57    const MACH_DIGITS: i32 = 53;
58    const MACH_MIN_EXP: i32 = -1021;
59    const MACH_MAX_EXP: i32 = 1024;
60
61    #[inline]
62    fn from_f64(x: f64) -> f64 {
63        x
64    }
65    #[inline]
66    fn tol() -> f64 {
67        2.220446049250313e-16
68    } // max(EPSILON, 1e-18)
69    #[inline]
70    fn fnul() -> f64 {
71        85.92135864716212
72    } // 10 + 6*(DIG - 3)
73    #[inline]
74    fn rl() -> f64 {
75        21.784271729432426
76    } // 1.2*DIG + 3
77    #[inline]
78    fn elim() -> f64 {
79        700.9217936944459
80    } // 2.303*(K*R1M5 - 3)
81    #[inline]
82    fn alim() -> f64 {
83        664.8716455337102
84    } // ELIM + max(-2.303*R1M5*(DIGITS-1), -41.45)
85
86    #[cfg(feature = "std")]
87    #[inline]
88    fn fma(self, a: f64, b: f64) -> f64 {
89        Float::mul_add(self, a, b)
90    }
91
92    #[cfg(not(feature = "std"))]
93    #[inline]
94    fn fma(self, a: f64, b: f64) -> f64 {
95        self * a + b
96    }
97}
98
99// Derived constants are written at full f64 precision to document the exact
100// formula results; the compiler rounds to f32 at compile time.
101#[allow(clippy::excessive_precision)]
102impl BesselFloat for f32 {
103    const MACH_EPSILON: f32 = 1.1920929e-7;
104    const MACH_TINY: f32 = 1.1754944e-38;
105    const MACH_HUGE: f32 = 3.4028235e+38;
106    const MACH_DIGITS: i32 = 24;
107    const MACH_MIN_EXP: i32 = -125;
108    const MACH_MAX_EXP: i32 = 128;
109
110    #[inline]
111    fn from_f64(x: f64) -> f32 {
112        x as f32
113    }
114    #[inline]
115    fn tol() -> f32 {
116        1.1920929e-7
117    } // max(EPSILON, 1e-18)
118    #[inline]
119    fn fnul() -> f32 {
120        33.542139401629406
121    } // 10 + 6*(DIG - 3)
122    #[inline]
123    fn rl() -> f32 {
124        11.308427880325882
125    } // 1.2*DIG + 3
126    #[inline]
127    fn elim() -> f32 {
128        79.75001000176859
129    } // 2.303*(K*R1M5 - 3)
130    #[inline]
131    fn alim() -> f32 {
132        63.80475216144317
133    } // ELIM + max(-2.303*R1M5*(DIGITS-1), -41.45)
134
135    #[cfg(feature = "std")]
136    #[inline]
137    fn fma(self, a: f32, b: f32) -> f32 {
138        Float::mul_add(self, a, b)
139    }
140
141    #[cfg(not(feature = "std"))]
142    #[inline]
143    fn fma(self, a: f32, b: f32) -> f32 {
144        self * a + b
145    }
146}
complex_bessel/machine.rs

complex_bessel/
machine.rs
