Crate complex_bessel

Expand description

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⁽¹⁾, H⁽²⁾, 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

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_ν⁽¹⁾(z), Hankel first kind	exp(−iz) · H_ν⁽¹⁾(z)
`hankel2(ν, z)` `hankel2_scaled(ν, z)` `hankel2_seq(ν, z, n, scaling)`	H_ν⁽²⁾(z), Hankel second kind	exp(iz) · H_ν⁽²⁾(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.

§Exponential scaling

Each _scaled variant multiplies the result by an exponential factor that cancels the asymptotic growth (or decay), keeping values in a representable floating-point range. The exact factor is listed in the “Scaled variant returns” column above. Use scaled functions whenever |z| or |Im(z)| is large enough that the unscaled result would overflow or underflow.

§Negative orders

All functions (single-value and _seq variants) accept any real order, including negative values. DLMF reflection formulas are applied automatically:

J: J_−ν(z) = cos(νπ) J_ν(z) − sin(νπ) Y_ν(z) (DLMF 10.2.3)
Y: Y_−ν(z) = sin(νπ) J_ν(z) + cos(νπ) Y_ν(z) (DLMF 10.2.3)
I: I_−ν(z) = I_ν(z) + (2/π) sin(νπ) K_ν(z) (DLMF 10.27.2)
K: K_−ν(z) = K_ν(z) (even in ν, DLMF 10.27.3)
H⁽¹⁾: H⁽¹⁾_−ν(z) = exp(νπi) H⁽¹⁾_ν(z) (DLMF 10.4.6)
H⁽²⁾: H⁽²⁾_−ν(z) = exp(−νπi) H⁽²⁾_ν(z) (DLMF 10.4.6)

For integer orders, simplified identities are used (e.g., J_−n(z) = (−1)ⁿ J_n(z)).

§Function variants

§Consecutive orders

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⁽¹⁾/H⁽²⁾ × 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"] }

Re-exports§

pub use machine::BesselFloat;
pub use types::BesselResult;
pub use types::Accuracy;
pub use types::AiryResult;
pub use types::Error;
pub use types::Scaling;

Modules§

machine: Machine constants and the BesselFloat trait.
types: Core types for Bessel function computation.

Functions§

airy: Airy function Ai(z).
airy_raw: Airy function Ai(z) with precision status.
airy_scaled: Scaled Airy function: exp(ζ) · Ai(z), where ζ = (2/3) z√z.
airyprime: Derivative of the Airy function, Ai’(z).
airyprime_raw: Derivative of the Airy function Ai’(z) with precision status.
airyprime_scaled: Scaled derivative of the Airy function: exp(ζ) · Ai'(z), where ζ = (2/3) z√z.
besseli: Modified Bessel function of the first kind, I_ν(z).
besseli_scaled: Scaled modified Bessel function of the first kind: exp(-|Re(z)|) · I_ν(z).
besseli_seq: Compute I_{ν+j}(z) for j = 0, 1, …, n−1 in a single call.
besselj: Bessel function of the first kind, J_ν(z).
besselj_scaled: Scaled Bessel function of the first kind: exp(-|Im(z)|) · J_ν(z).
besselj_seq: Compute J_{ν+j}(z) for j = 0, 1, …, n−1 in a single call.
besselk: Modified Bessel function of the second kind, K_ν(z).
besselk_scaled: Scaled modified Bessel function of the second kind: exp(z) · K_ν(z).
besselk_seq: Compute K_{ν+j}(z) for j = 0, 1, …, n−1 in a single call.
bessely: Bessel function of the second kind, Y_ν(z).
bessely_scaled: Scaled Bessel function of the second kind: exp(-|Im(z)|) · Y_ν(z).
bessely_seq: Compute Y_{ν+j}(z) for j = 0, 1, …, n−1 in a single call.
biry: Airy function of the second kind, Bi(z).
biry_raw: Airy function of the second kind Bi(z) with precision status.
biry_scaled: Scaled Airy function of the second kind: exp(-|Re(ζ)|) · Bi(z), where ζ = (2/3) z√z.
biryprime: Derivative of the Airy function of the second kind, Bi’(z).
biryprime_raw: Derivative of the Airy function of the second kind Bi’(z) with precision status.
biryprime_scaled: Scaled derivative of the Airy function of the second kind: exp(-|Re(ζ)|) · Bi'(z), where ζ = (2/3) z√z.
hankel1: Hankel function of the first kind, H_ν^(1)(z).
hankel2: Hankel function of the second kind, H_ν^(2)(z).
hankel1_scaled: Scaled Hankel function of the first kind: exp(-iz) · H_ν^(1)(z).
hankel1_seq: Compute H_{ν+j}^(1)(z) for j = 0, 1, …, n−1 in a single call.
hankel2_scaled: Scaled Hankel function of the second kind: exp(iz) · H_ν^(2)(z).
hankel2_seq: Compute H_{ν+j}^(2)(z) for j = 0, 1, …, n−1 in a single call.