Complex-bessel — Pure Rust Bessel, Hankel, and Airy functions

announcements
1

Hi everyone,

I'd like to share a crate I recently published: complex-bessel.

It's a pure Rust implementation of the Amos (TOMS 644) algorithm for computing Bessel, Hankel, and Airy functions of complex argument and real order. I previously relied on a Fortran-based wrapper crate, but ran into difficulties getting gfortran to work with the MSVC toolchain on Windows. That led me to rewrite the algorithm entirely in Rust.

Features:

  • Bessel functions J, Y, I, K and Hankel functions H1, H2
  • Airy functions Ai, Bi
  • No Fortran/C FFI dependencies
  • no_std support

Feedback and issues are welcome!

5 Likes
2

This is pretty impressive. There's no T::mul_add in your crate and I'm curious: I think rustc will never convert a*b + c to fma, but ifort and icc can and will under certain settings. Are there any places (e.g. zseri) where fused multiply-add can save some flops?

1 Like
3

Thanks for the insightful suggestion! You clearly have deep expertise in this area.

You're right, rustc won't auto-convert a*b+c to FMA like ifort/icc can. I'm planning to add explicit mul_add support:

  • Add a mul_add method to the internal BesselFloat trait
  • std builds → hardware FMA via Float::mul_add
  • no_std builds → fallback to plain a * b + c (to avoid the slow libm::fma software emulation)
  • Apply it across scalar and complex arithmetic throughout the codebase

This should bring the Rust implementation closer to what Fortran gets automatically, both in precision and performance. Thanks again for pointing this out!

4

I've applied FMA on the dev branch based on your suggestion. Tested with my benchmark suite — seeing about 4~18% speedup depending on the function, with no change in accuracy. Thanks for the great tip!

1 Like
5

Correct, it explicitly doesn't. But you might want to follow

which is about adding an opt-in version that can, if you're ok with the corresponding nondeterminism.

2 Likes
6

Thanks for the pointer! That's exactly what I need — right now I'm maintaining a BesselFloat trait that branches between hardware mul_add (std) and plain a * b + c (no_std) to avoid the slow software FMA fallback. float_mul_add_relaxed would let me drop that conditional logic entirely and just use a single code path. I'll keep an eye on the tracking issue.

1 Like