fortran-bessels is a Modern Fortran (2008+) port of Bessels.jl, built for one thing: evaluating Bessel functions as fast and as accurately as possible in pure Fortran. The fixed-order kernels beat the gfortran intrinsics by roughly 2× and the netlib specfun reference by orders of magnitude, while matching them to a few units in the last place (ULP).

Every function is an elemental real(BK) function (BK = real64), so it accepts a scalar or a whole array interchangeably — no wrapper loops, no temporaries.

Function reference

First kind

Bessel Functions of the First Kind: J — \( J_0,\; J_1,\; J_n \) — oscillatory, finite at the origin

Second kind

Bessel Functions of the Second Kind: Y — \( Y_0,\; Y_1,\; Y_\nu \) — oscillatory, singular at the origin

Modified, first kind

Modified Bessel Functions of the First Kind: I — \( I_0,\; I_1 \) — exponentially growing

Modified, second kind

Modified Bessel Functions of the Second Kind: K — \( K_0,\; K_1 \) — exponentially decaying

Hankel

Hankel Functions — \( H^{(1)}_\nu,\; H^{(2)}_\nu \) — complex-valued

Support

Gamma Function — \( \Gamma(x) \)

Quick start

use bessels, only: besselj0, besselk1, bk
 
real(BK) :: x(4)
x = [0.5_bk, 1.0_bk, 2.0_bk, 4.0_bk]
 
print *, besselj0(x)        ! elemental: acts on the whole array
print *, besselk1(2.0_bk)

Building

The canonical build is fpm:

fpm test --profile release --flag "-march=native"

See the project README for the full benchmark table and the manual gfortran recipe.

Warning: The variable-order routines bessely(nu,x), besselh, hankelh1, and hankelh2 are validated for integer order \( \nu \) only. Non-integer order has known port bugs in the variable-order Y/Hankel paths and will return incorrect values. See roadmap item 07.

Roadmap

Still to come: real-order \( I_\nu,\; K_\nu \), spherical Bessels, and the Airy functions. See the implementation roadmap.