From 69a2222d53c340ed9831606d0bc255a3b69354e7 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Sun, 9 Mar 2025 19:33:31 +0800 Subject: [PATCH 1/6] Add special::xsf::chbevl Introduces the xsf namespace without tainting the special:: namespace. Introduces the chbevl function, which will be needed for functions such as i0, i1, k0, k1, rgamma and shichi. --- sci-rs/src/special/mod.rs | 4 ++++ sci-rs/src/special/xsf/chbevl.rs | 28 ++++++++++++++++++++++++++++ sci-rs/src/special/xsf/mod.rs | 2 ++ 3 files changed, 34 insertions(+) create mode 100644 sci-rs/src/special/xsf/chbevl.rs create mode 100644 sci-rs/src/special/xsf/mod.rs diff --git a/sci-rs/src/special/mod.rs b/sci-rs/src/special/mod.rs index 5cfce8c9..0f8cc627 100644 --- a/sci-rs/src/special/mod.rs +++ b/sci-rs/src/special/mod.rs @@ -2,3 +2,7 @@ mod combinatorics; pub use combinatorics::*; + +// Name is from special/xsf folder, which is Scipy has designated as X special functions (written +// in C++) that are not exposed to Python. We keep these set of functions as being crate internal. +pub(crate) mod xsf; diff --git a/sci-rs/src/special/xsf/chbevl.rs b/sci-rs/src/special/xsf/chbevl.rs new file mode 100644 index 00000000..9dc7b4f4 --- /dev/null +++ b/sci-rs/src/special/xsf/chbevl.rs @@ -0,0 +1,28 @@ +use num_traits::real::Real; + +/// This function evaluates the series y = Sum{i=0..N} coef[i] T_i(x/2) of Chebyshev polynomials Ti +/// at argument x/2. Returns 0 for empty coef. +/// +/// Coefficients are stored in the reverse order, i.e.: the zero order term is the last in the +/// array. +pub(crate) fn chbevl(x: T, coef: &[T]) -> T +where + T: Real + core::ops::Mul, +{ + // Summation over the empty set is defined to be zero. + if coef.is_empty() { + return T::zero(); + } + + let (b0, _, b2): (T, T, T) = coef.iter().fold( + ( + // Safety:: 0 len is checked above. + *unsafe { coef.first().unwrap_unchecked() }, + T::zero(), + T::zero(), + ), + |acc, &e| (x * acc.0 - acc.1 + e, acc.0, acc.1), + ); + + (b0 - b2) / unsafe { T::from(2.).unwrap_unchecked() } +} diff --git a/sci-rs/src/special/xsf/mod.rs b/sci-rs/src/special/xsf/mod.rs new file mode 100644 index 00000000..8f5747c3 --- /dev/null +++ b/sci-rs/src/special/xsf/mod.rs @@ -0,0 +1,2 @@ +mod chbevl; +pub(crate) use chbevl::*; From c1d498595770ad1d0b5ffb520f42eef49f2f19c3 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Sun, 9 Mar 2025 21:34:59 +0800 Subject: [PATCH 2/6] Add Bessel as a trait: i0 for f64 Implement bessel as a trait for f64 with the modified i0 bessel function first. Can easily be extended for Vec and Array as necessary later. --- sci-rs/src/special/bessel.rs | 14 +++++ sci-rs/src/special/mod.rs | 4 ++ sci-rs/src/special/xsf/i0.rs | 108 ++++++++++++++++++++++++++++++++++ sci-rs/src/special/xsf/mod.rs | 2 + 4 files changed, 128 insertions(+) create mode 100644 sci-rs/src/special/bessel.rs create mode 100644 sci-rs/src/special/xsf/i0.rs diff --git a/sci-rs/src/special/bessel.rs b/sci-rs/src/special/bessel.rs new file mode 100644 index 00000000..52ba4efe --- /dev/null +++ b/sci-rs/src/special/bessel.rs @@ -0,0 +1,14 @@ +/// All [functions located in the `Faster versions of common Bessel +/// functions.`]() +pub trait Bessel { + /// Modified Bessel function of order 0. + /// + /// ## Notes + /// * The range is partitioned into the two intervals [0, 8] and (8, infinity). + /// * [Scipy has this as a + /// ufunc](), + /// as a supposed wrapper over the Cephes routine. We try to define it over reasonable types in + /// the impl. + /// + fn i0(&self) -> Self; +} diff --git a/sci-rs/src/special/mod.rs b/sci-rs/src/special/mod.rs index 0f8cc627..ef410ae8 100644 --- a/sci-rs/src/special/mod.rs +++ b/sci-rs/src/special/mod.rs @@ -3,6 +3,10 @@ mod combinatorics; pub use combinatorics::*; +/// Adds the [Bessel] trait. +mod bessel; +pub use bessel::Bessel; + // Name is from special/xsf folder, which is Scipy has designated as X special functions (written // in C++) that are not exposed to Python. We keep these set of functions as being crate internal. pub(crate) mod xsf; diff --git a/sci-rs/src/special/xsf/i0.rs b/sci-rs/src/special/xsf/i0.rs new file mode 100644 index 00000000..c771e6e6 --- /dev/null +++ b/sci-rs/src/special/xsf/i0.rs @@ -0,0 +1,108 @@ +use super::super::xsf; +use super::super::Bessel; +use num_traits::{real::Real, Pow}; + +// Chebyshev coefficients for exp(-x) I0(x) +// in the interval [0,8]. +// +// lim(x->0){ exp(-x) I0(x) } = 1. +const I0_A_F64: [f64; 30] = [ + -4.41534164647933937950E-18, + 3.33079451882223809783E-17, + -2.43127984654795469359E-16, + 1.71539128555513303061E-15, + -1.16853328779934516808E-14, + 7.67618549860493561688E-14, + -4.85644678311192946090E-13, + 2.95505266312963983461E-12, + -1.72682629144155570723E-11, + 9.67580903537323691224E-11, + -5.18979560163526290666E-10, + 2.65982372468238665035E-9, + -1.30002500998624804212E-8, + 6.04699502254191894932E-8, + -2.67079385394061173391E-7, + 1.11738753912010371815E-6, + -4.41673835845875056359E-6, + 1.64484480707288970893E-5, + -5.75419501008210370398E-5, + 1.88502885095841655729E-4, + -5.76375574538582365885E-4, + 1.63947561694133579842E-3, + -4.32430999505057594430E-3, + 1.05464603945949983183E-2, + -2.37374148058994688156E-2, + 4.93052842396707084878E-2, + -9.49010970480476444210E-2, + 1.71620901522208775349E-1, + -3.04682672343198398683E-1, + 6.76795274409476084995E-1, +]; + +// Chebyshev coefficients for exp(-x) sqrt(x) I0(x) +// in the inverted interval [8,infinity]. +// +// lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi). +const I0_B_F64: [f64; 25] = [ + -7.23318048787475395456E-18, + -4.83050448594418207126E-18, + 4.46562142029675999901E-17, + 3.46122286769746109310E-17, + -2.82762398051658348494E-16, + -3.42548561967721913462E-16, + 1.77256013305652638360E-15, + 3.81168066935262242075E-15, + -9.55484669882830764870E-15, + -4.15056934728722208663E-14, + 1.54008621752140982691E-14, + 3.85277838274214270114E-13, + 7.18012445138366623367E-13, + -1.79417853150680611778E-12, + -1.32158118404477131188E-11, + -3.14991652796324136454E-11, + 1.18891471078464383424E-11, + 4.94060238822496958910E-10, + 3.39623202570838634515E-9, + 2.26666899049817806459E-8, + 2.04891858946906374183E-7, + 2.89137052083475648297E-6, + 6.88975834691682398426E-5, + 3.36911647825569408990E-3, + 8.04490411014108831608E-1, +]; + +impl Bessel for f64 { + fn i0(&self) -> Self { + let x = self.abs(); + if x <= 8. { + let y = (x / 2.) - 2.; + return x.exp() * xsf::chbevl(y, &I0_A_F64); + } + x.exp() * xsf::chbevl(32. / x - 2., &I0_B_F64) / x.sqrt() + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn i0_f64() { + let result: f64 = (0.).i0(); + let exp = 1.; + assert_relative_eq!(result, exp, epsilon = 1e-6); + let result: f64 = (1.).i0(); + let exp = 1.2660658777520082; + assert_relative_eq!(result, exp, epsilon = 1e-6); + let result: f64 = 0.213.i0(); + let exp = 1.0113744522192416; + assert_relative_eq!(result, exp, epsilon = 1e-6); + let result: f64 = 5.0.i0(); + let exp = 27.239871823604442; + assert_relative_eq!(result, exp, epsilon = 1e-6); + let result: f64 = 30.546.i0(); + let exp = 1337209608661.4026; + assert_relative_eq!(result, exp, epsilon = 1e-6); + } +} diff --git a/sci-rs/src/special/xsf/mod.rs b/sci-rs/src/special/xsf/mod.rs index 8f5747c3..e259c768 100644 --- a/sci-rs/src/special/xsf/mod.rs +++ b/sci-rs/src/special/xsf/mod.rs @@ -1,2 +1,4 @@ mod chbevl; pub(crate) use chbevl::*; + +mod i0; From b4e30f5117a9e519129f9216e3845b29eb85a8f3 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Sun, 9 Mar 2025 21:43:22 +0800 Subject: [PATCH 3/6] Add modified bessel i0 for f32 --- sci-rs/src/special/xsf/i0.rs | 93 ++++++++++++++++++++++++++++++++++++ 1 file changed, 93 insertions(+) diff --git a/sci-rs/src/special/xsf/i0.rs b/sci-rs/src/special/xsf/i0.rs index c771e6e6..e3556d3c 100644 --- a/sci-rs/src/special/xsf/i0.rs +++ b/sci-rs/src/special/xsf/i0.rs @@ -82,6 +82,86 @@ impl Bessel for f64 { } } +// Chebyshev coefficients for exp(-x) I0(x) +// in the interval [0,8]. +// +// lim(x->0){ exp(-x) I0(x) } = 1. +const I0_A_F32: [f32; 30] = [ + -4.41534164647933937950E-18, + 3.33079451882223809783E-17, + -2.43127984654795469359E-16, + 1.71539128555513303061E-15, + -1.16853328779934516808E-14, + 7.67618549860493561688E-14, + -4.85644678311192946090E-13, + 2.95505266312963983461E-12, + -1.72682629144155570723E-11, + 9.67580903537323691224E-11, + -5.18979560163526290666E-10, + 2.65982372468238665035E-9, + -1.30002500998624804212E-8, + 6.04699502254191894932E-8, + -2.67079385394061173391E-7, + 1.11738753912010371815E-6, + -4.41673835845875056359E-6, + 1.64484480707288970893E-5, + -5.75419501008210370398E-5, + 1.88502885095841655729E-4, + -5.76375574538582365885E-4, + 1.63947561694133579842E-3, + -4.32430999505057594430E-3, + 1.05464603945949983183E-2, + -2.37374148058994688156E-2, + 4.93052842396707084878E-2, + -9.49010970480476444210E-2, + 1.71620901522208775349E-1, + -3.04682672343198398683E-1, + 6.76795274409476084995E-1, +]; + +// Chebyshev coefficients for exp(-x) sqrt(x) I0(x) +// in the inverted interval [8,infinity]. +// +// lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi). +const I0_B_F32: [f32; 25] = [ + -7.23318048787475395456E-18, + -4.83050448594418207126E-18, + 4.46562142029675999901E-17, + 3.46122286769746109310E-17, + -2.82762398051658348494E-16, + -3.42548561967721913462E-16, + 1.77256013305652638360E-15, + 3.81168066935262242075E-15, + -9.55484669882830764870E-15, + -4.15056934728722208663E-14, + 1.54008621752140982691E-14, + 3.85277838274214270114E-13, + 7.18012445138366623367E-13, + -1.79417853150680611778E-12, + -1.32158118404477131188E-11, + -3.14991652796324136454E-11, + 1.18891471078464383424E-11, + 4.94060238822496958910E-10, + 3.39623202570838634515E-9, + 2.26666899049817806459E-8, + 2.04891858946906374183E-7, + 2.89137052083475648297E-6, + 6.88975834691682398426E-5, + 3.36911647825569408990E-3, + 8.04490411014108831608E-1, +]; + +impl Bessel for f32 { + fn i0(&self) -> Self { + let x = self.abs(); + if x <= 8. { + let y = (x / 2.) - 2.; + return x.exp() * xsf::chbevl(y, &I0_A_F32); + } + x.exp() * xsf::chbevl(32. / x - 2., &I0_B_F32) / x.sqrt() + } +} + #[cfg(test)] mod tests { use super::*; @@ -105,4 +185,17 @@ mod tests { let exp = 1337209608661.4026; assert_relative_eq!(result, exp, epsilon = 1e-6); } + + #[test] + fn i0_f32() { + let result: f32 = (0.).i0(); + let exp = 1.; + assert_relative_eq!(result, exp, epsilon = 1e-6); + let result: f32 = (1.).i0(); + let exp = 1.2660658777520082; + assert_relative_eq!(result, exp, epsilon = 1e-6); + let result: f32 = 0.213.i0(); + let exp = 1.0113744522192416; + assert_relative_eq!(result, exp, epsilon = 1e-6); + } } From ddd0527721cc9892af5621436485f690064d61cb Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Sun, 9 Mar 2025 21:58:43 +0800 Subject: [PATCH 4/6] Add exponentially scaled modified bessel i0e to Bessel trait --- sci-rs/src/special/bessel.rs | 10 +++++++ sci-rs/src/special/xsf/i0.rs | 56 ++++++++++++++++++++++++++++++++++++ 2 files changed, 66 insertions(+) diff --git a/sci-rs/src/special/bessel.rs b/sci-rs/src/special/bessel.rs index 52ba4efe..ab53269a 100644 --- a/sci-rs/src/special/bessel.rs +++ b/sci-rs/src/special/bessel.rs @@ -11,4 +11,14 @@ pub trait Bessel { /// the impl. /// fn i0(&self) -> Self; + + /// Exponentially scaled modified Bessel function of order 0. + /// + /// ## Notes + /// * The range is partitioned into the two intervals [0, 8] and (8, infinity). + /// * [Scipy has this as a + /// ufunc](), + /// as a supposed wrapper over the Cephes routine. We try to define it over reasonable types in + /// the impl. + fn i0e(&self) -> Self; } diff --git a/sci-rs/src/special/xsf/i0.rs b/sci-rs/src/special/xsf/i0.rs index e3556d3c..eae8ef48 100644 --- a/sci-rs/src/special/xsf/i0.rs +++ b/sci-rs/src/special/xsf/i0.rs @@ -80,6 +80,15 @@ impl Bessel for f64 { } x.exp() * xsf::chbevl(32. / x - 2., &I0_B_F64) / x.sqrt() } + + fn i0e(&self) -> Self { + let x = self.abs(); + if x <= 8. { + let y = (x / 2.) - 2.; + return xsf::chbevl(y, &I0_A_F64); + } + xsf::chbevl(32. / x - 2., &I0_B_F64) / x.sqrt() + } } // Chebyshev coefficients for exp(-x) I0(x) @@ -160,6 +169,15 @@ impl Bessel for f32 { } x.exp() * xsf::chbevl(32. / x - 2., &I0_B_F32) / x.sqrt() } + + fn i0e(&self) -> Self { + let x = self.abs(); + if x <= 8. { + let y = (x / 2.) - 2.; + return xsf::chbevl(y, &I0_A_F32); + } + xsf::chbevl(32. / x - 2., &I0_B_F32) / x.sqrt() + } } #[cfg(test)] @@ -198,4 +216,42 @@ mod tests { let exp = 1.0113744522192416; assert_relative_eq!(result, exp, epsilon = 1e-6); } + + #[test] + fn i0e_f64() { + let result: f64 = (0.).i0e(); + let exp = 1.; + assert_relative_eq!(result, exp, epsilon = 1e-6); + let result: f64 = (1.).i0e(); + let exp = 0.46575960759364043; + assert_relative_eq!(result, exp, epsilon = 1e-6); + let result: f64 = 0.213.i0e(); + let exp = 0.8173484705849442; + assert_relative_eq!(result, exp, epsilon = 1e-6); + let result: f64 = 5.0.i0e(); + let exp = 0.18354081260932834; + assert_relative_eq!(result, exp, epsilon = 1e-6); + let result: f64 = 30.546.i0e(); + let exp = 0.0724836816695565; + assert_relative_eq!(result, exp, epsilon = 1e-6); + } + + #[test] + fn i0e_f32() { + let result: f32 = (0.).i0e(); + let exp = 1.; + assert_relative_eq!(result, exp, epsilon = 1e-6); + let result: f32 = (1.).i0e(); + let exp = 0.46575960759364043; + assert_relative_eq!(result, exp, epsilon = 1e-6); + let result: f32 = 0.213.i0e(); + let exp = 0.8173484705849442; + assert_relative_eq!(result, exp, epsilon = 1e-6); + let result: f32 = 5.0.i0e(); + let exp = 0.18354081260932834; + assert_relative_eq!(result, exp, epsilon = 1e-6); + let result: f32 = 30.546.i0e(); + let exp = 0.0724836816695565; + assert_relative_eq!(result, exp, epsilon = 1e-6); + } } From a2b886081fa818ab9b8bad051e0e382472c29be1 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Sun, 9 Mar 2025 22:04:14 +0800 Subject: [PATCH 5/6] Silence clippy excessive precisions in Bessel implementation We keep the extra precision to keep it easy to maintain the implementations (currently just mirrored) between the f32 and f64 versions. --- sci-rs/src/special/xsf/i0.rs | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/sci-rs/src/special/xsf/i0.rs b/sci-rs/src/special/xsf/i0.rs index eae8ef48..56ea50e7 100644 --- a/sci-rs/src/special/xsf/i0.rs +++ b/sci-rs/src/special/xsf/i0.rs @@ -6,6 +6,7 @@ use num_traits::{real::Real, Pow}; // in the interval [0,8]. // // lim(x->0){ exp(-x) I0(x) } = 1. +#[allow(clippy::excessive_precision)] const I0_A_F64: [f64; 30] = [ -4.41534164647933937950E-18, 3.33079451882223809783E-17, @@ -43,6 +44,7 @@ const I0_A_F64: [f64; 30] = [ // in the inverted interval [8,infinity]. // // lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi). +#[allow(clippy::excessive_precision)] const I0_B_F64: [f64; 25] = [ -7.23318048787475395456E-18, -4.83050448594418207126E-18, @@ -95,6 +97,7 @@ impl Bessel for f64 { // in the interval [0,8]. // // lim(x->0){ exp(-x) I0(x) } = 1. +#[allow(clippy::excessive_precision)] const I0_A_F32: [f32; 30] = [ -4.41534164647933937950E-18, 3.33079451882223809783E-17, @@ -132,6 +135,7 @@ const I0_A_F32: [f32; 30] = [ // in the inverted interval [8,infinity]. // // lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi). +#[allow(clippy::excessive_precision)] const I0_B_F32: [f32; 25] = [ -7.23318048787475395456E-18, -4.83050448594418207126E-18, From c01737880988d058b7eef34fbf64ace84436da46 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Mon, 10 Mar 2025 09:40:59 +0800 Subject: [PATCH 6/6] Add Bessel trait to Vec and ArrayD --- sci-rs/src/special/bessel.rs | 54 ++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/sci-rs/src/special/bessel.rs b/sci-rs/src/special/bessel.rs index ab53269a..3c2394b8 100644 --- a/sci-rs/src/special/bessel.rs +++ b/sci-rs/src/special/bessel.rs @@ -1,3 +1,6 @@ +use ndarray::ArrayD; +use num_traits::real::Real; + /// All [functions located in the `Faster versions of common Bessel /// functions.`]() pub trait Bessel { @@ -22,3 +25,54 @@ pub trait Bessel { /// the impl. fn i0e(&self) -> Self; } + +#[cfg(feature = "std")] +impl Bessel for Vec +where + T: Bessel, +{ + fn i0(&self) -> Self { + self.iter().map(Bessel::i0).collect() + } + + fn i0e(&self) -> Self { + self.iter().map(Bessel::i0e).collect() + } +} + +#[cfg(feature = "std")] +impl Bessel for ArrayD +where + T: Bessel, +{ + fn i0(&self) -> Self { + self.map(Bessel::i0) + } + + fn i0e(&self) -> Self { + self.map(Bessel::i0e) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[cfg(feature = "std")] + #[test] + fn i0_vec_f64() { + let inp = vec![0., 1., 0.213, 5., 30.546]; + let result = vec![ + 1., + 1.2660658777520082, + 1.0113744522192416, + 27.239871823604442, + 1337209608661.4026, + ]; + assert_eq!(result.len(), inp.i0().len()); + for (&r, e) in result.iter().zip(inp.i0()) { + assert_relative_eq!(r, e, epsilon = 1e-6); + } + } +}