peroxide/special/
function.rs

1use crate::special::lanczos::{gamma_approx, ln_gamma_approx};
2use std::f64::consts::PI;
3
4/// Gaussian function
5///
6/// `N(x|μ,σ) = 1/√(2πσ^2) exp(-(x-μ)^2/(2σ^2))`
7pub fn gaussian(x: f64, mu: f64, sigma: f64) -> f64 {
8    1f64 / ((2f64 * PI).sqrt() * sigma) * (-0.5 * ((x - mu) / sigma).powi(2)).exp()
9}
10
11/// Gamma function
12///
13/// # Description
14/// Use Lanczos approximation to implement Gamma function ($g=5, n=7$)
15///
16/// # References
17/// * [Robert Munafo, Coefficients for the Lanczos Approximation to the Gamma Function](https://mrob.com/pub/ries/lanczos-gamma.html)
18/// * [Paul Godfrey, A note on the computation of the convergent Lanczos complex Gamma approximation (web page), 2001.](http://my.fit.edu/~gabdo/gamma.txt)
19pub fn gamma(x: f64) -> f64 {
20    gamma_approx(x)
21}
22
23/// Logarithm Gamma function
24///
25/// # Description
26/// Use Lanczos approximation to implement Gamma function ($g=5, n=7$)
27///
28/// # References
29/// * [Robert Munafo, Coefficients for the Lanczos Approximation to the Gamma Function](https://mrob.com/pub/ries/lanczos-gamma.html)
30/// * [Paul Godfrey, A note on the computation of the convergent Lanczos complex Gamma approximation (web page), 2001.](http://my.fit.edu/~gabdo/gamma.txt)
31pub fn ln_gamma(x: f64) -> f64 {
32    ln_gamma_approx(x)
33}
34
35/// Pochhammer symbol
36pub fn poch(x: f64, n: usize) -> f64 {
37    let mut s = 1f64;
38    for i in 0..n {
39        s *= x + i as f64;
40    }
41    s
42}
43
44// /// Digamma function
45// ///
46// /// Wrapper of `digamma` function of `special` crate
47// pub fn digamma(x: f64) -> f64 {
48//     x.digamma()
49// }
50
51/// Regularized incomplete gamma integral (Lower)
52///
53/// Wrapper of `gammp` function of `puruspe` crate
54pub fn inc_gamma(a: f64, x: f64) -> f64 {
55    puruspe::gammp(a, x)
56}
57
58/// Inverse of regularized incomplete gamma integral (Lower)
59///
60/// Wrapper of `invgammp` function of `puruspe` crate
61pub fn inv_inc_gamma(p: f64, a: f64) -> f64 {
62    puruspe::invgammp(p, a)
63}
64
65/// Error function
66///
67/// Wrapper of `erf` function of `puruspe` crate
68pub fn erf(x: f64) -> f64 {
69    puruspe::erf(x)
70}
71
72/// Complement error function
73///
74/// Wrapper of `erfc` function of `puruspe` crate
75pub fn erfc(x: f64) -> f64 {
76    puruspe::erfc(x)
77}
78
79/// Inverse error function
80///
81/// Wrapper of `inverf` function of `puruspe` crate
82pub fn inv_erf(x: f64) -> f64 {
83    puruspe::inverf(x)
84}
85
86/// Inverse complementary error function
87///
88/// Wrapper of `inverfc` function of `puruspe` crate
89pub fn inv_erfc(p: f64) -> f64 {
90    puruspe::inverfc(p)
91}
92
93/// Beta function
94///
95/// Wrapper of `beta` function of `puruspe` crate
96pub fn beta(a: f64, b: f64) -> f64 {
97    puruspe::beta(a, b)
98}
99
100/// Regularized incomplete Beta function
101///
102/// Wrapper of `betai` function of `puruspe` crate
103pub fn inc_beta(a: f64, b: f64, x: f64) -> f64 {
104    puruspe::betai(a, b, x)
105}
106
107/// Inverse regularized incomplete beta function
108///
109/// Wrapper of `invbetai` function of `puruspe` crate
110pub fn inv_inc_beta(p: f64, a: f64, b: f64) -> f64 {
111    puruspe::invbetai(p, a, b)
112}
113
114/// Phi (CDF for Normal Dist)
115///
116/// $$\Phi(x) = \frac{1}{2}\left[1 + \text{erf}\left(\frac{x}{\sqrt{2}}\right) \right]$$
117pub fn phi(x: f64) -> f64 {
118    0.5 * (1f64 + erf(x / 2f64.sqrt()))
119}
120
121/// The principal branch of the Lambert W function, W_0(`z`).
122///
123/// Returns [`NAN`](f64::NAN) if the given input is smaller than -1/e (≈ -0.36787944117144233).
124///
125/// Use [`Precise`](LambertWAccuracyMode::Precise) for 50 bits of accuracy and the [`Simple`](LambertWAccuracyMode::Simple) mode
126/// for only 24 bits, but with faster execution time.
127///
128/// Wrapper of the `lambert_w_0` and `sp_lambert_w_0` functions of the `puruspe` crate.
129///
130/// # Reference
131///
132/// [Toshio Fukushima, Precise and fast computation of Lambert W function by piecewise minimax rational function approximation with variable transformation](https://www.researchgate.net/publication/346309410_Precise_and_fast_computation_of_Lambert_W_function_by_piecewise_minimax_rational_function_approximation_with_variable_transformation)
133pub fn lambert_w0(z: f64, mode: LambertWAccuracyMode) -> f64 {
134    match mode {
135        LambertWAccuracyMode::Precise => puruspe::lambert_w0(z),
136        LambertWAccuracyMode::Simple => puruspe::sp_lambert_w0(z),
137    }
138}
139
140/// The secondary branch of the Lambert W function, W_-1(`z`).
141///
142/// Returns [`NAN`](f64::NAN) if the given input is positive or smaller than -1/e (≈ -0.36787944117144233).
143///
144/// Use [`Precise`](LambertWAccuracyMode::Precise) for 50 bits of accuracy and the [`Simple`](LambertWAccuracyMode::Simple) mode
145/// for only 24 bits, but with faster execution time.
146///
147/// Wrapper of the `lambert_w_m1` and `sp_lambert_w_m1` functions of the `puruspe` crate.
148///
149/// # Reference
150///
151/// [Toshio Fukushima, Precise and fast computation of Lambert W function by piecewise minimax rational function approximation with variable transformation](https://www.researchgate.net/publication/346309410_Precise_and_fast_computation_of_Lambert_W_function_by_piecewise_minimax_rational_function_approximation_with_variable_transformation)
152pub fn lambert_wm1(z: f64, mode: LambertWAccuracyMode) -> f64 {
153    match mode {
154        LambertWAccuracyMode::Precise => puruspe::lambert_wm1(z),
155        LambertWAccuracyMode::Simple => puruspe::sp_lambert_wm1(z),
156    }
157}
158
159/// Decides the accuracy mode of the Lambert W functions.
160#[derive(Debug, Clone, Copy, PartialEq, Eq)]
161pub enum LambertWAccuracyMode {
162    /// Faster, 24 bits of accuracy.
163    Simple,
164    /// Slower, 50 bits of accuracy.
165    Precise,
166}
peroxide/special/function.rs

peroxide/special/
function.rs
