gamma(a, x) and loggamma(a, x) throw for every real shape a > 51

Open
  • Is something wrong?
#567 0 comments 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

Assessment

Difficulty
4/5
Estimated time
3-5 days
Newbie friendliness
25/100
Issue type
Bug
Clarity
Clearly specified
Activity status
Active
Domain
backend

Research direction

The failure comes from _gamma and _loggamma in src/gamma_inc.jl, which route through _expint in src/expint.jl, where the order limit |ν| > 50 is enforced. Start by reproducing loggamma(52, 1e5) and gamma(52, 10.0), then read the existing gamma_inc and gamma_inc_cf paths the issue proposes reusing. Done means both functions return finite results for a > 51 that match a high-precision reference, with existing tests still passing. The author says they plan to open a PR, so check for one before starting.

Written by the indexing model from the issue text.

Description

gamma(a, x) and loggamma(a, x) throw for every real shape a > 51, at every x:

using SpecialFunctions

loggamma(51, 1e5)    # fine: -99424.35322663146
loggamma(52, 1e5)    # ArgumentError
gamma(52, 10.0)      # ArgumentError
loggamma(100, 1e5)   # ArgumentError; log Γ(100, 1e5) ≈ -98860.2193884875

On SpecialFunctions 2.9.1 (master at d8ffcc1b456171884112a157b7be1bc4707f77be behaves the same), Julia 1.10.12:

loggamma(51, 100000.0) => -99424.35322663146
loggamma(52, 100000.0) => ArgumentError: Unsupported order |ν| > 50 off the positive real axis
gamma(52, 10.0) => ArgumentError: Unsupported order |ν| > 50 off the positive real axis
loggamma(100, 100000.0) => ArgumentError: Unsupported order |ν| > 50 off the positive real axis

The boundary is exactly a = 51: loggamma(51.0, 10.0) works, loggamma(51.0001, 10.0) throws.

Cause

Both functions are evaluated through the exponential integral for every argument — gamma(a, x) = x^a E_{1-a}(x) in _gamma, and a log x + log(expintx(1 - a, x)) - x in _loggamma (src/gamma_inc.jl). _expint refuses orders |ν| > 50 that are not positive real (src/expint.jl, "Unsupported order |ν| > 50 off the positive real axis"), and ν = 1 - a is a negative real below -50 as soon as a > 51. The docstrings promise "arbitrary real or complex a and x".

The route is also the wrong one for large shapes in principle, independent of the order limit: Γ(a, x) overflows Float64 for small x once a > 171, and underflows for x ≫ a, long before log Γ(a, x) leaves the representable range — Γ(200, 700) ≈ 2.06e262 is representable, but Γ(200) = Inf and Q(200, 700) ≈ 0.0 are not; log Γ(170, 1210) ≈ -10.22 while Q(170, 1210) is below floatmin. expintx cannot carry those magnitudes either.

Where it shows up

StatsFuns.gammalogccdf falls back to loggamma(k, x) - loggamma(k) once the survival probability underflows (src/distrs/gamma.jl), so Distributions.logccdf(Gamma(100.0), 1e5) throws this error instead of returning -99219.35359385707 (the exact integer-shape identity log Q(k, x) = -x + log(sum(x^j / j! for j in 0:k-1)) at 512 bits). The same failure reaches logccdf of Erlang, Chisq (ν > 102) and Chi, and logcdf of InverseGamma and Poisson (x ≥ 51 when the probability underflows) — every log tail of a gamma-family distribution with shape above 51.

Proposed fix

For real a > 0 and finite x ≥ 0, compute both functions from the incomplete gamma ratios that gamma_inc already provides for exactly this region: Γ(a, x) = Q(a, x) Γ(a) and log Γ(a, x) = log Γ(a) + log Q(a, x) (with log1p(-P) when Q is the complement), and when Q(a, x) underflows, log Γ(a, x) = a log x - x + log(CF) with the Legendre continued fraction that gamma_inc_cf already evaluates (DLMF 8.9.2) — it is O(1/x) there, so it never underflows. gamma(a, x) takes Γ(a) Q(a, x) while both are representable and exp(loggamma(a, x)) otherwise. Every other argument (non-positive or complex a, negative or complex x) keeps the expint route. I have this implemented and tested against a 256-bit quadrature of Γ(a, x) (which agrees with MPFR's mpfr_gamma_inc to ~1e-50 where that is fast enough to run) over shapes from 0.1 to 1e8, including both overflow regimes; the relative error stays below 3e-15 everywhere, and for a ≤ 51 it is never worse than the current expint route. I will open a PR against master.

expint(ν, z) itself keeps its order limit; the same identity (E_ν(z) = z^{ν-1} Γ(1 - ν, z)) would lift it for real ν < -50 on the positive real axis, but that is a separate change.

Dominant language
Julia
Stars
384
Forks
115
Avg merge
2d 9h
Merged PRs (30d)
2

Getting set up

This project ships no dev container, Dockerfile or contributing guide, so setting up is up to you: start from its README, and see our first-contribution guide for the general steps.

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

More from JuliaMath/SpecialFunctions.jl

All issues in JuliaMath/SpecialFunctions.jl

Similar issues

More Julia issues

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.