Modified Bessel Function of the first kind of order zero

[jbe]

@stefan-k

In argmin I created a new trait ArgminFloat with all the trait bounds that I found reasonable for my use case and implemented this trait for all types which fulfill these trait bounds. I think this is similar to the approach in rustfft. You can find the implementation here.

So practically a "trait alias" (see also nightly feature) for all the traits you need. It seems to be common practice then, that each provider of floating point functions will provide their own trait, which specifies which functionality for floating point types is needed for a particular library.

I could create one for Float + NumAssignOps, I guess. But maybe desugaring x += y to x = x + y is not a big problem in my case. However, += (AddAssign) works by reference on the left hand side and might thus be faster for some implementations of Floats, I guess? That would be a good reason to provide my own trait.

@stefan-k

In linfa they did it in a similar way, but they implemented their Float trait on f32 and f64 explicitly and they added a convenient cast method.

I think there might be an advantage to implement the trait on f32 and f64 explicitly (and maybe even make it a sealed trait?). If you later need to add more supertraits to your own trait, this would break semver compatibility. See also what @afetisov said above:

@afetisov

(new methods on primitives are regularly added, are they also added to num-traits?)

I didn't consider semver compatibility in that matter.

Edit: Maybe a reasonable solution would be to simply document that new bounds might be added to whichever Float trait (alias) is being used, and that this doesn't break semver compatibility (i.e. you can only rely on f32 and f64 implementing it).

I see three and a half ways to go (which may be different for the API and internal/private functions):

• Not use generics and only work on f32 and f64 directly
  • either by duplicating code manually
  • or through macros;
• Using a trait alias (or trait with blanket implementation on stable Rust) that combines something like num_traits::float::Float with any other traits that are needed for a particular use case;
• Using an own (possibly sealed?) trait that's directly implemented for f32 and f64.

Hmmmm… lots of possible trouble. Not sure what I should do.

Side question:

@stefan-k

argmin

Does it also contain a (multivariate) Estimation of Distribution Algorithm as one possible strategy, and what do you think about those algorithms?

---

@stefan-k

A simple float!(..) macro avoids having to write F::from(...).unwrap() everywhere (which still isn't great, but at least it leads to fewer distracting unwraps).

If I understand it right, the macro requires the type to be named F (source). I think that should be documented (currently doesn't seem to be documented). Maybe I'll create a private macro for my purpose. Something like:

use num_traits::Float;

macro_rules! flt {
    ($x:expr) => {
        Flt::from($x).unwrap()
    };
}

#[allow(non_snake_case)]
pub fn bessel_I0<Flt: Float>(x: Flt) -> Flt {
    let base = x * x / flt!(4);
    let mut addend = Flt::one();
    let mut sum = Flt::one();
    for j in 1.. {
        addend = addend * base / flt!(j * j);
        let old = sum;
        sum = sum + addend;
        if sum == old || !sum.is_finite() {
            break;
        }
    }
    sum
}

fn main() {
    for x in [-1.0, 0.0, 0.5, 1.0, 2.0, 13.4, 500.0, 1000.0] {
        println!("I_0({}) = {}", x, bessel_I0(x));
    }
}

(Playground)