Merge pull request #60 from jorenham/orthopoly-api-consistency

Author: jorenham

src/ffi.rs

@@ -137,6 +137,7 @@ unsafe extern "C" {
     pub(crate) fn math_gegenbauer_derivative(n: c_uint, lambda: f64, x: f64, k: c_uint) -> f64;
 
     // boost/math/special_functions/hermite.hpp
+    #[cfg(test)]
     pub(crate) fn math_hermite(n: c_uint, x: f64) -> f64;
 
     // boost/math/special_functions/heuman_lambda.hpp

src/math/mod.rs

@@ -97,8 +97,12 @@
 //!   - [`laguerre_assoc`]
 //!   - [`laguerre_assoc_next`]
 //! - Hermite Polynomials
-//!   - [`hermite`]
-//!   - [`hermite_next`]
+//!   - [`hermite_h`]
+//!   - [`hermite_h_derivative`]
+//!   - [`hermite_h_next`]
+//!   - [`hermite_he`]
+//!   - [`hermite_he_derivative`]
+//!   - [`hermite_he_next`]
 //! - Gegenbauer Polynomials
 //!   - [`gegenbauer`]
 //!   - [`gegenbauer_derivative`]

src/math/special_functions/chebyshev.rs

@@ -2,6 +2,39 @@
 
 use crate::ffi;
 
+/// Chebyshev polynomial of the 1st kind T<sub>n</sub>(x).
+///
+/// Defined as <i>T<sub>n</sub>(</i>cos<i>(θ)) = </i>cos<i>(n θ)</i>.
+///
+/// See [`chebyshev_t_prime`] for the derivative, and [`chebyshev_u`] for the second kind.
+///
+/// Corresponds to `boost::math::chebyshev_t` in C++.
+/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_poly/chebyshev.html>
+pub fn chebyshev_t(n: u32, x: f64) -> f64 {
+    unsafe { ffi::math_chebyshev_t(n, x) }
+}
+
+/// Derivative of [`chebyshev_t`].
+///
+/// Corresponds to `boost::math::chebyshev_t_prime` in C++.
+/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_poly/chebyshev.html>
+#[doc(alias = "chebyshev_t_derivative")]
+pub fn chebyshev_t_prime(n: u32, x: f64) -> f64 {
+    unsafe { ffi::math_chebyshev_t_prime(n, x) }
+}
+
+/// Chebyshev polynomial of the 2nd kind U<sub>n</sub>(x).
+///
+/// Defined as <i>U<sub>n</sub>(</i>cos<i>(θ)) = </i>sin<i>((n+1) θ) / </i>sin<i>(θ)</i>.
+///
+/// See [`chebyshev_t`] for the first kind.
+///
+/// Corresponds to `boost::math::chebyshev_u` in C++.
+/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_poly/chebyshev.html>
+pub fn chebyshev_u(n: u32, x: f64) -> f64 {
+    unsafe { ffi::math_chebyshev_u(n, x) }
+}
+
 /// Recurrence relation for Chebyshev polynomials
 ///
 /// *T<sub>n+1</sub>(x) = 2 x T<sub>n</sub>(x) - T<sub>n-1</sub>(x)*
@@ -20,8 +53,8 @@ use crate::ffi;
 /// let t1 = chebyshev_t(1, x); // x
 /// let t2 = chebyshev_t(2, x); // 2x² - 1
 /// let t3 = chebyshev_t(3, x); // 4x³ - 3x
-/// assert_relative_eq!(chebyshev_next(&x, &t1, &t0), t2);
-/// assert_relative_eq!(chebyshev_next(&x, &t2, &t1), t3);
+/// assert_relative_eq!(chebyshev_next(x, t1, t0), t2);
+/// assert_relative_eq!(chebyshev_next(x, t2, t1), t3);
 /// ```
 ///
 /// [`chebyshev_u`] recurrence:
@@ -34,46 +67,15 @@ use crate::ffi;
 /// let u1 = chebyshev_u(1, x); // 2x
 /// let u2 = chebyshev_u(2, x); // 4x² - 1
 /// let u3 = chebyshev_u(3, x); // 8x³ - 4x
-/// assert_relative_eq!(chebyshev_next(&x, &u1, &u0), u2);
-/// assert_relative_eq!(chebyshev_next(&x, &u2, &u1), u3);
+/// assert_relative_eq!(chebyshev_next(x, u1, u0), u2);
+/// assert_relative_eq!(chebyshev_next(x, u2, u1), u3);
 /// ```
-#[allow(non_snake_case)]
 #[inline(always)]
-pub fn chebyshev_next(x: &f64, Tn: &f64, Tn_1: &f64) -> f64 {
-    2.0 * x * Tn - Tn_1
-}
-
-/// Chebyshev polynomial of the 1st kind T<sub>n</sub>(x).
-///
-/// Defined as <i>T<sub>n</sub>(</i>cos<i>(θ)) = </i>cos<i>(n θ)</i>.
-///
-/// See [`chebyshev_t_prime`] for the derivative, and [`chebyshev_u`] for the second kind.
-///
-/// Corresponds to `boost::math::chebyshev_t` in C++.
-/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_poly/chebyshev.html>
-pub fn chebyshev_t(n: u32, x: f64) -> f64 {
-    unsafe { ffi::math_chebyshev_t(n, x) }
-}
-
-/// Derivative of [`chebyshev_t`].
-///
-/// Corresponds to `boost::math::chebyshev_t_prime` in C++.
-/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_poly/chebyshev.html>
-#[doc(alias = "chebyshev_t_derivative")]
-pub fn chebyshev_t_prime(n: u32, x: f64) -> f64 {
-    unsafe { ffi::math_chebyshev_t_prime(n, x) }
-}
-
-/// Chebyshev polynomial of the 2nd kind U<sub>n</sub>(x).
-///
-/// Defined as <i>U<sub>n</sub>(</i>cos<i>(θ)) = </i>sin<i>((n+1) θ) / </i>sin<i>(θ)</i>.
-///
-/// See [`chebyshev_t`] for the first kind.
-///
-/// Corresponds to `boost::math::chebyshev_u` in C++.
-/// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_poly/chebyshev.html>
-pub fn chebyshev_u(n: u32, x: f64) -> f64 {
-    unsafe { ffi::math_chebyshev_u(n, x) }
+#[allow(non_snake_case)]
+#[doc(alias = "chebyshev_t_next")]
+#[doc(alias = "chebyshev_u_next")]
+pub fn chebyshev_next(x: f64, Tn: f64, Tn_prev: f64) -> f64 {
+    2.0 * x * Tn - Tn_prev
 }
 
 #[cfg(test)]

src/math/special_functions/hermite.rs

@@ -1,48 +1,202 @@
 //! boost/math/special_functions/hermite.hpp
 
-use crate::ffi;
-use core::ffi::c_uint;
+/// Recurrence relation for [`hermite_h`]
+///
+/// *H<sub>n+1</sub>(x) = 2xH<sub>n</sub>(x) - 2nH<sub>n-1</sub>(x)*
+///
+/// # Examples
+///
+/// ```
+/// # use boost::math::{hermite_h, hermite_h_next};
+/// let x = 0.42;
+/// let p0 = hermite_h(0, x); // 1
+/// let p1 = hermite_h(1, x); // 2x
+/// let p2 = hermite_h(2, x); // 4x² - 2
+/// let p3 = hermite_h(3, x); // 8x³ - 12x
+/// assert_eq!(hermite_h_next(1, x, p1, p0), p2);
+/// assert_eq!(hermite_h_next(2, x, p2, p1), p3);
+/// ```
+#[inline(always)]
+pub fn hermite_h_next(n: u32, x: f64, pn: f64, pn_prev: f64) -> f64 {
+    2.0 * hermite_he_next(n, x, pn, pn_prev)
+}
 
 /// Hermite Polynomial *H<sub>n</sub>(x)*
 ///
-/// Corresponds to `boost::math::hermite(n, x)`.
+/// Note that this is the  "physicist's" Hermite polynomial.
+///
+/// This is a pure rust implementation equivalent to the `boost::math::hermite(n, x)` C++ function.
 /// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_poly/hermite.html>
-pub fn hermite(n: u32, x: f64) -> f64 {
-    unsafe { ffi::math_hermite(n as c_uint, x) }
+pub fn hermite_h(n: u32, x: f64) -> f64 {
+    // Implement Hermite polynomials via recurrence:
+    let (mut p0, mut p1) = (1.0, 2.0 * x);
+    if n == 0 {
+        p0
+    } else {
+        for c in 1..n {
+            (p0, p1) = (p1, hermite_h_next(c, x, p1, p0));
+        }
+        p1
+    }
 }
 
-/// Recurrence relation for [`hermite`]
+/// *k*-th derivative of the Hermite polynomial *H<sub>n</sub>(x)*
 ///
-/// *H<sub>n+1</sub>(x) = 2xH<sub>n</sub>(x) - 2nH<sub>n-1</sub>(x)*
+/// This function does not exist in the Boost Math C++ library.
+#[doc(alias = "hermite_h_prime")]
+pub fn hermite_h_derivative(n: u32, x: f64, k: u32) -> f64 {
+    if n < k {
+        0.0
+    } else {
+        let mut p = hermite_h(n - k, x);
+        for m in 0..k {
+            p *= 2.0 * (n - m) as f64
+        }
+        p
+    }
+}
+
+/// Recurrence relation for [`hermite_he`]
+///
+/// *He<sub>n+1</sub>(x) = x He<sub>n</sub>(x) - n He<sub>n-1</sub>(x)*
 ///
 /// # Examples
 ///
 /// ```
-/// # use boost::math::{hermite, hermite_next};
+/// # use boost::math::{hermite_he, hermite_he_next};
 /// let x = 0.42;
-/// let h0 = hermite(0, x); // 1
-/// let h1 = hermite(1, x); // 2x
-/// let h2 = hermite(2, x); // 4x² - 2
-/// let h3 = hermite(3, x); // 8x³ - 12x
-/// assert_eq!(hermite_next(1, &x, &h1, &h0), h2);
-/// assert_eq!(hermite_next(2, &x, &h2, &h1), h3);
+/// let p0 = hermite_he(0, x); // 1
+/// let p1 = hermite_he(1, x); // x
+/// let p2 = hermite_he(2, x); // x² - 1
+/// let p3 = hermite_he(3, x); // x³ - 3x
+/// assert_eq!(hermite_he_next(1, x, p1, p0), p2);
+/// assert_eq!(hermite_he_next(2, x, p2, p1), p3);
 /// ```
-#[allow(non_snake_case)]
 #[inline(always)]
-pub fn hermite_next(n: u32, x: &f64, Hn: &f64, Hn_1: &f64) -> f64 {
-    2.0 * (x * Hn - (n as f64) * Hn_1)
+pub fn hermite_he_next(n: u32, x: f64, pn: f64, pn_prev: f64) -> f64 {
+    x * pn - (n as f64) * pn_prev
+}
+
+/// Monic Hermite Polynomial *He<sub>n</sub>(x)*
+///
+/// Note that this is the "probabilist's" Hermite polynomial, which is monic (leading coefficient
+/// is 1). It is related to the "physicist's" Hermite polynomial ([`hermite_h`]) by
+///
+/// *He<sub>n</sub>(x) = 2<sup>-n/2</sup> H<sub>n</sub>(x / √2)*
+///
+/// This function does not exist in the Boost Math C++ library.
+pub fn hermite_he(n: u32, x: f64) -> f64 {
+    let (mut p0, mut p1) = (1.0, x);
+    if n == 0 {
+        p0
+    } else {
+        for c in 1..n {
+            (p0, p1) = (p1, hermite_he_next(c, x, p1, p0));
+        }
+        p1
+    }
+}
+
+/// *k*-th derivative of the Hermite polynomial *He<sub>n</sub>(x)*
+///
+/// This function does not exist in the Boost Math C++ library.
+#[doc(alias = "hermite_he_prime")]
+pub fn hermite_he_derivative(n: u32, x: f64, k: u32) -> f64 {
+    if n < k {
+        0.0
+    } else {
+        let mut p = hermite_he(n - k, x);
+        for m in 0..k {
+            p *= (n - m) as f64
+        }
+        p
+    }
 }
 
 #[cfg(test)]
 mod tests {
     use super::*;
+    use core::f64::consts::FRAC_1_SQRT_2;
+
+    const ATOL: f64 = 1e-15;
+    const RTOL: f64 = 1e-12;
+
+    fn hermite_h_ffi(n: u32, x: f64) -> f64 {
+        unsafe { crate::ffi::math_hermite(n as core::ffi::c_uint, x) }
+    }
+
+    #[test]
+    fn test_hermite_h() {
+        for n in 0..20 {
+            for x in -10..=10 {
+                let x = x as f64 * 0.1;
+                let h = hermite_h_ffi(n, x);
+                assert_relative_eq!(hermite_h(n, x), h, epsilon = ATOL, max_relative = RTOL);
+            }
+        }
+    }
+
+    #[test]
+    fn test_hermite_he() {
+        for n in 0..20 {
+            let c = FRAC_1_SQRT_2.powi(n as i32);
+            for x in -10..=10 {
+                let x = x as f64 * 0.1;
+                let h = c * hermite_h(n, FRAC_1_SQRT_2 * x);
+                assert_relative_eq!(hermite_he(n, x), h, epsilon = ATOL, max_relative = RTOL);
+            }
+        }
+    }
 
     #[test]
-    fn test_hermite() {
-        assert_relative_eq!(hermite(0, 1.0), 1.0);
-        assert_relative_eq!(hermite(1, 1.0), 2.0);
-        assert_relative_eq!(hermite(2, 1.0), 2.0);
-        assert_relative_eq!(hermite(3, 1.0), -4.0);
-        assert_relative_eq!(hermite(4, 1.0), -20.0);
+    fn test_hermite_h_derivative() {
+        for n in 0..20 {
+            for x in -10..=10 {
+                let x = x as f64 * 0.1;
+
+                let d0 = hermite_h_derivative(n, x, 0);
+                assert_relative_eq!(d0, hermite_h(n, x));
+
+                let d1 = hermite_h_derivative(n, x, 1);
+                if n >= 1 {
+                    assert_relative_eq!(d1, (2 * n) as f64 * hermite_h(n - 1, x));
+                } else {
+                    assert_eq!(d1, 0.0)
+                }
+
+                let d2 = hermite_h_derivative(n, x, 2);
+                if n >= 2 {
+                    assert_relative_eq!(d2, (4 * n * (n - 1)) as f64 * hermite_h(n - 2, x));
+                } else {
+                    assert_eq!(d2, 0.0)
+                }
+            }
+        }
+    }
+
+    #[test]
+    fn test_hermite_he_derivative() {
+        for n in 0..20 {
+            for x in -10..=10 {
+                let x = x as f64 * 0.1;
+
+                let d0 = hermite_he_derivative(n, x, 0);
+                assert_relative_eq!(d0, hermite_he(n, x));
+
+                let d1 = hermite_he_derivative(n, x, 1);
+                if n >= 1 {
+                    assert_relative_eq!(d1, n as f64 * hermite_he(n - 1, x));
+                } else {
+                    assert_eq!(d1, 0.0)
+                }
+
+                let d2 = hermite_he_derivative(n, x, 2);
+                if n >= 2 {
+                    assert_relative_eq!(d2, (n * (n - 1)) as f64 * hermite_he(n - 2, x));
+                } else {
+                    assert_eq!(d2, 0.0)
+                }
+            }
+        }
     }
 }

src/math/special_functions/jacobi.rs

@@ -3,20 +3,20 @@
 use crate::ffi;
 use core::ffi::c_uint;
 
-/// Jacobi Polynomial *P<sub>n</sub><sup>(α, β)</sup>(x)*
-///
-/// Corresponds to `boost::math::jacobi(n, alpha, beta, x)`.
+/// Jacobi Polynomial *P<sub>n</sub><sup>(α,β)</sup>(x)*
 ///
+/// Corresponds to `boost::math::jacobi(n, alpha, beta, x)` in C++.
 /// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_poly/jacobi.html>
 pub fn jacobi(n: u32, alpha: f64, beta: f64, x: f64) -> f64 {
     unsafe { ffi::math_jacobi(n as c_uint, alpha, beta, x) }
 }
 
 /// *k*-th derivative of [`jacobi`] with respect to `x`
 ///
-/// Corresponds to `boost::math::jacobi_derivative(n, alpha, beta, x, k)`.
-///
+/// Corresponds to `boost::math::jacobi_derivative(n, alpha, beta, x, k)` in C++.
 /// <https://boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/sf_poly/jacobi.html>
+#[doc(alias = "jacobi_prime")]
+#[doc(alias = "jacobi_double_prime")]
 pub fn jacobi_derivative(n: u32, alpha: f64, beta: f64, x: f64, k: u32) -> f64 {
     unsafe { ffi::math_jacobi_derivative(n as c_uint, alpha, beta, x, k as c_uint) }
 }