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, Ai′, Bi, Bi′
Consecutive orders — _seq and _seq_into variants evaluate at orders ν, ν+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 — no allocator needed; only the Vec-returning _seq variants require std

Accuracy & performance

Relative error against mpmath for complex-bessel (Rust), Amos TOMS 644 (Fortran), SciPy (Python), and SpecialFunctions.jl (Julia).

Evaluation time per call.

The charts are for f64. The test repository has the methodology, full results, and f32 accuracy.

Quick start

[dependencies]
complex-bessel = "0.5"

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)besselj_seq_into(ν, z, scaling, out)`	Jν(z), Bessel first kind	exp(−\|Im(z)\|) · Jν(z)
`bessely(ν, z)bessely_scaled(ν, z)bessely_seq(ν, z, n, scaling)bessely_seq_into(ν, z, scaling, out)`	Yν(z), Bessel second kind	exp(−\|Im(z)\|) · Yν(z)
`besseli(ν, z)besseli_scaled(ν, z)besseli_seq(ν, z, n, scaling)besseli_seq_into(ν, z, scaling, out)`	Iν(z), modified first kind	exp(−\|Re(z)\|) · Iν(z)
`besselk(ν, z)besselk_scaled(ν, z)besselk_seq(ν, z, n, scaling)besselk_seq_into(ν, z, scaling, out)`	Kν(z), modified second kind	exp(z) · Kν(z)
`hankel1(ν, z)hankel1_scaled(ν, z)hankel1_seq(ν, z, n, scaling)hankel1_seq_into(ν, z, scaling, out)`	Hν(1)(z), Hankel first kind	exp(−iz) · Hν(1)(z)
`hankel2(ν, z)hankel2_scaled(ν, z)hankel2_seq(ν, z, n, scaling)hankel2_seq_into(ν, z, scaling, out)`	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 _seq_into variants (besselj_seq_into, besselk_seq_into, …) write the sequence into a caller-provided slice and return a SeqInfo with Accuracy; they never allocate. 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 the accuracy status, use a _seq/_seq_into function (Bessel) or a _raw function (Airy):

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

let z = Complex::new(1.0, 2.0);
let mut out = [Complex::new(0.0, 0.0); 1];
let info = besselk_seq_into(0.0, z, Scaling::Unscaled, &mut out).unwrap();
assert!(matches!(info.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, H, Y (except at negative half-integer orders) and for J, I of negative non-integer order; n < 1 (empty `out`)
`Overflow`	\|z\| or \|ν\| too large (or \|z\| too small) for a finite result
`TotalPrecisionLoss`	\|z\| or \|ν\| too large for meaningful computation
`ConvergenceFailure`	internal algorithm did not converge

Using without `std`

Turn off the default std feature:

complex-bessel = { version = "0.5", default-features = false }

All functions except the six _seq variants (which return a Vec) remain available, and none of them allocate. Math functions then come from libm instead of the platform library.

Improvements over Amos

Scaled Y at large |Im z| — Amos's ZBESY drops a dominant term as an underflow (large ν, |Im z| > ~350); complex-bessel keeps it.
No spurious overflow in K and H — double-precision ZUOIK reports Overflow for representable values near 1e300; complex-bessel follows the single-precision CUOIK test.
Complex division — Smith's algorithm instead of ZDIV, for better f32 accuracy.
Airy limit in uniform expansions — at |z| ≈ 7e8–1e9 and large ν, ZUNI2/ZUNK2 use an undefined result when the Airy argument exceeds its range limit; complex-bessel evaluates the Airy function there anyway.
Small components on the axes — where a value combines two real functions of very different size (e.g. H = J + iY on the real axis), Amos's small component can lose all its digits; complex-bessel recomputes it from the real functions.
Phase of scaled K and H — in the analytic continuation, Amos's ZS1S2 rounds the phase of exp(−2i Im z), a relative error up to ε·|z|; complex-bessel applies this factor separately, keeping the phase exact.

The last two show up in the componentwise relative error (of the real and imaginary parts separately):

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.5.0