Special functions¶
Table of contents
Binomial function¶
-
template<typename T1, typename T2>
constexpr common_t<T1, T2> binomial_coef(const T1 n, const T2 k) noexcept¶ Compile-time binomial coefficient.
- Parameters
n – integral-valued input.
k – integral-valued input.
- Returns
computes the Binomial coefficient
\[ \binom{n}{k} = \frac{n!}{k!(n-k)!} \]also known as ‘nchoosek‘.
-
template<typename T1, typename T2>
constexpr common_return_t<T1, T2> log_binomial_coef(const T1 n, const T2 k) noexcept¶ Compile-time log binomial coefficient.
- Parameters
n – integral-valued input.
k – integral-valued input.
- Returns
computes the log Binomial coefficient
\[ \ln \frac{n!}{k!(n-k)!} = \ln \Gamma(n+1) - [ \ln \Gamma(k+1) + \ln \Gamma(n-k+1) ] \]
Beta function¶
-
template<typename T1, typename T2>
constexpr common_return_t<T1, T2> beta(const T1 a, const T2 b) noexcept¶ Compile-time beta function.
- Parameters
a – a real-valued input.
b – a real-valued input.
- Returns
the beta function using
\[ \text{B}(\alpha,\beta) := \int_0^1 t^{\alpha - 1} (1-t)^{\beta - 1} dt = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)} \]where \( \Gamma \) denotes the gamma function.
-
template<typename T1, typename T2>
constexpr common_return_t<T1, T2> lbeta(const T1 a, const T2 b) noexcept¶ Compile-time log-beta function.
- Parameters
a – a real-valued input.
b – a real-valued input.
- Returns
the log-beta function using
\[ \ln \text{B}(\alpha,\beta) := \ln \int_0^1 t^{\alpha - 1} (1-t)^{\beta - 1} dt = \ln \Gamma(\alpha) + \ln \Gamma(\beta) - \ln \Gamma(\alpha + \beta) \]where \( \Gamma \) denotes the gamma function.
Gamma function¶
-
template<typename T>
constexpr return_t<T> tgamma(const T x) noexcept¶ Compile-time gamma function.
- Parameters
x – a real-valued input.
- Returns
computes the true gamma function
\[ \Gamma(x) = \int_0^\infty y^{x-1} \exp(-y) dy \]using a polynomial form:\[ \Gamma(x+1) \approx (x+g+0.5)^{x+0.5} \exp(-x-g-0.5) \sqrt{2 \pi} \left[ c_0 + \frac{c_1}{x+1} + \frac{c_2}{x+2} + \cdots + \frac{c_n}{x+n} \right] \]where the value \( g \) and the coefficients \( (c_0, c_1, \ldots, c_n) \) are taken from Paul Godfrey, whose note can be found here: http://my.fit.edu/~gabdo/gamma.txt
-
template<typename T>
constexpr return_t<T> lgamma(const T x) noexcept¶ Compile-time log-gamma function.
- Parameters
x – a real-valued input.
- Returns
computes the log-gamma function
\[ \ln \Gamma(x) = \ln \int_0^\infty y^{x-1} \exp(-y) dy \]using a polynomial form:\[ \Gamma(x+1) \approx (x+g+0.5)^{x+0.5} \exp(-x-g-0.5) \sqrt{2 \pi} \left[ c_0 + \frac{c_1}{x+1} + \frac{c_2}{x+2} + \cdots + \frac{c_n}{x+n} \right] \]where the value \( g \) and the coefficients \( (c_0, c_1, \ldots, c_n) \) are taken from Paul Godfrey, whose note can be found here: http://my.fit.edu/~gabdo/gamma.txt
-
template<typename T1, typename T2>
constexpr return_t<T1> lmgamma(const T1 a, const T2 p) noexcept¶ Compile-time log multivariate gamma function.
- Parameters
a – a real-valued input.
p – integral-valued input.
- Returns
computes log-multivariate gamma function via recursion
\[ \Gamma_p(a) = \pi^{(p-1)/2} \Gamma(a) \Gamma_{p-1}(a-0.5) \]where \( \Gamma_1(a) = \Gamma(a) \).
Incomplete integral functions¶
-
template<typename T>
constexpr return_t<T> erf(const T x) noexcept¶ Compile-time Gaussian error function.
- Parameters
x – a real-valued input.
- Returns
computes the Gaussian error function
\[ \text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x \exp( - t^2) dt \]using a continued fraction representation:\[ \text{erf}(x) = \frac{2x}{\sqrt{\pi}} \exp(-x^2) \dfrac{1}{1 - 2x^2 + \dfrac{4x^2}{3 - 2x^2 + \dfrac{8x^2}{5 - 2x^2 + \dfrac{12x^2}{7 - 2x^2 + \ddots}}}} \]
-
template<typename T1, typename T2, typename T3>
constexpr common_return_t<T1, T2, T3> incomplete_beta(const T1 a, const T2 b, const T3 z) noexcept¶ Compile-time regularized incomplete beta function.
- Parameters
a – a real-valued, non-negative input.
b – a real-valued, non-negative input.
z – a real-valued, non-negative input.
- Returns
computes the regularized incomplete beta function,
\[ \frac{\text{B}(z;\alpha,\beta)}{\text{B}(\alpha,\beta)} = \frac{1}{\text{B}(\alpha,\beta)}\int_0^z t^{a-1} (1-t)^{\beta-1} dt \]using a continued fraction representation, found in the Handbook of Continued Fractions for Special Functions, and a modified Lentz method.\[ \frac{\text{B}(z;\alpha,\beta)}{\text{B}(\alpha,\beta)} = \frac{z^{\alpha} (1-t)^{\beta}}{\alpha \text{B}(\alpha,\beta)} \dfrac{a_1}{1 + \dfrac{a_2}{1 + \dfrac{a_3}{1 + \dfrac{a_4}{1 + \ddots}}}} \]where \( a_1 = 1 \) and\[ a_{2m+2} = - \frac{(\alpha + m)(\alpha + \beta + m)}{(\alpha + 2m)(\alpha + 2m + 1)}, \ m \geq 0 \]\[ a_{2m+1} = \frac{m(\beta - m)}{(\alpha + 2m - 1)(\alpha + 2m)}, \ m \geq 1 \]The Lentz method works as follows: let \( f_j \) denote the value of the continued fraction up to the first \( j \) terms; \( f_j \) is updated as follows:\[ c_j = 1 + a_j / c_{j-1}, \ \ d_j = 1 / (1 + a_j d_{j-1}) \]\[ f_j = c_j d_j f_{j-1} \]
-
template<typename T1, typename T2>
constexpr common_return_t<T1, T2> incomplete_gamma(const T1 a, const T2 x) noexcept¶ Compile-time regularized lower incomplete gamma function.
- Parameters
a – a real-valued, non-negative input.
x – a real-valued, non-negative input.
- Returns
the regularized lower incomplete gamma function evaluated at (
a,x),\[ \frac{\gamma(a,x)}{\Gamma(a)} = \frac{1}{\Gamma(a)} \int_0^x t^{a-1} \exp(-t) dt \]Whenais not too large, the value is computed using the continued fraction representation of the upper incomplete gamma function, \( \Gamma(a,x) \), using\[ \Gamma(a,x) = \Gamma(a) - \dfrac{x^a\exp(-x)}{a - \dfrac{ax}{a + 1 + \dfrac{x}{a + 2 - \dfrac{(a+1)x}{a + 3 + \dfrac{2x}{a + 4 - \ddots}}}}} \]where \( \gamma(a,x) \) and \( \Gamma(a,x) \) are connected via\[ \frac{\gamma(a,x)}{\Gamma(a)} + \frac{\Gamma(a,x)}{\Gamma(a)} = 1 \]When \( a > 10 \), a 50-point Gauss-Legendre quadrature scheme is employed.
Inverse incomplete integral functions¶
-
template<typename T>
constexpr return_t<T> erf_inv(const T p) noexcept¶ Compile-time inverse Gaussian error function.
- Parameters
p – a real-valued input with values in the unit-interval.
- Returns
Computes the inverse Gaussian error function, a value \( x \) such that
\[ f(x) := \text{erf}(x) - p \]is equal to zero, for a givenp. GCE-Math finds this root using Halley’s method:\[ x_{n+1} = x_n - \frac{f(x_n)/f'(x_n)}{1 - 0.5 \frac{f(x_n)}{f'(x_n)} \frac{f''(x_n)}{f'(x_n)} } \]where\[ \frac{\partial}{\partial x} \text{erf}(x) = \exp(-x^2), \ \ \frac{\partial^2}{\partial x^2} \text{erf}(x) = -2x\exp(-x^2) \]
-
template<typename T1, typename T2, typename T3>
constexpr common_t<T1, T2, T3> incomplete_beta_inv(const T1 a, const T2 b, const T3 p) noexcept¶ Compile-time inverse incomplete beta function.
- Parameters
a – a real-valued, non-negative input.
b – a real-valued, non-negative input.
p – a real-valued input with values in the unit-interval.
- Returns
Computes the inverse incomplete beta function, a value \( x \) such that
\[ f(x) := \frac{\text{B}(x;\alpha,\beta)}{\text{B}(\alpha,\beta)} - p \]equal to zero, for a givenp. GCE-Math finds this root using Halley’s method:\[ x_{n+1} = x_n - \frac{f(x_n)/f'(x_n)}{1 - 0.5 \frac{f(x_n)}{f'(x_n)} \frac{f''(x_n)}{f'(x_n)} } \]where\[ \frac{\partial}{\partial x} \left(\frac{\text{B}(x;\alpha,\beta)}{\text{B}(\alpha,\beta)}\right) = \frac{1}{\text{B}(\alpha,\beta)} x^{\alpha-1} (1-x)^{\beta-1} \]\[ \frac{\partial^2}{\partial x^2} \left(\frac{\text{B}(x;\alpha,\beta)}{\text{B}(\alpha,\beta)}\right) = \frac{1}{\text{B}(\alpha,\beta)} x^{\alpha-1} (1-x)^{\beta-1} \left( \frac{\alpha-1}{x} - \frac{\beta-1}{1 - x} \right) \]
-
template<typename T1, typename T2>
constexpr common_return_t<T1, T2> incomplete_gamma_inv(const T1 a, const T2 p) noexcept¶ Compile-time inverse incomplete gamma function.
- Parameters
a – a real-valued, non-negative input.
p – a real-valued input with values in the unit-interval.
- Returns
Computes the inverse incomplete gamma function, a value \( x \) such that
\[ f(x) := \frac{\gamma(a,x)}{\Gamma(a)} - p \]equal to zero, for a givenp. GCE-Math finds this root using Halley’s method:\[ x_{n+1} = x_n - \frac{f(x_n)/f'(x_n)}{1 - 0.5 \frac{f(x_n)}{f'(x_n)} \frac{f''(x_n)}{f'(x_n)} } \]where\[ \frac{\partial}{\partial x} \left(\frac{\gamma(a,x)}{\Gamma(a)}\right) = \frac{1}{\Gamma(a)} x^{a-1} \exp(-x) \]\[ \frac{\partial^2}{\partial x^2} \left(\frac{\gamma(a,x)}{\Gamma(a)}\right) = \frac{1}{\Gamma(a)} x^{a-1} \exp(-x) \left( \frac{a-1}{x} - 1 \right) \]