`zeta(s)` and `eta(s)` return `NaN` or `Inf` for every real s < `-170.62`

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

Maintainers usually reply within 1 day

Nobody has claimed this yet.

Assessment

Difficulty
4/5
Estimated time
3-5 days
Newbie friendliness
48/100
Issue type
Bug
Clarity
Mostly clear
Activity status
Active
Tech stack
julia
Domain
backend

Research direction

Start in src/gamma.jl at _zeta(s::ComplexOrReal{Float64}) and inspect both the real-axis branch and the @pg_horner expansion. Reproduce the reported zeta values in Julia, then verify the relevant behavior around the gamma overflow threshold and negative even integers. Done means real inputs below -170.62 no longer incorrectly produce NaN or Inf, including the reported trivial zeros and finite values.

Written by the indexing model from the issue text.

Description

For example:

julia> using SpecialFunctions

julia> zeta(-170.0)
-0.0

julia> zeta(-172.0)   # true value 0 (trivial zero): https://www.wolframalpha.com/input?i=zeta%28-172%29
NaN

julia> zeta(-180.0)   # true value 0 (trivial zero): https://www.wolframalpha.com/input?i=zeta%28-180%29
NaN

julia> zeta(-200.0)   # true value 0 (trivial zero): https://www.wolframalpha.com/input?i=zeta%28-200%29
NaN

julia> zeta(-250.0)   # true value 0 (trivial zero): https://www.wolframalpha.com/input?i=zeta%28-250%29
NaN

julia> zeta(-175.0)   # true value 7.4526910304300896e177: https://www.wolframalpha.com/input?i=zeta%28-175%29
Inf

julia> zeta(-171.0)   # true value 1.2819489863517886e172: https://www.wolframalpha.com/input?i=zeta%28-171%29
Inf

Possible Cause (as considered by an AI)

src/gamma.jl, _zeta(s::ComplexOrReal{Float64}):

    elseif real(s) < 0.5
        absim = abs(imag(s))
        if abs(real(s)) + absim < 1e-3 # Taylor series for small |s|
            ...
        end
        if absim > 12 # amplitude of sinpi(s/2) ≈ exp(imag(s)*π/2)
            # avoid overflow/underflow (issue #128)
            lg = loggamma(1 - s)
            ...
        else
            return zeta(1 - s) * gamma(1 - s) * sinpi(s*0.5) * twoπ^s * invπ
        end

The else branch is the one from issue #128 that "avoids overflow" only when absim > 12.
On the real axis gamma(1 - s) is Inf for all s < -170.62, sinpi(s*0.5) is ∓0.0 at negative even integers, and twoπ^s is a tiny denormal-range factor that would have brought the product back into range had the multiplication been done in log space.

For the second regime, the @pg_horner macro (same file) expands to a Horner scheme whose coefficients are m*(m+1)*… products; for m ≈ 5.7e19 the third such product exceeds floatmax.

Dominant language
Julia
Stars
384
Forks
115
Avg merge
1d 13h
Merged PRs (30d)
5

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.