From 5268c7c87d80b16e308e76cd2266b97c4b43652a Mon Sep 17 00:00:00 2001 From: Quytelda Kahja Date: Thu, 20 Mar 2025 12:35:28 -0700 Subject: [PATCH] Fix premature truncation of the recurrence sum in zeta(s,z). Skipping the recurrence sum or truncating early causes the function to break inside a region determined by shape of the cutoff. --- src/gamma.jl | 101 ++++++++++++++++++++++++++------------------------ test/gamma.jl | 5 +++ 2 files changed, 57 insertions(+), 49 deletions(-) diff --git a/src/gamma.jl b/src/gamma.jl index 5004404d..74f3ae56 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -262,6 +262,7 @@ function _zeta(s::T, z::T) where {T<:ComplexOrReal{Float64}} end m = s - 1 + minus_s = -s ζ = zero(T) # Algorithm is just the m-th derivative of digamma formula above, @@ -270,67 +271,69 @@ function _zeta(s::T, z::T) where {T<:ComplexOrReal{Float64}} # Note: we multiply by -(-1)^m m! in polygamma below, so this factor is # pulled out of all of our derivatives. - cutoff = 7 + real(m) + abs(imag(m)) # TODO: this cutoff is too conservative? - if x < cutoff - # shift using recurrence formula - xf = floor(x) - nx = Int(xf) - n = ceil(Int, cutoff - nx) - minus_s = -s - if nx < 0 # x < 0 - # need to use (-z)^(-s) recurrence to be correct for real z < 0 - # [the general form of the recurrence term is (z^2)^(-s/2)] - minus_z = -z - ζ += pow_oftype(ζ, minus_z, minus_s) # ν = 0 term - if xf != z - ζ += pow_oftype(ζ, z - nx, minus_s) - # real(z - nx) > 0, so use correct branch cut - # otherwise, if xf==z, then the definition skips this term - end - # do loop in different order, depending on the sign of s, - # so that we are looping from largest to smallest summands and - # can halt the loop early if possible; see issue #15946 - # FIXME: still slow for small m, large Im(z) - if real(s) > 0 - for ν in -nx-1:-1:1 - ζₒ= ζ - ζ += pow_oftype(ζ, minus_z - ν, minus_s) - ζ == ζₒ && break # prevent long loop for large -x > 0 - end - else - for ν in 1:-nx-1 - ζₒ= ζ - ζ += pow_oftype(ζ, minus_z - ν, minus_s) - ζ == ζₒ && break # prevent long loop for large -x > 0 - end - end - else # x ≥ 0 && z != 0 - ζ += pow_oftype(ζ, z, minus_s) - end - # loop order depends on sign of s, as above + # If x < 0 the series begins with p=-⌊x⌋ terms where real(z+k) < 0. + # Since we're using the recurrence ((z+k)^2)^(s/2), we reflect + # these to (-z-k)^(-s). + nx = floor(Int, x) + p = max(0, -nx) + if p > 0 + # do loop in different order, depending on the sign of s, + # so that we are looping from largest to smallest summands and + # can halt the loop early if possible; see issue #15946 + # FIXME: still slow for small m, large Im(z) + minus_z = -z if real(s) > 0 - for ν in max(1,1-nx):n-1 - ζₒ= ζ - ζ += pow_oftype(ζ, z + ν, minus_s) - ζ == ζₒ && break # prevent long loop for large m + for k in p-1:-1:0 + ζₒ = ζ + ζ += pow_oftype(ζ, minus_z - k, minus_s) + ζ == ζₒ && break # prevent long loop for large -x > 0 end else - for ν in n-1:-1:max(1,1-nx) - ζₒ= ζ - ζ += pow_oftype(ζ, z + ν, minus_s) - ζ == ζₒ && break # prevent long loop for large m + for k in 0:1:p-1 + ζₒ = ζ + ζ += pow_oftype(ζ, minus_z - k, minus_s) + ζ == ζₒ && break # prevent long loop for large -x > 0 end end - z += n + z += p # Shift according to recurrence + end + + # real(z) ≥ 0 is guaranteed here. + # If any term must be dropped, it will be this one. + if z != 0 + ζ += pow_oftype(ζ, z, minus_s) end + # Add all terms where ⌊real(z+k)⌋ < 7 + real(s) + |imag(s)|, using + # at least enough terms to balance any reflected terms added + # previously. The loop order depends on sign of s, as above. + n = max(1+p, ceil(Int, 7 + real(m) + abs(imag(m)) - (nx+p))) + if real(s) > 0 + for k in 1:1:n-1 + ζₒ = ζ + ζ += pow_oftype(ζ, z + k, minus_s) + ζ == ζₒ && break # prevent long loop for large m + end + else + for k in n-1:-1:1 + ζₒ = ζ + ζ += pow_oftype(ζ, z + k, minus_s) + ζ == ζₒ && break # prevent long loop for large m + end + end + z += n + + # Euler-Maclaurin integral and tail sum t = inv(z) w = isa(t, Real) ? conj(oftype(ζ, t))^m : oftype(ζ, t)^m ζ += w * (inv(m) + 0.5*t) t *= t # 1/z^2 - ζ += w*t * @pg_horner(t,m,0.08333333333333333,-0.008333333333333333,0.003968253968253968,-0.004166666666666667,0.007575757575757576,-0.021092796092796094,0.08333333333333333,-0.4432598039215686,3.0539543302701198) - + ζ += w * t * @pg_horner(t, m, 0.08333333333333333, -0.008333333333333333, + 0.003968253968253968, -0.004166666666666667, + 0.007575757575757576, -0.021092796092796094, + 0.08333333333333333, -0.4432598039215686, + 3.0539543302701198) return ζ end diff --git a/test/gamma.jl b/test/gamma.jl index 8c970786..832a1f83 100644 --- a/test/gamma.jl +++ b/test/gamma.jl @@ -277,6 +277,11 @@ end # issue #450 @test SpecialFunctions.cotderiv(0, 2.0) == Inf @test_throws DomainError SpecialFunctions.cotderiv(-1, 2.0) + + # issue #488 + # TODO: Perhaps these bounds can be tightened in the future. + @test 3e-5 > relerr(zeta(-5.75, 0.5), 0.00140748175562420497363476203726333231826481355014969602507003223784179195) + @test 5e-5 > relerr(zeta(-8-1im, 0.5), -0.00345211118818533736386710113396098188185995501107179962865430034343404+0.0109171284162538012544319013865107198958348806377836496833453787991565im) end @testset "logabsbinomial" begin