Change Bessel interface to include for more types

By using as a trait, we can target between f32, f64, Vec<F> and even
Array<F, Dimensions>.

A bug with the f32 implementation of i0 was found, and is ignored in the
tests.

Author: SpookyYomo

sci-rs/Cargo.toml

@@ -28,6 +28,7 @@ std = ['nalgebra/std', 'nalgebra/macros', 'rustfft', 'alloc']
 plot = ['std']
 
 [dependencies]
+special-fun = { version = "0.3.0", default-features = false }
 num-traits = { version = "0.2.15", default-features = false }
 itertools = { version = "0.13.0", default-features = false }
 nalgebra = { version = "0.33.2", default-features = false }

sci-rs/src/signal/windows/kaiser.rs

@@ -1,5 +1,5 @@
 use super::{extend, len_guard, truncate};
-use crate::special::i0;
+use crate::special::Bessel;
 use num_traits::{real::Real, Float};
 
 #[cfg(feature = "alloc")]
@@ -49,8 +49,7 @@ where
 impl<F, W> GetWindow<W> for Kaiser<F>
 where
     F: Real,
-    W: Real,
-    f64: From<W>,
+    W: Real + Bessel,
 {
     /// Return a [Kaiser] window.
     ///
@@ -189,11 +188,12 @@ where
         //     .collect();
         let w: Vec<W> = n
             .map(|ni| {
-                i0(beta
+                (beta
                     * (W::one()
                         - ((W::from(ni).unwrap() - alpha) / alpha).powf(W::from(2).unwrap()))
                     .sqrt())
-                    / i0(beta)
+                .i0()
+                    / (beta.i0())
             })
             .collect();
         return truncate(w, needs_trunc);

sci-rs/src/special/faster_bessel.rs

@@ -1,31 +1,101 @@
+use ndarray::{Array, Dimension};
 use num_traits::real::Real;
 
-/// 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](<https://docs.scipy.org/doc/scipy/reference/special.html#special-functions-scipy-special>),
-/// as a supposed wrapper over the Cephes routine.
-pub fn i0<F>(x: F) -> F
+#[cfg(feature = "alloc")]
+use alloc::vec::Vec;
+
+/// All [functions located in the `Faster versions of common Bessel
+/// functions.`](<https://docs.scipy.org/doc/scipy/reference/special.html#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](<https://docs.scipy.org/doc/scipy/reference/special.html#special-functions-scipy-special>),
+    /// as a supposed wrapper over the Cephes routine. We try to define it over reasonable types in
+    /// 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](<https://docs.scipy.org/doc/scipy/reference/special.html#special-functions-scipy-special>),
+    /// as a supposed wrapper over the Cephes routine. We try to define it over reasonable types in
+    /// the impl.
+    fn i0e(&self) -> Self;
+}
+
+impl Bessel for f32 {
+    // Known to yield wrong result.
+    fn i0(&self) -> Self {
+        return unsafe { special_fun::unsafe_cephes_single::i0f(*self) };
+    }
+
+    fn i0e(&self) -> Self {
+        return unsafe { special_fun::unsafe_cephes_single::i0ef(*self) };
+    }
+}
+
+impl Bessel for f64 {
+    fn i0(&self) -> Self {
+        return unsafe { special_fun::unsafe_cephes_double::i0(*self) };
+    }
+
+    fn i0e(&self) -> Self {
+        return unsafe { special_fun::unsafe_cephes_double::i0e(*self) };
+    }
+}
+
+#[cfg(feature = "alloc")]
+impl<F> Bessel for Vec<F>
 where
-    F: Real,
-    f64: From<F>,
+    F: Real + super::Bessel,
 {
-    return F::from(unsafe { special_fun::unsafe_cephes_double::i0(f64::from(x)) }).unwrap();
+    fn i0(&self) -> Self {
+        self.iter().map(|f| f.i0()).collect()
+    }
+
+    fn i0e(&self) -> Self {
+        self.iter().map(|f| f.i0e()).collect()
+    }
 }
 
-/// 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](<https://docs.scipy.org/doc/scipy/reference/special.html#special-functions-scipy-special>),
-/// as a supposed wrapper over the Cephes routine.
-pub fn i0e<F>(x: F) -> F
+#[cfg(feature = "alloc")]
+impl<F, D> Bessel for Array<F, D>
 where
-    F: Real,
-    f64: From<F>,
+    F: Real + super::Bessel,
+    D: Dimension,
 {
-    return F::from(unsafe { special_fun::unsafe_cephes_double::i0e(f64::from(x)) }).unwrap();
+    fn i0(&self) -> Self {
+        self.map(|f| f.i0())
+    }
+
+    fn i0e(&self) -> Self {
+        self.map(|f| f.i0e())
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use approx::assert_abs_diff_eq;
+
+    #[test]
+    #[ignore = "upstream"]
+    fn i0single() {
+        let inp: f32 = 0.213;
+        // upstream (https://www.moshier.net/#Cephes) has to fix
+        assert_abs_diff_eq!(1.0113744522192416, inp.i0());
+        assert_abs_diff_eq!(0.8173484705849442, inp.i0e());
+    }
+
+    #[test]
+    fn i0double() {
+        let inp: f64 = 0.213;
+        assert_abs_diff_eq!(1.0113744522192416, inp.i0());
+        assert_abs_diff_eq!(0.8173484705849442, inp.i0e());
+    }
 }