exponential_integral/
lib.rs

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
//! The exponential integral, often written $\text{Ei}$,
//! equal to the the integral of an exponentiated input over the input itself:
//! $\text{Ei}(t) = \int_{-\infty}^{t} \frac{ e^{u} }{ u } \text{d}u$
//!
//! Inspired by [GSL's implementation](https://github.com/ampl/gsl/blob/ff49e28bdffb893a1c0f6e3eff151296e0e71f82/specfunc/expint.c#L8).

#![no_std]
#![expect(non_snake_case, reason = "Proper mathematical names")]

pub mod chebyshev;
mod constants;
mod implementation;

#[cfg(test)]
mod test;

use {
    core::fmt,
    sigma_types::{Finite, Negative, NonZero, Positive},
};

#[cfg(feature = "error")]
use sigma_types::NonNegative;

/// An approximate value alongside an estimate of its own approximation error.
/// # Original C code
/// ```c
/// struct gsl_sf_result_struct {
///   double val;
///   double err;
/// };
/// typedef struct gsl_sf_result_struct gsl_sf_result;
/// ```
#[expect(clippy::exhaustive_structs, reason = "Simple structure")]
#[derive(Clone, Copy, Debug, PartialEq, PartialOrd)]
pub struct Approx {
    /// Estimate of the approximation error for `value`.
    #[cfg(feature = "error")]
    pub error: NonNegative<Finite<f64>>,
    /// Approximate value.
    pub value: Finite<f64>,
}

impl fmt::Display for Approx {
    #[inline]
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let Self {
            #[cfg(feature = "error")]
            ref error,
            ref value,
        } = *self;
        #[cfg(feature = "error")]
        {
            write!(f, "{value} +/- {error}")
        }
        #[cfg(not(feature = "error"))]
        {
            write!(f, "{value}")
        }
    }
}

/// An approximate value alongside an estimate of its own approximation error.
#[non_exhaustive]
#[derive(Clone, Copy, Debug, PartialEq, PartialOrd)]
pub enum Error {
    /// Argument was less than the safe minimum.
    ArgumentTooNegative(Negative<Finite<f64>>),
    /// Argument was less than the safe maximum.
    ArgumentTooPositive(Positive<Finite<f64>>),
}

impl fmt::Display for Error {
    #[inline]
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match *self {
            Self::ArgumentTooNegative(arg) => write!(
                f,
                "Argument too large (negative): minimum is {}, but {arg} was supplied",
                constants::NXMAX,
            ),
            Self::ArgumentTooPositive(arg) => write!(
                f,
                "Argument too large (positive): maximum is {}, but {arg} was supplied",
                constants::XMAX,
            ),
        }
    }
}

/// # Original C code
/// ```c
/// int gsl_sf_expint_E1_e(const double x, gsl_sf_result * result)
/// {
///   return expint_E1_impl(x, result, 0);
/// }
/// ```
///
/// # Errors
/// See `Error`.
#[inline]
pub fn E1(
    x: NonZero<Finite<f64>>,
    #[cfg(feature = "precision")] max_precision: usize,
) -> Result<Approx, Error> {
    implementation::E1(
        x,
        #[cfg(feature = "precision")]
        max_precision,
    )
}

/// # Original C code
/// ```c
/// int gsl_sf_expint_Ei_e(const double x, gsl_sf_result * result)
/// {
///   /* CHECK_POINTER(result) */
///
///   {
///     int status = gsl_sf_expint_E1_e(-x, result);
///     result->val = -result->val;
///     return status;
///   }
/// }
/// ```
///
/// # Errors
/// See `Error`.
#[inline(always)]
pub fn Ei(
    x: NonZero<Finite<f64>>,
    #[cfg(feature = "precision")] max_precision: usize,
) -> Result<Approx, Error> {
    #![expect(
        clippy::arithmetic_side_effects,
        reason = "property-based testing ensures this never happens"
    )]

    E1(
        -x,
        #[cfg(feature = "precision")]
        max_precision,
    )
    .map(|mut approx| {
        approx.value = -approx.value;
        approx
    })
}