complex-bessel 0.1.0

Pure Rust implementation of complex Bessel functions (J, Y, I, K, H, Airy) based on Amos Algorithm 644
Documentation

complex-bessel

Pure Rust implementation of complex Bessel functions based on Amos Algorithm 644 (ACM TOMS 644).

crates.io docs.rs license CI

Features

  • f32 & f64 — all functions accept Complex<f64> or Complex<f32>
  • Complete function set — J, Y, I, K, H(1), H(2), Ai, Bi
  • Consecutive orders_seq variants return ν, ν+1, …, ν+n−1 in one call
  • Exponential scaling_scaled variants prevent overflow/underflow
  • Negative orders — supports ν < 0 via DLMF reflection formulas (not in Amos)
  • no_std support — 3-tier: bare no_std (no allocator), alloc, std (default)

Quick start

[dependencies]
complex-bessel = "0.1"
use complex_bessel::*;
use num_complex::Complex;

let z = Complex::new(1.0, 2.0);

let j = besselj(0.5, z).unwrap();
let k = besselk(1.0, z).unwrap();
let h = hankel1(0.0, z).unwrap();

// Scaled versions prevent overflow/underflow
let k_scaled = besselk_scaled(1.0, z).unwrap();

// Airy functions
let ai = airy(z).unwrap();
let ai_prime = airyprime(z).unwrap();

Functions

Function Description Scaled variant returns
besselj(ν, z)besselj_scaled(ν, z)besselj_seq(ν, z, n, scaling) Jν(z), Bessel first kind exp(−|Im(z)|) · Jν(z)
bessely(ν, z)bessely_scaled(ν, z)bessely_seq(ν, z, n, scaling) Yν(z), Bessel second kind exp(−|Im(z)|) · Yν(z)
besseli(ν, z)besseli_scaled(ν, z)besseli_seq(ν, z, n, scaling) Iν(z), modified first kind exp(−|Re(z)|) · Iν(z)
besselk(ν, z)besselk_scaled(ν, z)besselk_seq(ν, z, n, scaling) Kν(z), modified second kind exp(z) · Kν(z)
hankel1(ν, z)hankel1_scaled(ν, z)hankel1_seq(ν, z, n, scaling) Hν(1)(z), Hankel first kind exp(−iz) · Hν(1)(z)
hankel2(ν, z)hankel2_scaled(ν, z)hankel2_seq(ν, z, n, scaling) Hν(2)(z), Hankel second kind exp(iz) · Hν(2)(z)
airy(z)airy_scaled(z)airy_raw(z, scaling) Ai(z), Airy first kind exp(ζ) · Ai(z)
airyprime(z)airyprime_scaled(z)airyprime_raw(z, scaling) Ai′(z), derivative of Airy first kind exp(ζ) · Ai′(z)
biry(z)biry_scaled(z)biry_raw(z, scaling) Bi(z), Airy second kind exp(−|Re(ζ)|) · Bi(z)
biryprime(z)biryprime_scaled(z)biryprime_raw(z, scaling) Bi′(z), derivative of Airy second kind exp(−|Re(ζ)|) · Bi′(z)

where ζ = (2/3) z√z.

Function variants

The _seq variants (besselj_seq, besselk_seq, …) compute values at consecutive orders ν, ν+1, …, ν+n−1 in a single call and return a BesselResult that includes an Accuracy field. The _raw Airy variants (airy_raw, biry_raw, …) similarly return an AiryResult with Accuracy. All other functions return only the computed value without Accuracy.

Accuracy:

Status Meaning
Normal Full machine precision
Reduced More than half of significant digits may be lost.Occurs only when |z| or ν exceeds ~32767

Reduced is extremely rare in practice. SciPy's Bessel wrappers also silently discard the equivalent Amos IERR=3 flag by default.

To check accuracy status, use a _seq function (Bessel) or a _raw function (Airy):

use complex_bessel::*;
use num_complex::Complex;

let z = Complex::new(1.0, 2.0);
let result = besselk_seq(0.0, z, 1, Scaling::Unscaled).unwrap();
assert!(matches!(result.status, Accuracy::Normal));

let result = airy_raw(z, Scaling::Unscaled).unwrap();
assert!(matches!(result.status, Accuracy::Normal));

Error handling

All functions return Result<_, Error>. The four error variants are:

Variant Cause
InvalidInput z = 0 for K/Y/H, n < 1
Overflow |z| or ν too large (or too small) for finite result
TotalPrecisionLoss Complete loss of significant digits; |z| or ν too large
ConvergenceFailure Internal algorithm did not converge

Error implements Display always and std::error::Error with the std feature.

no_std support

Cargo features Available API
default-features = false 24 single-value functions
features = ["alloc"] + 6 _seq variants + BesselResult
features = ["std"] (default) + impl Error for Error

The 24 single-value functions include 12 Bessel (J/Y/I/K/H(1)/H(2) × unscaled/scaled) and 12 Airy (Ai/Ai'/Bi/Bi' × unscaled/scaled/raw).

# Bare no_std — no allocator needed:
complex-bessel = { version = "0.1", default-features = false }

# no_std with alloc (adds _seq functions and BesselResult):
complex-bessel = { version = "0.1", default-features = false, features = ["alloc"] }

Accuracy & performance

Precision of complex-bessel (Rust), AMOS/TOMS 644 (Fortran), and SciPy (Python) across all functions.

Accuracy — relative error vs mpmath

Since this library is a Rust translation of AMOS/TOMS 644, the two should produce nearly identical results.

Rust vs Fortran — relative error

Median evaluation time per function call.

Performance — median time per call

For methodology, datasets, and full benchmark results, see the test repository.

License

Licensed under either of

at your option.