diff --git a/Project.toml b/Project.toml
index 7ddc510..32a703b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,13 +3,12 @@ uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 version = "0.3.28"
 
 [deps]
+Gamma = "a0844989-3bd2-4988-8bea-c9407ab0941b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-OpenLibm_jll = "05823500-19ac-5b8b-9628-191a04bc5112"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
-SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"
-julia = "1.3"
+Gamma = "1.1"
+julia = "1.10"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/src/HypergeometricFunctions.jl b/src/HypergeometricFunctions.jl
index 557dcd5..376ee2b 100644
--- a/src/HypergeometricFunctions.jl
+++ b/src/HypergeometricFunctions.jl
@@ -1,6 +1,6 @@
 module HypergeometricFunctions
 
-using LinearAlgebra, SpecialFunctions
+using LinearAlgebra, Gamma
 
 export _₁F₁, _₂F₁, _₃F₂, pFq
 
diff --git a/src/specialfunctions.jl b/src/specialfunctions.jl
index 31472c0..e4ccc32 100644
--- a/src/specialfunctions.jl
+++ b/src/specialfunctions.jl
@@ -1,5 +1,3 @@
-using OpenLibm_jll
-
 @inline errcheck(x, y, tol) = isfinite(x) && isfinite(y) && (norm(x-y) > max(norm(x), norm(y))*tol)
 
 kind2string(::Val{0}) = "₀"
@@ -45,6 +43,9 @@ end
 
 ogamma(x::Number) = (isinteger(x) && x<0) ? zero(float(x)) : inv(gamma(x))
 
+unsafe_gamma(z) = gamma(z)
+unsafe_gamma(x::Real) = (x == -Inf || (isinteger(x) && x < 0)) ? oftype(float(x), NaN) : gamma(x)
+
 """
     @clenshaw(x, c...)
 
@@ -154,17 +155,6 @@ speciallogseries(x::Float64) = @clenshaw(5.0x, 1.0087391788544393911192, 1.22047
 speciallogseries(x::ComplexF64) = @evalpoly(x, 1.0000000000000000000000, 5.9999999999999999999966e-01, 4.2857142857142857142869e-01, 3.3333333333333333333347e-01, 2.7272727272727272727292e-01, 2.3076923076923076923072e-01, 1.9999999999999999999996e-01, 1.7647058823529411764702e-01, 1.5789473684210526315786e-01, 1.4285714285714285714283e-01, 1.3043478260869565217384e-01, 1.2000000000000000000000e-01, 1.1111111111111111111109e-01, 1.0344827586206896551722e-01, 9.6774193548387096774217e-02, 9.0909090909090909090938e-02, 8.5714285714285714285696e-02, 8.1081081081081081081064e-02, 7.6923076923076923076907e-02, 7.3170731707317073170688e-02)
 
 
-const libm = OpenLibm_jll.libopenlibm
-
-unsafe_gamma(x::Float64) = ccall((:tgamma , libm), Float64, (Float64, ), x)
-unsafe_gamma(x::Float32) = ccall((:tgammaf, libm), Float32, (Float32, ), x)
-function unsafe_gamma(x::BigFloat)
-    z = BigFloat()
-    ccall((:mpfr_gamma, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, Base.MPFR.ROUNDING_MODE[])
-    return z
-end
-unsafe_gamma(z) = gamma(z)
-
 """
     @lanczosratio(z, ϵ, c₀, c...)
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 4382e4c..6e7ff46 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,8 +1,9 @@
-using HypergeometricFunctions, SpecialFunctions, Test
+using HypergeometricFunctions, Gamma, Test
 import LinearAlgebra: norm
 import HypergeometricFunctions: iswellpoised, isalmostwellpoised, M, U,
                                 pochhammer, _₂F₁general, _₂F₁general2,
-                                pFqdrummond, pFqweniger, pFq2string
+                                pFqdrummond, pFqweniger, pFq2string,
+                                unsafe_gamma
 
 const rtol = 1.0e-3
 const NumberType = Float64
@@ -23,11 +24,26 @@ const NumberType = Float64
     @test pochhammer(3, 1:5) == [3, 12, 60, 360, 2520]
 end
 
+@testset "unsafe gamma" begin
+    @test isnan(unsafe_gamma(-1.0))
+    @test isnan(unsafe_gamma(-1.0f0))
+    @test isnan(unsafe_gamma(BigFloat(-1)))
+    @test isnan(unsafe_gamma(-Inf))
+    @test isnan(unsafe_gamma(BigFloat(-Inf)))
+    @test unsafe_gamma(0.0) == Inf
+    @test unsafe_gamma(BigFloat(0)) == Inf
+    @test unsafe_gamma(-0.0) == -Inf
+    @test signbit(unsafe_gamma(-0.0))
+    @test unsafe_gamma(BigFloat(-0.0)) == -Inf
+    @test signbit(unsafe_gamma(BigFloat(-0.0)))
+    @test unsafe_gamma(1.2+0.4im) == gamma(1.2+0.4im)
+end
+
 @testset "Hypergeometric Functions" begin
     @testset "_₂F₁ vs _₂F₁general vs _₂F₁general2" begin
         e = exp(1.0)
         regression_max_accumulated_error = 2^12*eps() # discovered by running the test
-        for z in (.9rand(Float64, 10), 10rand(ComplexF64, 10))
+        for z in (range(-1.0, 0.875; length = 16), range(-10, 10; length = 10) .+ im*range(-10, 10; length = 10)')
             j = 1
             for (a,b,c) in ((√2/e, 1.3, 1.3), (1.2, √3, 1.2), (-0.4, 0.4, 0.5),
                             (-0.3, 1.3, 0.5), (0.35, 0.85, 0.5), (0.5, 0.5, 1.5),