complex-bessel

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

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, alloc, std (default)

Quick start

[dependencies]
complex-bessel = "0.4"

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

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`	non-finite z or ν, z = 0 for K/Y/H, n < 1
`Overflow`	\|z\| or \|ν\| too large (or \|z\| too small) for finite result
`TotalPrecisionLoss`	\|z\| or \|ν\| too large for meaningful computation
`ConvergenceFailure`	Internal algorithm did not converge

`no_std` support

Cargo features	Available API
`default-features = false`	24 single-value functions
`features = ["alloc"]`	full API (+ 6 `_seq` variants)
`features = ["std"]` (default)	full API, with platform-native math

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.4", default-features = false }

# no_std + alloc — full API:
complex-bessel = { version = "0.4", default-features = false, features = ["alloc"] }

Differences from TOMS 644

Scaled Y — keeps a term that ZBESY wrongly drops as an underflow (large ν, |Im z| > ~350).
K/H overflow test — follows single-precision CUOIK; double ZUOIK reports spurious Overflow near 1e300.
Complex division — Smith's algorithm instead of ZDIV, for better f32 accuracy.

Accuracy & performance

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

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

Median evaluation time per function call.

The charts are for f64. For methodology, datasets, full benchmark results, and f32 accuracy, see the test repository.

License

Licensed under either of

Apache License, Version 2.0 (LICENSE-APACHE or http://www.apache.org/licenses/LICENSE-2.0)
MIT License (LICENSE-MIT or http://opensource.org/licenses/MIT)

at your option.

complex-bessel 0.4.0