exponential_integral/
lib.rs

1//! The exponential integral, often written $\text{Ei}$,
2//! equal to the the integral of an exponentiated input over the input itself:
3//! $\text{Ei}(t) = \int_{-\infty}^{t} \frac{ e^{u} }{ u } \text{d}u$
4//!
5//! Inspired by [GSL's implementation](https://github.com/ampl/gsl/blob/ff49e28bdffb893a1c0f6e3eff151296e0e71f82/specfunc/expint.c#L8).
6
7#![no_std]
8#![expect(non_snake_case, reason = "Proper mathematical names")]
9
10pub mod chebyshev;
11mod constants;
12mod implementation;
13
14pub mod neg {
15    //! Inputs less than 0.
16
17    use {
18        crate::{Approx, constants, implementation::neg, pos},
19        core::fmt,
20        sigma_types::{Finite, Negative},
21    };
22
23    /// Argument too large (negative): minimum is `constants::NXMAX`, just under -710.
24    #[non_exhaustive]
25    #[derive(Clone, Copy, Debug, PartialEq, PartialOrd)]
26    pub struct HugeArgument(pub Negative<Finite<f64>>);
27
28    impl fmt::Display for HugeArgument {
29        #[inline]
30        fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
31            let Self(ref arg) = *self;
32            write!(
33                f,
34                "Argument too large (negative): minimum is {}, but {arg} was supplied",
35                constants::NXMAX,
36            )
37        }
38    }
39
40    /// E1 on inputs less than 0.
41    /// # Errors
42    /// If `x` is so large that floating-point operations will fail down the line (absolute value of just over 710).
43    #[inline]
44    pub fn E1(
45        x: Negative<Finite<f64>>,
46        #[cfg(feature = "precision")] max_precision: usize,
47    ) -> Result<Approx, HugeArgument> {
48        neg::E1(
49            x,
50            #[cfg(feature = "precision")]
51            max_precision,
52        )
53    }
54
55    /// Ei on inputs less than 0.
56    /// # Errors
57    /// If `x` is so large that floating-point operations will fail down the line (absolute value of just over 710).
58    #[inline(always)]
59    pub fn Ei(
60        x: Negative<Finite<f64>>,
61        #[cfg(feature = "precision")] max_precision: usize,
62    ) -> Result<Approx, HugeArgument> {
63        #![expect(
64            clippy::arithmetic_side_effects,
65            reason = "property-based testing ensures this never happens"
66        )]
67
68        pos::E1(
69            -x,
70            #[cfg(feature = "precision")]
71            max_precision,
72        )
73        .map(|mut approx| {
74            approx.value = -approx.value;
75            approx
76        })
77        .map_err(|pos::HugeArgument(arg)| HugeArgument(-arg))
78    }
79}
80
81pub mod pos {
82    //! Inputs greater than 0.
83
84    use {
85        crate::{Approx, constants, implementation::pos, neg},
86        core::fmt,
87        sigma_types::{Finite, Positive},
88    };
89
90    /// Argument too large (positive): maximum is `constants::XMAX`, just over 710.
91    #[non_exhaustive]
92    #[derive(Clone, Copy, Debug, PartialEq, PartialOrd)]
93    pub struct HugeArgument(pub Positive<Finite<f64>>);
94
95    impl fmt::Display for HugeArgument {
96        #[inline]
97        fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
98            let Self(ref arg) = *self;
99            write!(
100                f,
101                "Argument too large (positive): maximum is {}, but {arg} was supplied",
102                constants::XMAX,
103            )
104        }
105    }
106
107    /// E1 on inputs less than 0.
108    /// # Errors
109    /// If `x` is so large that floating-point operations will fail down the line (absolute value of just over 710).
110    #[inline]
111    pub fn E1(
112        x: Positive<Finite<f64>>,
113        #[cfg(feature = "precision")] max_precision: usize,
114    ) -> Result<Approx, HugeArgument> {
115        pos::E1(
116            x,
117            #[cfg(feature = "precision")]
118            max_precision,
119        )
120    }
121
122    /// Ei on inputs less than 0.
123    /// # Errors
124    /// If `x` is so large that floating-point operations will fail down the line (absolute value of just over 710).
125    #[inline(always)]
126    pub fn Ei(
127        x: Positive<Finite<f64>>,
128        #[cfg(feature = "precision")] max_precision: usize,
129    ) -> Result<Approx, HugeArgument> {
130        #![expect(
131            clippy::arithmetic_side_effects,
132            reason = "property-based testing ensures this never happens"
133        )]
134
135        neg::E1(
136            -x,
137            #[cfg(feature = "precision")]
138            max_precision,
139        )
140        .map(|mut approx| {
141            approx.value = -approx.value;
142            approx
143        })
144        .map_err(|neg::HugeArgument(arg)| HugeArgument(-arg))
145    }
146}
147
148#[cfg(test)]
149mod test;
150
151use {
152    core::fmt,
153    sigma_types::{Finite, Negative, NonZero, Positive},
154};
155
156#[cfg(feature = "error")]
157use sigma_types::NonNegative;
158
159/// An approximate value alongside an estimate of its own approximation error.
160/// # Original C code
161/// ```c
162/// struct gsl_sf_result_struct {
163///   double val;
164///   double err;
165/// };
166/// typedef struct gsl_sf_result_struct gsl_sf_result;
167/// ```
168#[expect(clippy::exhaustive_structs, reason = "Simple structure")]
169#[derive(Clone, Copy, Debug, PartialEq, PartialOrd)]
170pub struct Approx {
171    /// Estimate of the approximation error for `value`.
172    #[cfg(feature = "error")]
173    pub error: NonNegative<Finite<f64>>,
174    /// Approximate value.
175    pub value: Finite<f64>,
176}
177
178impl fmt::Display for Approx {
179    #[inline]
180    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
181        let Self {
182            #[cfg(feature = "error")]
183            ref error,
184            ref value,
185        } = *self;
186        #[cfg(feature = "error")]
187        {
188            write!(f, "{value} +/- {error}")
189        }
190        #[cfg(not(feature = "error"))]
191        {
192            write!(f, "{value}")
193        }
194    }
195}
196
197/// An approximate value alongside an estimate of its own approximation error.
198#[non_exhaustive]
199#[derive(Clone, Copy, Debug, PartialEq, PartialOrd)]
200pub enum Error {
201    /// Argument was less than the safe minimum.
202    ArgumentTooNegative(Negative<Finite<f64>>),
203    /// Argument was less than the safe maximum.
204    ArgumentTooPositive(Positive<Finite<f64>>),
205}
206
207impl fmt::Display for Error {
208    #[inline]
209    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
210        match *self {
211            Self::ArgumentTooNegative(arg) => fmt::Display::fmt(&neg::HugeArgument(arg), f),
212            Self::ArgumentTooPositive(arg) => fmt::Display::fmt(&pos::HugeArgument(arg), f),
213        }
214    }
215}
216
217/// # Original C code
218/// ```c
219/// int gsl_sf_expint_E1_e(const double x, gsl_sf_result * result)
220/// {
221///   return expint_E1_impl(x, result, 0);
222/// }
223/// ```
224///
225/// # Errors
226/// If `x` is so large that floating-point operations will fail down the line (absolute value of just over 710).
227#[inline]
228pub fn E1(
229    x: NonZero<Finite<f64>>,
230    #[cfg(feature = "precision")] max_precision: usize,
231) -> Result<Approx, Error> {
232    implementation::E1(
233        x,
234        #[cfg(feature = "precision")]
235        max_precision,
236    )
237}
238
239/// # Original C code
240/// ```c
241/// int gsl_sf_expint_Ei_e(const double x, gsl_sf_result * result)
242/// {
243///   /* CHECK_POINTER(result) */
244///
245///   {
246///     int status = gsl_sf_expint_E1_e(-x, result);
247///     result->val = -result->val;
248///     return status;
249///   }
250/// }
251/// ```
252///
253/// # Errors
254/// If `x` is so large that floating-point operations will fail down the line (absolute value of just over 710).
255#[inline(always)]
256pub fn Ei(
257    x: NonZero<Finite<f64>>,
258    #[cfg(feature = "precision")] max_precision: usize,
259) -> Result<Approx, Error> {
260    #![expect(
261        clippy::arithmetic_side_effects,
262        reason = "property-based testing ensures this never happens"
263    )]
264
265    E1(
266        -x,
267        #[cfg(feature = "precision")]
268        max_precision,
269    )
270    .map(|mut approx| {
271        approx.value = -approx.value;
272        approx
273    })
274}
exponential_integral/lib.rs

exponential_integral/
lib.rs
