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, 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

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

Functions with a branch cut along the negative real axis (J and I of non-integer order, Y, K, H, and the factor exp(ζ) of the scaled Ai and Ai′) take their values there from above the cut (arg z = π), as in Amos: an imaginary part of −0 is treated like +0.

§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 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) (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)). For half-integer orders, cos(νπ) = 0, so only one function is computed: J_−ν(z) = −sin(νπ) Y_ν(z) and Y_−ν(z) = sin(νπ) J_ν(z).

§Accuracy on the axes

Amos’s algorithm is accurate normwise: the error is small relative to |f(z)|. On the real and imaginary axes some values combine two real functions of very different size, e.g. H⁽¹⁾_ν(x) = J_ν(x) + iY_ν(x) with J ≪ Y for x < ν, and Amos’s small component can then be off entirely. The crate recomputes such a component from the two real functions:

H⁽¹⁾, H⁽²⁾ on the positive and negative real axis
Y on the negative real axis (also scaled)
K on the imaginary axis
H⁽²⁾ on the positive and H⁽¹⁾ on the negative imaginary axis (also scaled)

For Y and H⁽²⁾ on the negative real axis this includes parts for |ν| < x: Re Y_ν(−x) = cos(νπ) Y_ν(x) for orders close to a half-integer, and, for |ν| < 18 next to their zeros, parts that are Y_|ν| alone: Re Y_ν(−x) for ν ≥ 0, and parts at or close to integer and half-integer orders, e.g. Im H⁽²⁾_n(−x) = (−1)ⁿ⁺¹ Y_n(x). J and Y are exactly real on the positive real axis.

A remaining limitation: in the other functions, a small component that is a nearly vanishing cos(νπ) or sin(νπ) times one real function, for orders close to an integer or a half-integer, is still only accurate relative to |f(z)|; its relative error grows roughly like ε/|cos(νπ)| (or ε/|sin(νπ)|). An example is Re J_ν(−x) = cos(νπ) J_ν(x).

Only values exactly on an axis are recomputed. Next to an axis (e.g. Im z = 10⁻¹⁰·|z|) the same cancellation remains, less deep the farther z is from the axis, and the parts are accurate only relative to |f(z)|, as in Amos.

For f32, J_ν(x) on the real axis with ν ≳ 10 and 35 ≲ x ≲ 45 errs by up to about 7e-3·|H_ν(x)|, as in Amos’s single-precision CBESJ.

§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 _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.

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;
pub use types::SeqInfo;

Modules§

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

Functions§

airy: Airy function Ai(z).
airy_raw: Airy function Ai(z) with accuracy 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 accuracy 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.
besseli_seq_into: Compute I_ν+j(z) for j = 0, 1, …, n−1 into out, where n = out.len().
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.
besselj_seq_into: Compute J_ν+j(z) for j = 0, 1, …, n−1 into out, where n = out.len().
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.
besselk_seq_into: Compute K_ν+j(z) for j = 0, 1, …, n−1 into out, where n = out.len().
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.
bessely_seq_into: Compute Y_ν+j(z) for j = 0, 1, …, n−1 into out, where n = out.len().
biry: Airy function of the second kind, Bi(z).
biry_raw: Airy function of the second kind Bi(z) with accuracy 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 accuracy 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_ν⁽¹⁾(z).
hankel2: Hankel function of the second kind, H_ν⁽²⁾(z).
hankel1_scaled: Scaled Hankel function of the first kind, exp(−iz) · H_ν⁽¹⁾(z).
hankel1_seq: Compute H_ν+j⁽¹⁾(z) for j = 0, 1, …, n−1 in a single call.
hankel1_seq_into: Compute H_ν+j⁽¹⁾(z) for j = 0, 1, …, n−1 into out, where n = out.len().
hankel2_scaled: Scaled Hankel function of the second kind, exp(iz) · H_ν⁽²⁾(z).
hankel2_seq: Compute H_ν+j⁽²⁾(z) for j = 0, 1, …, n−1 in a single call.
hankel2_seq_into: Compute H_ν+j⁽²⁾(z) for j = 0, 1, …, n−1 into out, where n = out.len().