From 5f38083bfc0a0513d27584c1719c73a16e6049e6 Mon Sep 17 00:00:00 2001 From: ArtemiszenN Date: Wed, 5 Mar 2025 20:01:20 +0800 Subject: [PATCH 1/6] Add Bessel i0 for f64 Bessel i0 function written with reference to Numerical Recipes in C++, currently only with f64. Authored-by: ArtemiszenN Signed-off-by: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> --- sci-rs/src/special/i0.rs | 61 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 sci-rs/src/special/i0.rs diff --git a/sci-rs/src/special/i0.rs b/sci-rs/src/special/i0.rs new file mode 100644 index 00000000..449adecb --- /dev/null +++ b/sci-rs/src/special/i0.rs @@ -0,0 +1,61 @@ +fn poly(cof: &[f64], x: f64) -> f64 { + cof[..cof.len() - 1] + .iter() + .rev() + .fold(*cof.last().expect("Coefficients are empty!"), |acc, e| { + acc * x + e + }) +} + +const I0P: [f64; 14] = [ + 9.999999999999997e-1, + 2.466405579426905e-1, + 1.478980363444585e-2, + 3.826993559940360e-4, + 5.395676869878828e-6, + 4.700912200921704e-8, + 2.733894920915608e-10, + 1.115830108455192e-12, + 3.301093025084127e-15, + 7.209167098020555e-18, + 1.166898488777214e-20, + 1.378948246502109e-23, + 1.124884061857506e-26, + 5.498556929587117e-30, +]; + +const I0Q: [f64; 5] = [ + 4.463598170691436e-1, + 1.702205745042606e-3, + 2.792125684538934e-6, + 2.369902034785866e-9, + 8.965900179621208e-13, +]; + +const I0PP: [f64; 5] = [ + 1.192273748120670e-1, + 1.947452015979746e-1, + 7.629241821600588e-2, + 8.474903580801549e-3, + 2.023821945835647e-4, +]; + +const I0QQ: [f64; 6] = [ + 2.962898424533095e-1, + 4.866115913196384e-1, + 1.938352806477617e-1, + 2.261671093400046e-2, + 6.450448095075585e-4, + 1.529835782400450e-6, +]; + +pub fn i0(x: f64) -> f64 { + let ax = f64::abs(x); + match ax < 15.0 { + true => poly(&I0P, x * x) / poly(&I0Q, 225. - (x * x)), + false => { + f64::exp(ax) * poly(&I0PP, 1.0 - 15.0 / ax) + / (poly(&I0QQ, 1.0 - 15.0 / ax) * f64::sqrt(ax)) + } + } +} From e8bf485b391d3ced1872f060c2cece5a9655593a Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Sat, 8 Mar 2025 00:03:53 +0800 Subject: [PATCH 2/6] Merge i0(f64) code into crate-only function Currently code is not ready for generic float, and requires more work. --- sci-rs/src/special/i0.rs | 28 +++++++++++++++++++++++++++- sci-rs/src/special/mod.rs | 4 ++++ 2 files changed, 31 insertions(+), 1 deletion(-) diff --git a/sci-rs/src/special/i0.rs b/sci-rs/src/special/i0.rs index 449adecb..f6a2a136 100644 --- a/sci-rs/src/special/i0.rs +++ b/sci-rs/src/special/i0.rs @@ -1,3 +1,10 @@ +//! This file provides the modified Bessel function of order zero of the argument. +//! It has been written with reference to Numerical Recipes in C++ by William H. Press, +//! Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. +//! The function as provisioned here are not *yet* intended to be exposed out of the crate. +//! See: https://github.com/scipy/scipy/pull/21297, scipy/special/xsf/cephes/i0.h for a different +//! implementation. + fn poly(cof: &[f64], x: f64) -> f64 { cof[..cof.len() - 1] .iter() @@ -49,7 +56,7 @@ const I0QQ: [f64; 6] = [ 1.529835782400450e-6, ]; -pub fn i0(x: f64) -> f64 { +pub(crate) fn i0(x: f64) -> f64 { let ax = f64::abs(x); match ax < 15.0 { true => poly(&I0P, x * x) / poly(&I0Q, 225. - (x * x)), @@ -59,3 +66,22 @@ pub fn i0(x: f64) -> f64 { } } } + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + fn i0_f64() { + let result: f64 = i0(1.); + let exp = 1.2660658777520082; + assert_relative_eq!(result, exp, epsilon = 1e-6); + let result: f64 = i0(0.213); + let exp = 1.0113744522192416; + assert_relative_eq!(result, exp, epsilon = 1e-6); + let result: f64 = i0(30.546); + let exp = 1337209608661.4026; + assert_relative_eq!(result, exp, epsilon = 1e-6); + } +} diff --git a/sci-rs/src/special/mod.rs b/sci-rs/src/special/mod.rs index 5cfce8c9..f73493e9 100644 --- a/sci-rs/src/special/mod.rs +++ b/sci-rs/src/special/mod.rs @@ -1,4 +1,8 @@ /// Various combinatoric functions for integer-types. mod combinatorics; +/// Modified Bessel function. +mod i0; + pub use combinatorics::*; +pub(crate) use i0::*; From 807741e21b706dcb5ca7ff3b0a1f895316e5831c Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Sat, 8 Mar 2025 14:11:30 +0800 Subject: [PATCH 3/6] Make i0(f64) code more rust-like We also add unsafe bounds since the function is currently only intended to be crate-internal. --- sci-rs/src/special/i0.rs | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/sci-rs/src/special/i0.rs b/sci-rs/src/special/i0.rs index f6a2a136..d01d531b 100644 --- a/sci-rs/src/special/i0.rs +++ b/sci-rs/src/special/i0.rs @@ -5,11 +5,13 @@ //! See: https://github.com/scipy/scipy/pull/21297, scipy/special/xsf/cephes/i0.h for a different //! implementation. -fn poly(cof: &[f64], x: f64) -> f64 { - cof[..cof.len() - 1] - .iter() +/// Evaluates a polynomial with coeffecients evaluated at `x`. +unsafe fn poly(cof: &[f64], x: f64) -> f64 { + cof.iter() + .take(cof.len() - 1) .rev() - .fold(*cof.last().expect("Coefficients are empty!"), |acc, e| { + // SAFETY: poly is called only by const arrays + .fold(unsafe { *cof.last().unwrap_unchecked() }, |acc, e| { acc * x + e }) } @@ -57,13 +59,12 @@ const I0QQ: [f64; 6] = [ ]; pub(crate) fn i0(x: f64) -> f64 { - let ax = f64::abs(x); + let ax = x.abs(); match ax < 15.0 { - true => poly(&I0P, x * x) / poly(&I0Q, 225. - (x * x)), - false => { - f64::exp(ax) * poly(&I0PP, 1.0 - 15.0 / ax) - / (poly(&I0QQ, 1.0 - 15.0 / ax) * f64::sqrt(ax)) - } + true => unsafe { poly(&I0P, x * x) / poly(&I0Q, 225. - (x * x)) }, + false => unsafe { + ax.exp() * poly(&I0PP, 1.0 - 15.0 / ax) / (poly(&I0QQ, 1.0 - 15.0 / ax) * ax.sqrt()) + }, } } From 9138e5e2e8e3dafc76c63d6e61eb2d1523c5482f Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Sat, 8 Mar 2025 14:41:56 +0800 Subject: [PATCH 4/6] Change Bessel into trait This will allow for f32 implementation for generics later. --- sci-rs/src/special/i0.rs | 42 ++++++++++++++++++++++++++++------------ 1 file changed, 30 insertions(+), 12 deletions(-) diff --git a/sci-rs/src/special/i0.rs b/sci-rs/src/special/i0.rs index d01d531b..3cda9c73 100644 --- a/sci-rs/src/special/i0.rs +++ b/sci-rs/src/special/i0.rs @@ -5,6 +5,33 @@ //! See: https://github.com/scipy/scipy/pull/21297, scipy/special/xsf/cephes/i0.h for a different //! implementation. +/// All [functions located in the `Faster versions of common Bessel +/// functions.`]() +pub(crate) 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; +} + +impl Bessel for f64 { + fn i0(&self) -> Self { + let x = *self; + let ax = x.abs(); + match ax < 15.0 { + true => unsafe { poly(&I0P, x * x) / poly(&I0Q, 225. - (x * x)) }, + false => unsafe { + ax.exp() * poly(&I0PP, 1.0 - 15.0 / ax) / (poly(&I0QQ, 1.0 - 15.0 / ax) * ax.sqrt()) + }, + } + } +} + /// Evaluates a polynomial with coeffecients evaluated at `x`. unsafe fn poly(cof: &[f64], x: f64) -> f64 { cof.iter() @@ -58,15 +85,6 @@ const I0QQ: [f64; 6] = [ 1.529835782400450e-6, ]; -pub(crate) fn i0(x: f64) -> f64 { - let ax = x.abs(); - match ax < 15.0 { - true => unsafe { poly(&I0P, x * x) / poly(&I0Q, 225. - (x * x)) }, - false => unsafe { - ax.exp() * poly(&I0PP, 1.0 - 15.0 / ax) / (poly(&I0QQ, 1.0 - 15.0 / ax) * ax.sqrt()) - }, - } -} #[cfg(test)] mod tests { @@ -75,13 +93,13 @@ mod tests { #[test] fn i0_f64() { - let result: f64 = i0(1.); + let result: f64 = (1.).i0(); let exp = 1.2660658777520082; assert_relative_eq!(result, exp, epsilon = 1e-6); - let result: f64 = i0(0.213); + let result: f64 = 0.213.i0(); let exp = 1.0113744522192416; assert_relative_eq!(result, exp, epsilon = 1e-6); - let result: f64 = i0(30.546); + let result: f64 = 30.546.i0(); let exp = 1337209608661.4026; assert_relative_eq!(result, exp, epsilon = 1e-6); } From e441b7aa33372e987ae7a3c8512d6cde9ba950e2 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Sat, 8 Mar 2025 14:42:24 +0800 Subject: [PATCH 5/6] Implement Bessel for f32 too This is very hacky and just copies all the const arrays from f64 to f32. Also not as accurate for very large arguments such as x>~20 being passed to i0. (Correspondingly similar unit test for ~30 is not being tested in f32). --- sci-rs/src/special/i0.rs | 76 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 74 insertions(+), 2 deletions(-) diff --git a/sci-rs/src/special/i0.rs b/sci-rs/src/special/i0.rs index 3cda9c73..f84acd8b 100644 --- a/sci-rs/src/special/i0.rs +++ b/sci-rs/src/special/i0.rs @@ -5,8 +5,11 @@ //! See: https://github.com/scipy/scipy/pull/21297, scipy/special/xsf/cephes/i0.h for a different //! implementation. +use num_traits::real::Real; + /// All [functions located in the `Faster versions of common Bessel /// functions.`]() +/// Note strongly that Bessel has tremendous inacurracy for f32 compared to f64. pub(crate) trait Bessel { /// Modified Bessel function of order 0. /// @@ -32,13 +35,31 @@ impl Bessel for f64 { } } +impl Bessel for f32 { + // Known to yield wrong result. + fn i0(&self) -> Self { + let x = *self; + let ax = x.abs(); + match ax < 15.0 { + true => unsafe { poly(&I0P_F32, x * x) / poly(&I0Q_F32, 225. - (x * x)) }, + false => unsafe { + ax.exp() * poly(&I0PP_F32, 1.0 - 15.0 / ax) + / (poly(&I0QQ_F32, 1.0 - 15.0 / ax) * ax.sqrt()) + }, + } + } +} + /// Evaluates a polynomial with coeffecients evaluated at `x`. -unsafe fn poly(cof: &[f64], x: f64) -> f64 { +unsafe fn poly(cof: &[T], x: T) -> T +where + T: Real + std::ops::Mul, +{ cof.iter() .take(cof.len() - 1) .rev() // SAFETY: poly is called only by const arrays - .fold(unsafe { *cof.last().unwrap_unchecked() }, |acc, e| { + .fold(unsafe { *cof.last().unwrap_unchecked() }, |acc, &e| { acc * x + e }) } @@ -60,6 +81,23 @@ const I0P: [f64; 14] = [ 5.498556929587117e-30, ]; +const I0P_F32: [f32; 14] = [ + 9.999999999999997e-1, + 2.466405579426905e-1, + 1.478980363444585e-2, + 3.826993559940360e-4, + 5.395676869878828e-6, + 4.700912200921704e-8, + 2.733894920915608e-10, + 1.115830108455192e-12, + 3.301093025084127e-15, + 7.209167098020555e-18, + 1.166898488777214e-20, + 1.378948246502109e-23, + 1.124884061857506e-26, + 5.498556929587117e-30, +]; + const I0Q: [f64; 5] = [ 4.463598170691436e-1, 1.702205745042606e-3, @@ -68,6 +106,14 @@ const I0Q: [f64; 5] = [ 8.965900179621208e-13, ]; +const I0Q_F32: [f32; 5] = [ + 4.463598170691436e-1, + 1.702205745042606e-3, + 2.792125684538934e-6, + 2.369902034785866e-9, + 8.965900179621208e-13, +]; + const I0PP: [f64; 5] = [ 1.192273748120670e-1, 1.947452015979746e-1, @@ -76,6 +122,14 @@ const I0PP: [f64; 5] = [ 2.023821945835647e-4, ]; +const I0PP_F32: [f32; 5] = [ + 1.192273748120670e-1, + 1.947452015979746e-1, + 7.629241821600588e-2, + 8.474903580801549e-3, + 2.023821945835647e-4, +]; + const I0QQ: [f64; 6] = [ 2.962898424533095e-1, 4.866115913196384e-1, @@ -85,6 +139,14 @@ const I0QQ: [f64; 6] = [ 1.529835782400450e-6, ]; +const I0QQ_F32: [f32; 6] = [ + 2.962898424533095e-1, + 4.866115913196384e-1, + 1.938352806477617e-1, + 2.261671093400046e-2, + 6.450448095075585e-4, + 1.529835782400450e-6, +]; #[cfg(test)] mod tests { @@ -103,4 +165,14 @@ mod tests { let exp = 1337209608661.4026; assert_relative_eq!(result, exp, epsilon = 1e-6); } + + #[test] + fn i0_f32() { + 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 90b9585df008eb5c4b6f4e2a7efc9cfb91c23012 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Sat, 8 Mar 2025 14:51:22 +0800 Subject: [PATCH 6/6] Use core::ops::Mul over std::ops::Mul for Bessel Unit test caught this issue --- sci-rs/src/special/i0.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sci-rs/src/special/i0.rs b/sci-rs/src/special/i0.rs index f84acd8b..e4d9ac12 100644 --- a/sci-rs/src/special/i0.rs +++ b/sci-rs/src/special/i0.rs @@ -53,7 +53,7 @@ impl Bessel for f32 { /// Evaluates a polynomial with coeffecients evaluated at `x`. unsafe fn poly(cof: &[T], x: T) -> T where - T: Real + std::ops::Mul, + T: Real + core::ops::Mul, { cof.iter() .take(cof.len() - 1)