fix overflow/underflow bug in zeta(z) (#129)

Author: stevengj

src/gamma.jl

@@ -402,21 +402,32 @@ function zeta(s::ComplexOrReal{Float64})
         end
         return NaN*zero(s) # NaN + NaN*im
     elseif real(s) < 0.5
-        if abs(real(s)) + abs(imag(s)) < 1e-3 # Taylor series for small |s|
+        absim = abs(imag(s))
+        if abs(real(s)) + absim < 1e-3 # Taylor series for small |s|
             return @evalpoly(s, -0.5,
                              -0.918938533204672741780329736405617639861,
                              -1.0031782279542924256050500133649802190,
                              -1.00078519447704240796017680222772921424,
                              -0.9998792995005711649578008136558752359121)
         end
-        return zeta(1 - s) * gamma(1 - s) * sinpi(s*0.5) * (2π)^s / π
+        if absim > 12 # amplitude of sinpi(s/2) ≈ exp(imag(s)*π/2)
+            # avoid overflow/underflow (issue #128)
+            lg = lgamma(1 - s)
+            ln2pi = 1.83787706640934548356 # log(2pi) to double precision
+            rehalf = real(s)*0.5
+            return zeta(1 - s) * exp(lg + absim*(pi/2) + s*ln2pi) * (0.5/π) *
+                Complex(sinpi(rehalf), copysign(cospi(rehalf), imag(s)))
+        else
+            return zeta(1 - s) * gamma(1 - s) * sinpi(s*0.5) * (2π)^s / π
+        end
     end
 
     m = s - 1
 
     # shift using recurrence formula:
     #   n is a semi-empirical cutoff for the Stirling series, based
     #   on the error term ~ (|m|/n)^18 / n^real(m)
+    # FIXME: slow for large imag(s) and small real(s)
     n = ceil(Int,6 + 0.7*abs(imag(s-1))^inv(1 + real(m)*0.05))
     ζ = one(s)
     for ν = 2:n

test/runtests.jl

@@ -59,7 +59,7 @@ relerrc(z, x) = max(relerr(real(z),real(x)), relerr(imag(z),imag(x)))
     @test erfcx(1.8e88) ≈ erfcx(big(1.8e88)) rtol=4*eps()
 
     @test_throws MethodError erf(big(1.0)*im)
-    @test_throws MethodError erfi(big(1.0))    
+    @test_throws MethodError erfi(big(1.0))
 end
 
 @testset "sine and cosine integrals" begin
@@ -475,11 +475,17 @@ end
     @test 1e-12 > relerrc(zeta(1e-3+1e-3im), -0.5009189365276307665899456585255302329444338284981610162-0.0009209468912269622649423786878087494828441941303691216750im)
     @test 1e-13 > relerrc(zeta(1e-4 + 2e-4im), -0.5000918637469642920007659467492165281457662206388959645-0.0001838278317660822408234942825686513084009527096442173056im)
 
-    # Issue #7169: (TODO: better accuracy should be possible?)
-    @test 1e-9  > relerrc(zeta(0 + 99.69im), 4.67192766128949471267133846066040655597942700322077493021802+3.89448062985266025394674304029984849370377607524207984092848im)
+    # Issue #7169:
+    @test 1e-13  > relerrc(zeta(0 + 99.69im), 4.67192766128949471267133846066040655597942700322077493021802+3.89448062985266025394674304029984849370377607524207984092848im)
     @test 1e-12 > relerrc(zeta(3 + 99.69im), 1.09996958148566565003471336713642736202442134876588828500-0.00948220959478852115901654819402390826992494044787958181148im)
-    @test 1e-9  > relerrc(zeta(-3 + 99.69im), 10332.6267578711852982128675093428012860119184786399673520976+13212.8641740351391796168658602382583730208014957452167440726im)
+    @test 1e-13  > relerrc(zeta(-3 + 99.69im), 10332.6267578711852982128675093428012860119184786399673520976+13212.8641740351391796168658602382583730208014957452167440726im)
     @test 1e-13 > relerrc(zeta(2 + 99.69im, 1.3), 0.41617652544777996034143623540420694985469543821307918291931-0.74199610821536326325073784018327392143031681111201859489991im)
+
+    # issue #128
+    @test 1e-13 > relerrc(zeta(.4 + 453.0im), 5.595631794716693 - 4.994584420588448im)
+    @test 1e-10 > relerrc(zeta(.4 + 4053.0im), -0.1248993234383550+0.9195498409364987im)
+    @test 1e-13 > relerrc(zeta(.4 + 12.01im), 1.0233184799021265846512208845-0.8008078492939259287905322251im)
+    @test zeta(.4 + 12.01im) == conj(zeta(.4 - 12.01im))
 end
 
 @testset "vectorization of 2-arg functions" begin