Merge pull request #17 from jorenham/gamma

Author: jorenham

src/ffi.rs

@@ -40,9 +40,19 @@ unsafe extern "C" {
 
     // boost/math/special_functions/gamma.hpp
     pub fn math_tgamma(x: f64) -> f64;
-    pub fn math_lgamma(x: f64) -> f64;
-    pub fn math_gamma_p(a: f64, x: f64) -> f64;
+    pub fn math_tgamma1pm1(x: f64) -> f64;
+    pub fn math_tgamma_(a: f64, x: f64) -> f64;
+    pub fn math_tgamma_lower(a: f64, x: f64) -> f64;
+    pub fn math_tgamma_ratio(a: f64, b: f64) -> f64;
+    pub fn math_tgamma_delta_ratio(x: f64, delta: f64) -> f64;
+    pub fn math_lgamma(x: f64, sign: *mut c_int) -> f64;
     pub fn math_gamma_q(a: f64, x: f64) -> f64;
+    pub fn math_gamma_q_inv(a: f64, q: f64) -> f64;
+    pub fn math_gamma_q_inva(x: f64, q: f64) -> f64;
+    pub fn math_gamma_p(a: f64, x: f64) -> f64;
+    pub fn math_gamma_p_inv(a: f64, p: f64) -> f64;
+    pub fn math_gamma_p_inva(x: f64, p: f64) -> f64;
+    pub fn math_gamma_p_derivative(a: f64, x: f64) -> f64;
 
     // boost/math/special_functions/jacobi.hpp
     pub fn math_jacobi(n: c_uint, alpha: f64, beta: f64, x: f64) -> f64;

src/math/special_functions/gamma.rs

@@ -1,75 +1,272 @@
 //! boost/math/special_functions/gamma.hpp
 //!
-//! # TODO:
-//! - `tgamma1pm1`
-//! - `tgamma_lower`
-//! - `tgamma_delta_ratio`
-//! - `tgamma_ratio`
-//! - `gamma_p_derivative`
-//! - `gamma_p_inv`
-//! - `gamma_p_inva`
-//! - `gamma_q_inv`
-//! - `gamma_q_inva`
+//! Note that we deviate from the original "tgamma" names by dropping the "t" prefix.
 
 use crate::ffi;
+use core::ffi::c_int;
 
 /// Gamma function *Γ(x)*
 ///
+/// Corresponds to `boost::math::tgamma(x)` in C++.
 /// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_gamma/tgamma.html>
-pub fn tgamma(x: f64) -> f64 {
+pub fn gamma(x: f64) -> f64 {
     unsafe { ffi::math_tgamma(x) }
 }
 
-/// Log-Gamma function *ln |Γ(x)|*
+/// Accurate evaluation of `tgamma(x + 1) - 1` for very small `x`
 ///
-/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_gamma/lgamma.html>
-pub fn lgamma(x: f64) -> f64 {
-    unsafe { ffi::math_lgamma(x) }
+/// Internally the implementation does not make use of the addition and subtraction implied by the
+/// definition, leading to accurate results even for very small `x`.
+///
+/// See [`gamma`] for the gamma function itself.
+///
+/// Corresponds to `boost::math::tgamma1pm1(x)` in C++.
+/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_gamma/tgamma.html>
+pub fn gamma1pm1(x: f64) -> f64 {
+    unsafe { ffi::math_tgamma1pm1(x) }
 }
 
-/// Incomplete gamma function *P(a,x) = γ(a,x) / Γ(a)*
+/// Upper incomplete gamma function *Γ(a,x)*
+///
+/// See also:
+/// - [`gamma`]: Gamma function *Γ(x)*
+/// - [`gamma_lower`]: Lower incomplete gamma function *γ(a,x)*
+/// - [`gamma_q`]: Normalized upper incomplete gamma function *Q(a,x) = Γ(a,x) / Γ(a)*
 ///
+/// Corresponds to `boost::math::tgamma(a, x)` in C++.
 /// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html>
-pub fn gamma_p(a: f64, x: f64) -> f64 {
-    unsafe { ffi::math_gamma_p(a, x) }
+pub fn gamma_upper(a: f64, x: f64) -> f64 {
+    unsafe { ffi::math_tgamma_(a, x) }
 }
 
-/// Incomplete gamma function *Q(a,x) = Γ(a,x) / Γ(a)*
+/// Lower incomplete gamma function *γ(a,x)*
+///
+/// See also:
+/// - [`gamma`]: Gamma function *Γ(x)*
+/// - [`gamma_upper`]: Upper incomplete gamma function *Γ(a,x)*
+/// - [`gamma_p`]: Normalized lower incomplete gamma function *P(a,x) = γ(a,x) / Γ(a)*
 ///
+/// Corresponds to `boost::math::tgamma_lower(a, x)` in C++.
+/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html>
+pub fn gamma_lower(a: f64, x: f64) -> f64 {
+    unsafe { ffi::math_tgamma_lower(a, x) }
+}
+
+/// Ratio of two gamma functions *Γ(a) / Γ(b)*
+///
+/// See [`gamma`] for the gamma function itself.
+///
+/// Corresponds to `boost::math::tgamma_ratio(a, b)` in C++.
+/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_gamma/gamma_ratios.html>
+pub fn gamma_ratio(a: f64, b: f64) -> f64 {
+    unsafe { ffi::math_tgamma_ratio(a, b) }
+}
+
+/// Ratio of two gamma functions *Γ(x) / Γ(x + δ)*
+///
+/// See [`gamma`] for the gamma function itself.
+///
+/// Corresponds to `boost::math::tgamma_delta_ratio(x, delta)` in C++.
+/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_gamma/gamma_ratios.html>
+pub fn gamma_delta_ratio(x: f64, delta: f64) -> f64 {
+    unsafe { ffi::math_tgamma_delta_ratio(x, delta) }
+}
+
+/// Natural logarithm of the absolute value of the gamma function
+///
+/// The integer part of the tuple indicates the sign of the gamma function.
+///
+/// See [`gamma`] for the gamma function itself.
+///
+/// Corresponds to `boost::math::lgamma(x, *sign)` in C++.
+/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_gamma/lgamma.html>
+pub fn lgamma(x: f64) -> (f64, i32) {
+    let mut sign: c_int = 0;
+    let out = unsafe { ffi::math_lgamma(x, &mut sign) };
+    assert_ne!(sign, 0);
+    (out, sign)
+}
+
+/// Normalized upper incomplete gamma function *Q(a,x)*
+///
+/// *Q(a,x) = Γ(a,x) / Γ(a)*
+///
+/// See also:
+/// - [`gamma`]: Gamma function *Γ(x)*
+/// - [`gamma_upper`]: Upper incomplete gamma function *Γ(a,x)*
+/// - [`gamma_p`]: Normalized lower incomplete gamma function *P(a,x)*
+///
+/// Corresponds to `boost::math::gamma_q(a, x)` in C++.
 /// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html>
 pub fn gamma_q(a: f64, x: f64) -> f64 {
     unsafe { ffi::math_gamma_q(a, x) }
 }
 
+/// Inverse of [`gamma_q`] w.r.t. `x`
+///
+/// Corresponds to `boost::math::gamma_q_inv(a, p)` in C++.
+/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_gamma/igamma_inv.html>
+pub fn gamma_q_inv(a: f64, q: f64) -> f64 {
+    unsafe { ffi::math_gamma_q_inv(a, q) }
+}
+
+/// Inverse of [`gamma_q`] w.r.t. `a`
+///
+/// Corresponds to `boost::math::gamma_q_inva(x, p)` in C++.
+/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_gamma/igamma_inv.html>
+pub fn gamma_q_inva(x: f64, q: f64) -> f64 {
+    unsafe { ffi::math_gamma_q_inva(x, q) }
+}
+
+/// Normalized lower incomplete gamma function *P(a,x)*
+///
+/// *P(a,x) = γ(a,x) / Γ(a)*
+///
+/// See also:
+/// - [`gamma`]: Gamma function *Γ(x)*
+/// - [`gamma_lower`]: Lower incomplete gamma function *γ(a,x)*
+/// - [`gamma_q`]: Normalized upper incomplete gamma function *Q(a,x)*
+///
+/// Corresponds to `boost::math::gamma_p(a, x)` in C++.
+/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html>
+pub fn gamma_p(a: f64, x: f64) -> f64 {
+    unsafe { ffi::math_gamma_p(a, x) }
+}
+
+/// Inverse of [`gamma_p`] w.r.t. `x`
+///
+/// Corresponds to `boost::math::gamma_p_inv(a, p)` in C++.
+/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_gamma/igamma_inv.html>
+pub fn gamma_p_inv(a: f64, p: f64) -> f64 {
+    unsafe { ffi::math_gamma_p_inv(a, p) }
+}
+
+/// Inverse of [`gamma_p`] w.r.t. `a`
+///
+/// Corresponds to `boost::math::gamma_p_inva(x, p)` in C++.
+/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_gamma/igamma_inv.html>
+pub fn gamma_p_inva(x: f64, p: f64) -> f64 {
+    unsafe { ffi::math_gamma_p_inva(x, p) }
+}
+
+/// Derivative of the normalized lower incomplete gamma function
+///
+/// *P'(a,x) = e<sup>-x</sup> x<sup>a-1</sup> / Γ(a)*
+///
+/// Note that the derivative of the function [`gamma_q`] can be obtained by negating the result of
+/// this function, i.e., *Q'(a,x) = -P'(a,x)*.
+///
+/// Also note that this is the probability distribution function of the standard gamma distribution.
+///
+/// See [`gamma_p`] for the normalized lower incomplete gamma function.
+///
+/// Corresponds to `boost::math::gamma_p_derivative(a, x)` in C++.
+/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_gamma/gamma_derivatives.html>
+pub fn gamma_p_derivative(a: f64, x: f64) -> f64 {
+    unsafe { ffi::math_gamma_p_derivative(a, x) }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
 
-    const ATOL: f64 = f64::EPSILON;
+    const RTOL: f64 = f64::EPSILON;
+    const SQRT_PI: f64 = 1.772_453_850_905_516; // √π
+
+    #[test]
+    fn test_gamma() {
+        assert!(gamma(0.0).is_infinite());
+        assert_relative_eq!(gamma(-0.5), -2.0 * SQRT_PI, epsilon = RTOL);
+        assert_relative_eq!(gamma(0.5), SQRT_PI, epsilon = RTOL);
+        assert_relative_eq!(gamma(1.0), 1.0, epsilon = RTOL);
+        assert_relative_eq!(gamma(2.0), 1.0, epsilon = RTOL);
+        assert_relative_eq!(gamma(3.0), 2.0, epsilon = RTOL);
+    }
 
     #[test]
-    fn test_tgamma() {
-        assert_abs_diff_eq!(tgamma(1.0), 1.0, epsilon = ATOL);
-        assert_abs_diff_eq!(tgamma(0.5), core::f64::consts::PI.sqrt(), epsilon = ATOL);
-        assert_abs_diff_eq!(tgamma(5.0), 24.0, epsilon = ATOL);
+    fn test_gamma1pm1() {
+        assert_relative_eq!(gamma1pm1(-0.5), SQRT_PI - 1.0, epsilon = RTOL);
+        assert_relative_eq!(gamma1pm1(0.0), 0.0, epsilon = RTOL);
+        assert_relative_eq!(gamma1pm1(1.0), 0.0, epsilon = RTOL);
+        assert_relative_eq!(gamma1pm1(2.0), 1.0, epsilon = RTOL);
+        // results from Wolfram Alpha
+        assert_abs_diff_eq!(
+            gamma1pm1(-1e-14),
+            5.772_156_649_015_427_5e-15,
+            epsilon = 1e-30
+        );
+        assert_abs_diff_eq!(
+            gamma1pm1(1e-14),
+            -5.772_156_649_015_229_5e-15,
+            epsilon = 1e-30
+        );
+    }
+
+    #[test]
+    fn test_gamma_upper() {
+        assert!(gamma_upper(4.2, 0.5).is_finite());
+    }
+
+    #[test]
+    fn test_gamma_lower() {
+        assert!(gamma_lower(4.2, 0.5).is_finite());
+    }
+
+    #[test]
+    fn test_gamma_ratio() {
+        assert_relative_eq!(gamma_ratio(4.0, 6.0), 0.05, epsilon = RTOL);
+    }
+
+    #[test]
+    fn test_gamma_delta_ratio() {
+        assert_relative_eq!(gamma_delta_ratio(4.0, 2.0), 0.05, epsilon = RTOL);
     }
 
     #[test]
     fn test_lgamma() {
-        assert_abs_diff_eq!(lgamma(1.0), 0.0, epsilon = ATOL);
+        let (val, sign) = lgamma(0.5);
+        assert!(val.is_finite());
+        assert!(val > 0.0);
+        assert_eq!(sign, 1);
+
+        let (val, sign) = lgamma(-0.5);
+        assert!(val.is_finite());
+        assert!(val > 0.0);
+        assert_eq!(sign, -1);
+    }
+
+    #[test]
+    fn test_gamma_q() {
+        assert!(gamma_q(4.2, 0.5).is_finite());
+    }
+
+    #[test]
+    fn test_gamma_q_inv() {
+        assert!(gamma_q_inv(4.2, 0.5).is_finite());
+    }
+
+    #[test]
+    fn test_gamma_q_inva() {
+        assert!(gamma_q_inva(4.2, 0.5).is_finite());
     }
 
     #[test]
     fn test_gamma_p() {
-        assert_abs_diff_eq!(gamma_p(2.0, 0.0), 0.0, epsilon = ATOL);
+        assert!(gamma_p(4.2, 0.5).is_finite());
     }
+
     #[test]
-    fn test_gamma_q() {
-        assert_abs_diff_eq!(gamma_q(2.0, 0.0), 1.0, epsilon = ATOL);
+    fn test_gamma_p_inv() {
+        assert!(gamma_p_inv(4.2, 0.5).is_finite());
+    }
+
+    #[test]
+    fn test_gamma_p_inva() {
+        assert!(gamma_p_inva(4.2, 0.5).is_finite());
     }
 
     #[test]
-    fn test_gamma_p_q_relation() {
-        assert_abs_diff_eq!(gamma_p(2.0, 1.0) + gamma_q(2.0, 1.0), 1.0, epsilon = ATOL);
+    fn test_gamma_p_derivative() {
+        assert!(gamma_p_derivative(4.2, 0.5).is_finite());
     }
 }

wrapper.cpp

@@ -1,8 +1,7 @@
-// Simple C++ wrappers for the Boost Math functions
+// Simple C++ wrappers for the Boost Math functions for https://github.com/jorenham/boost-rust
 
-// Boost static config, see
-// https://www.boost.org/doc/libs/latest/boost/math/tools/user.hpp
 #ifndef BOOST_MATH_TOOLS_USER_HPP
+    // https://www.boost.org/doc/libs/latest/boost/math/tools/user.hpp
     #define BOOST_MATH_TOOLS_USER_HPP
     #define BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
     #define BOOST_MATH_DOMAIN_ERROR_POLICY errno_on_error
@@ -65,9 +64,19 @@ double math_rising_factorial(double x, int n) { return rising_factorial(x, n); }
 
 // boost/math/special_functions/gamma.hpp
 double math_tgamma(double x) { return tgamma(x); }
-double math_lgamma(double x) { return lgamma(x); }
-double math_gamma_p(double a, double x) { return gamma_p(a, x); }
+double math_tgamma1pm1(double x) { return tgamma1pm1(x); }
+double math_tgamma_(double a, double x) { return tgamma(a, x); }
+double math_tgamma_lower(double a, double x) { return tgamma_lower(a, x); }
+double math_tgamma_ratio(double a, double b) { return tgamma_ratio(a, b); }
+double math_tgamma_delta_ratio(double x, double delta) { return tgamma_delta_ratio(x, delta); }
+double math_lgamma(double x, int* sign) { return lgamma(x, sign); }
 double math_gamma_q(double a, double x) { return gamma_q(a, x); }
+double math_gamma_q_inv(double a, double q) { return gamma_q_inv(a, q); }
+double math_gamma_q_inva(double x, double q) { return gamma_q_inva(x, q); }
+double math_gamma_p(double a, double x) { return gamma_p(a, x); }
+double math_gamma_p_inv(double a, double p) { return gamma_p_inv(a, p); }
+double math_gamma_p_inva(double x, double p) { return gamma_p_inva(x, p); }
+double math_gamma_p_derivative(double a, double x) { return gamma_p_derivative(a, x); }
 
 // boost/math/special_functions/jacobi.hpp
 double math_jacobi(unsigned n, double alpha, double beta, double x) {