Function gamma_func 

pub fn gamma_func(x: f64) -> f64

Gamma function Γ(x) for real x

Adapted from the C++ libsparseir gamma_func (SpM-lab/libsparseir, backend/cxx/src/specfuncs.cpp at commit 4bc58ea), which follows gamma(::Float64) of Bessels.jl v0.2.8 (src/gamma.jl), itself adapted from the Cephes Mathematical Library by Stephen L. Moshier. For x > 0 it uses Stirling’s series above 11.5 and otherwise a rational approximation on [2, 3) reached with Γ(x + 1) = xΓ(x). It deviates from the C++ source in three places: the Stirling prefactor is sqrt(2π) as in Bessels.jl (the C++ code uses sqrt(π/2), which halves every result above 11.5); sin(πx) uses exact argument reduction, so that the poles are detected and accuracy holds next to them; and no input throws or panics.

A negative non-integer x < -1 uses the reflection formula Γ(x) = π / (sin(πx) Γ(1 - x)), evaluated as π / (sin(πx) |x| Γ(|x|)) so that 1 - x is never rounded. For -1 < x < 0 the recurrence Γ(x) = Γ(1 + x) / x is used instead, because |x| sin(πx) underflows for tiny |x|.

Wherever Γ(x) is a normal f64, the relative error is a few ulp.

Special values follow C99 tgamma:

• x = ±0 (pole): ±∞, the sign of zero selecting the side of the pole.
• x a negative integer (a pole without a signed limit) or x = -∞: NaN.
• x = +∞, and x > 171.62... where Γ(x) exceeds f64::MAX: +∞.
• Non-integer x < -184, where |Γ(x)| is below the smallest subnormal: ±0 carrying the sign of Γ(x), which is the sign of sin(πx).
• x = NaN: NaN.