exponential_integral/
chebyshev.rs

1//! Chebyshev series/polynomial approximation.
2
3use {
4    crate::Approx,
5    sigma_types::{Finite, Zero as _, usize::LessThan},
6};
7
8#[cfg(feature = "error")]
9use {crate::constants, sigma_types::NonNegative};
10
11/// Chebyshev series/polynomial approximation.
12/// # Original C code
13/// ```c
14/// struct cheb_series_struct {
15///   double * c;   /* coefficients                */
16///   int order;    /* order of expansion          */
17///   double a;     /* lower interval point        */
18///   double b;     /* upper interval point        */
19///   int order_sp; /* effective single precision order */
20/// };
21/// typedef struct cheb_series_struct cheb_series;
22///
23/// // ...
24///
25/// static inline int
26/// cheb_eval_e(const cheb_series * cs,
27///             const double x,
28///             gsl_sf_result * result)
29/// {
30///   int j;
31///   double d  = 0.0;
32///   double dd = 0.0;
33///
34///   double y  = (2.0*x - cs->a - cs->b) / (cs->b - cs->a);
35///   double two_x = 2.0 * y;
36///
37///   double e = 0.0;
38///
39///   for(j = cs->order; j>=1; j--) {
40///     double temp = d;
41///     d = two_x*d - dd + cs->c[j];
42///     e += fabs(two_x*temp) + fabs(dd) + fabs(cs->c[j]);
43///     dd = temp;
44///   }
45///
46///   {
47///     double temp = d;
48///     d = y*d - dd + 0.5 * cs->c[0];
49///     e += fabs(y*temp) + fabs(dd) + 0.5 * fabs(cs->c[0]);
50///   }
51///
52///   result->val = d;
53///   result->err = GSL_DBL_EPSILON * e + fabs(cs->c[cs->order]);
54///
55///   return GSL_SUCCESS;
56/// }
57/// ```
58#[inline]
59#[must_use]
60pub fn eval<const N_COEFFICIENTS: usize>(
61    coefficients: &[Finite<f64>; N_COEFFICIENTS],
62    x: Finite<f64>,
63    #[cfg(feature = "precision")] order: LessThan<{ N_COEFFICIENTS }>,
64) -> Approx {
65    #![expect(
66        clippy::arithmetic_side_effects,
67        reason = "property-based testing ensures this never happens"
68    )]
69
70    debug_assert!(
71        N_COEFFICIENTS > 0,
72        "Chebyshev series without any coefficients",
73    );
74
75    let two_x: Finite<f64> = Finite::new(2_f64) * x;
76
77    #[cfg(feature = "error")]
78    let mut e = NonNegative::<Finite<f64>>::ZERO;
79
80    let mut d = Finite::<f64>::ZERO;
81    let mut dd = Finite::<f64>::ZERO;
82
83    {
84        let mut j: LessThan<{ N_COEFFICIENTS }> = {
85            #[cfg(feature = "precision")]
86            {
87                order
88            }
89            #[cfg(not(feature = "precision"))]
90            {
91                LessThan::new(const { N_COEFFICIENTS - 1 })
92            }
93        };
94        while *j >= 1 {
95            // SAFETY:
96            // See the `debug_assert` above.
97            let coefficient = *unsafe { coefficients.get_unchecked(*j) };
98            let tmp = d;
99            d = ((two_x * d) - dd) + coefficient;
100            #[cfg(feature = "error")]
101            {
102                e += NonNegative::<Finite<f64>>::new((two_x * tmp).map(f64::abs))
103                    + NonNegative::<Finite<f64>>::new(dd.map(f64::abs))
104                    + NonNegative::<Finite<f64>>::new(coefficient.map(f64::abs));
105            }
106            dd = tmp;
107
108            j.map_mut(|u| *u -= 1);
109        }
110    }
111
112    {
113        #[cfg(feature = "error")]
114        let tmp = d;
115        // SAFETY:
116        // Sigma types ensure validity.
117        let coefficient = *unsafe { coefficients.get_unchecked(0) };
118        let half_coefficient = coefficient.map(|c| 0.5_f64 * c);
119        d = x * d - dd + half_coefficient;
120        #[cfg(feature = "error")]
121        {
122            e += NonNegative::<Finite<f64>>::new((x * tmp).map(f64::abs))
123                + NonNegative::<Finite<f64>>::new(dd.map(f64::abs))
124                + NonNegative::<Finite<f64>>::new(half_coefficient.map(f64::abs));
125        }
126    }
127
128    #[cfg(feature = "error")]
129    // SAFETY:
130    // See `debug_assert`s above.
131    let last_coefficient = *unsafe { coefficients.get_unchecked(const { N_COEFFICIENTS - 1 }) };
132
133    Approx {
134        value: d,
135        #[cfg(feature = "error")]
136        error: NonNegative::new(Finite::new(constants::GSL_DBL_EPSILON)) * e
137            + NonNegative::new(last_coefficient.map(f64::abs)),
138    }
139}
140
141/// Compile-time-compatible minimum of two large unsigned integers.
142#[inline]
143#[cfg_attr(not(test), expect(dead_code, reason = "TODO: REMOVE"))]
144#[cfg_attr(test, expect(clippy::single_call_fn, reason = "TODO: REMOVE"))]
145pub(crate) const fn min(a: usize, b: usize) -> usize {
146    if a.checked_sub(b).is_some() { b } else { a }
147}
exponential_integral/chebyshev.rs

exponential_integral/
chebyshev.rs
