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    /// Tolerance: max(MACH_EPSILON, 1e-18).
26    fn tol() -> Self;
27    /// Decimal digits: log10(2) * (DIGITS - 1).
28    fn dig() -> Self;
29    /// Large order threshold: 10 + 6*(DIG - 3).
30    fn fnul() -> Self;
31    /// Asymptotic region boundary: 1.2*DIG + 3.
32    fn rl() -> Self;
33    /// Underflow elimination threshold: 2.303*(K*R1M5 - 3), K = min(|MIN_EXP|, MAX_EXP).
34    fn elim() -> Self;
35    /// Overflow elimination threshold: ELIM + max(-2.303*R1M5*(DIGITS-1), -41.45).
36    fn alim() -> Self;
37}
38
39impl BesselFloat for f64 {
40    const MACH_EPSILON: f64 = 2.220446049250313e-16;
41    const MACH_TINY: f64 = 2.2250738585072014e-308;
42    const MACH_HUGE: f64 = 1.7976931348623157e+308;
43    const MACH_DIGITS: i32 = 53;
44    const MACH_MIN_EXP: i32 = -1021;
45    const MACH_MAX_EXP: i32 = 1024;
46
47    #[inline]
48    fn tol() -> f64 { 2.220446049250313e-16 }       // max(EPSILON, 1e-18)
49    #[inline]
50    fn dig() -> f64 { 15.653559774527022 }           // R1M5 * (DIGITS - 1)
51    #[inline]
52    fn fnul() -> f64 { 85.92135864716212 }           // 10 + 6*(DIG - 3)
53    #[inline]
54    fn rl() -> f64 { 21.784271729432426 }            // 1.2*DIG + 3
55    #[inline]
56    fn elim() -> f64 { 700.9217936944459 }           // 2.303*(K*R1M5 - 3)
57    #[inline]
58    fn alim() -> f64 { 664.8716455337102 }           // ELIM + max(-2.303*R1M5*(DIGITS-1), -41.45)
59}
60
61impl BesselFloat for f32 {
62    const MACH_EPSILON: f32 = 1.1920929e-7;
63    const MACH_TINY: f32 = 1.1754944e-38;
64    const MACH_HUGE: f32 = 3.4028235e+38;
65    const MACH_DIGITS: i32 = 24;
66    const MACH_MIN_EXP: i32 = -125;
67    const MACH_MAX_EXP: i32 = 128;
68
69    #[inline]
70    fn tol() -> f32 { 1.1920929e-7 }                // max(EPSILON, 1e-18)
71    #[inline]
72    fn dig() -> f32 { 6.923689900271568 }            // R1M5 * (DIGITS - 1)
73    #[inline]
74    fn fnul() -> f32 { 33.542139401629406 }          // 10 + 6*(DIG - 3)
75    #[inline]
76    fn rl() -> f32 { 11.308427880325882 }            // 1.2*DIG + 3
77    #[inline]
78    fn elim() -> f32 { 79.75001000176859 }           // 2.303*(K*R1M5 - 3)
79    #[inline]
80    fn alim() -> f32 { 63.80475216144317 }           // ELIM + max(-2.303*R1M5*(DIGITS-1), -41.45)
81}
complex_bessel/machine.rs

complex_bessel/
machine.rs
