sparse_ir_basis/sve/
types.rs

1//! Type definitions for SVE computation
2
3use simba::scalar::ComplexField;
4
5/// Working precision type for SVE computations
6///
7/// Values match the C-API constants defined in sparseir.h
8#[derive(Debug, Clone, Copy, PartialEq, Eq)]
9pub enum TworkType {
10    /// Use double precision (64-bit)
11    Float64 = 0, // SPIR_TWORK_FLOAT64
12    /// Use extended precision (128-bit double-double)
13    Float64X2 = 1, // SPIR_TWORK_FLOAT64X2
14    /// Automatically choose precision based on epsilon
15    Auto = -1, // SPIR_TWORK_AUTO
16}
17
18/// Determine safe epsilon and working precision
19///
20/// This function determines the safe epsilon value based on the working precision,
21/// and automatically selects the working precision if TworkType::Auto is specified.
22///
23/// # Arguments
24///
25/// * `epsilon` - Required accuracy (non-negative); `None` selects the best
26///   accuracy of the working precision
27/// * `twork` - Working precision type (Auto for automatic selection)
28///
29/// # Returns
30///
31/// Tuple of (safe_epsilon, actual_twork)
32///
33/// # Panics
34///
35/// Panics if epsilon is negative or NaN. [`compute_sve`](crate::sve::compute_sve)
36/// checks it first.
37pub(crate) fn safe_epsilon(epsilon: Option<f64>, twork: TworkType) -> (f64, TworkType) {
38    // Check for a negative or NaN epsilon (following the C++ implementation
39    // for negative values; NaN used to select the automatic accuracy)
40    if let Some(eps) = epsilon.filter(|eps| !(*eps >= 0.0)) {
41        panic!("eps_required must be non-negative, got {eps:?}");
42    }
43
44    // First, choose the working dtype based on the eps required
45    let twork_actual = match twork {
46        TworkType::Auto => match epsilon {
47            Some(eps) if eps >= 1e-8 => TworkType::Float64,
48            _ => TworkType::Float64X2, // MAX_DTYPE equivalent
49        },
50        other => other,
51    };
52
53    // Next, work out the actual epsilon.
54    // The precision floor is the smallest epsilon achievable with the chosen
55    // working type.  The returned epsilon is the *larger* of the user's
56    // request and this floor so that (a) we never promise more accuracy than
57    // the arithmetic can deliver and (b) a user who asks for *less* accuracy
58    // actually gets what they asked for.
59    let precision_floor = match twork_actual {
60        TworkType::Float64 => {
61            // This is technically a bit too low (the true value is about 1.5e-8),
62            // but it's not too far off and easier to remember for the user.
63            1e-8
64        }
65        TworkType::Float64X2 => {
66            // sqrt(Df64 epsilon) ≈ sqrt(2.465e-32) ≈ 1.57e-16
67            use crate::numeric::CustomNumeric;
68            crate::Df64::epsilon().sqrt().to_f64()
69        }
70        _ => 1e-8,
71    };
72    let safe_eps = epsilon.map_or(precision_floor, |eps| eps.max(precision_floor));
73
74    (safe_eps, twork_actual)
75}
76
77#[cfg(test)]
78mod tests {
79    use super::*;
80
81    #[test]
82    fn test_safe_epsilon_auto_float64() {
83        // epsilon=1e-7 > floor=1e-8 → safe_eps should honour the user's request
84        let (safe_eps, twork) = safe_epsilon(Some(1e-7), TworkType::Auto);
85        assert_eq!(twork, TworkType::Float64);
86        assert_eq!(safe_eps, 1e-7);
87    }
88
89    #[test]
90    fn test_safe_epsilon_auto_float64x2() {
91        // epsilon=1e-10 > floor≈1.57e-16 → safe_eps should be the user's epsilon
92        let (safe_eps, twork) = safe_epsilon(Some(1e-10), TworkType::Auto);
93        assert_eq!(twork, TworkType::Float64X2);
94        assert_eq!(safe_eps, 1e-10);
95    }
96
97    #[test]
98    fn test_safe_epsilon_explicit_precision() {
99        // epsilon=1e-7 > floor≈1.57e-16 → safe_eps should honour the user's epsilon
100        let (safe_eps, twork) = safe_epsilon(Some(1e-7), TworkType::Float64X2);
101        assert_eq!(twork, TworkType::Float64X2);
102        assert_eq!(safe_eps, 1e-7);
103    }
104
105    #[test]
106    #[should_panic(expected = "eps_required must be non-negative")]
107    fn test_negative_epsilon_panics() {
108        safe_epsilon(Some(-1.0), TworkType::Auto);
109    }
110
111    /// NaN is not an accuracy: it used to select the automatic accuracy.
112    /// compute_sve rejects it first; the public helper panics like it does
113    /// for a negative epsilon.
114    #[test]
115    #[should_panic(expected = "eps_required must be non-negative, got NaN")]
116    fn test_nan_epsilon_panics() {
117        safe_epsilon(Some(f64::NAN), TworkType::Auto);
118    }
119}
sparse_ir_basis/sve/types.rs

sparse_ir_basis/sve/
types.rs
