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 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
44impl BesselFloat for f64 {
45    const MACH_EPSILON: f64 = 2.220446049250313e-16;
46    const MACH_TINY: f64 = 2.2250738585072014e-308;
47    const MACH_HUGE: f64 = 1.7976931348623157e+308;
48    const MACH_DIGITS: i32 = 53;
49    const MACH_MIN_EXP: i32 = -1021;
50    const MACH_MAX_EXP: i32 = 1024;
51
52    #[inline]
53    fn from_f64(x: f64) -> f64 {
54        x
55    }
56    #[inline]
57    fn tol() -> f64 {
58        2.220446049250313e-16
59    } // max(EPSILON, 1e-18)
60    #[inline]
61    fn fnul() -> f64 {
62        85.92135864716212
63    } // 10 + 6*(DIG - 3)
64    #[inline]
65    fn rl() -> f64 {
66        21.784271729432426
67    } // 1.2*DIG + 3
68    #[inline]
69    fn elim() -> f64 {
70        700.9217936944459
71    } // 2.303*(K*R1M5 - 3)
72    #[inline]
73    fn alim() -> f64 {
74        664.8716455337102
75    } // ELIM + max(-2.303*R1M5*(DIGITS-1), -41.45)
76}
77
78// Derived constants are written at full f64 precision to document the exact
79// formula results; the compiler rounds to f32 at compile time.
80#[allow(clippy::excessive_precision)]
81impl BesselFloat for f32 {
82    const MACH_EPSILON: f32 = 1.1920929e-7;
83    const MACH_TINY: f32 = 1.1754944e-38;
84    const MACH_HUGE: f32 = 3.4028235e+38;
85    const MACH_DIGITS: i32 = 24;
86    const MACH_MIN_EXP: i32 = -125;
87    const MACH_MAX_EXP: i32 = 128;
88
89    #[inline]
90    fn from_f64(x: f64) -> f32 {
91        x as f32
92    }
93    #[inline]
94    fn tol() -> f32 {
95        1.1920929e-7
96    } // max(EPSILON, 1e-18)
97    #[inline]
98    fn fnul() -> f32 {
99        33.542139401629406
100    } // 10 + 6*(DIG - 3)
101    #[inline]
102    fn rl() -> f32 {
103        11.308427880325882
104    } // 1.2*DIG + 3
105    #[inline]
106    fn elim() -> f32 {
107        79.75001000176859
108    } // 2.303*(K*R1M5 - 3)
109    #[inline]
110    fn alim() -> f32 {
111        63.80475216144317
112    } // ELIM + max(-2.303*R1M5*(DIGITS-1), -41.45)
113}
complex_bessel/machine.rs

complex_bessel/
machine.rs
