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.2"

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

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