digamma and trigamma accuracy improvements (#451)

* digamma accuracy improvement

* changed `trigamma` and use `sinpi`

* `sincospi` replaced by `sinpi, cospi`

* use `tanpi` if available

* added test case for `cotderiv(0, ...)`

* code coverage

Author: KlausC

src/gamma.jl

@@ -28,15 +28,16 @@ function _digamma(z::ComplexOrReal{Float64})
     # argument," Computer Phys. Commun.  vol. 4, pp. 221–226 (1972).
     x = real(z)
     if x <= 0 # reflection formula
-        ψ = -π * cot(π*z)
+        ψ = -π * _cotpi(z)
         z = 1 - z
         x = real(z)
     else
         ψ = zero(z)
     end
-    if x < 7
+    X = 8
+    if x < X
         # shift using recurrence formula
-        n = 7 - floor(Int,x)
+        n = X - floor(Int,x)
         for ν = 1:n-1
             ψ -= inv(z + ν)
         end
@@ -56,6 +57,13 @@ function _digamma(x::BigFloat)
     return z
 end
 
+"""
+    _cotpi(x) = cot(π * x)
+
+Accurate for integer arguments
+"""
+_cotpi(x) = @static VERSION >= v"1.10.0-DEV.525" ? inv(tanpi(x)) : cospi(x) / sinpi(x)
+
 """
     trigamma(x)
 
@@ -67,12 +75,13 @@ function _trigamma(z::ComplexOrReal{Float64})
     # via the derivative of the Kölbig digamma formulation
     x = real(z)
     if x <= 0 # reflection formula
-        return (π * csc(π*z))^2 - trigamma(1 - z)
+        return (π / sinpi(z))^2 - trigamma(1 - z)
     end
     ψ = zero(z)
-    if x < 8
+    N = 10
+    if x < N
         # shift using recurrence formula
-        n = 8 - floor(Int,x)
+        n = N - floor(Int,x)
         ψ += inv(z)^2
         for ν = 1:n-1
             ψ += inv(z + ν)^2
@@ -132,12 +141,12 @@ const cotderiv_Q = [cotderiv_q(m) for m in 1:100]
 function cotderiv(m::Integer, z)
     isinf(imag(z)) && return zero(z)
     if m <= 0
-        m == 0 && return π * cot(π*z)
+        m == 0 && return π * _cotpi(z)
         throw(DomainError(m, "`m` must be nonnegative."))
     end
     if m <= length(cotderiv_Q)
         q = cotderiv_Q[m]
-        x = cot(π*z)
+        x = _cotpi(z)
         y = x*x
         s = q[1] + q[2] * y
         t = y

test/gamma.jl

@@ -12,6 +12,9 @@
             @test digamma(convert(elty, 7e-7)) ≈ convert(elty, -1428572.005785942019703646)
             @test digamma(convert(elty, -0.5)) ≈ convert(elty, .03648997397857652055902367)
             @test digamma(convert(elty, -1.1)) ≈ convert(elty,  10.15416395914385769902271)
+            # issue #450
+            @test digamma(convert(elty, 0)) == convert(elty,  -Inf)
+            @test digamma(convert(elty, -1)) == convert(elty,  -Inf)
 
             @test digamma(convert(elty, 0.1)) ≈ convert(elty, -10.42375494041108)
             @test digamma(convert(elty, 1/2)) ≈ convert(elty, -γ - log(4))
@@ -35,6 +38,8 @@
             @test trigamma(convert(elty, 4)) ≈ convert(elty, π^2/6 - 49/36)
             @test trigamma(convert(elty, 5)) ≈ convert(elty, π^2/6 - 205/144)
             @test trigamma(convert(elty, 10)) ≈ convert(elty, π^2/6 - 9778141/6350400)
+            @test trigamma(convert(elty, 0)) == convert(elty, Inf)
+            @test trigamma(convert(elty, -1)) == convert(elty, Inf)
         end
     end
 
@@ -268,6 +273,10 @@ end
     @test zeta(-6+13im) ≅ conj(zeta(-6-13im)) ≅ 133.4764526350263089084083707864441932569167866714712267139316498-54.15465727586582149098585229287107039070546786014930791081909684im
     @test 1e-12 > relerr(zeta(-2+13im, 3), 2.3621038290867825837364823054-3.9497600485207119519185591345im)
     @test 1e-12 > relerr(zeta(-2-13im, 3), 2.3621038290867825837364823054+3.9497600485207119519185591345im)
+
+    # issue #450
+    @test SpecialFunctions.cotderiv(0, 2.0) == Inf
+    @test_throws DomainError SpecialFunctions.cotderiv(-1, 2.0)
 end
 
 @testset "logabsbinomial" begin