diff --git a/Project.toml b/Project.toml
index 769ec33..994508e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -15,8 +15,9 @@ SpecialFunctions = "0.10, 1.0, 2"
 julia = "1.3"
 
 [extras]
+FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Random"]
+test = ["FiniteDifferences", "Test", "Random"]
diff --git a/test/runtests.jl b/test/runtests.jl
index 072c962..ccb9a65 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,14 +1,12 @@
 using DiffRules
+using FiniteDifferences
 using Test
 
 import SpecialFunctions, NaNMath, LogExpFunctions
 import Random
 Random.seed!(1)
 
-function finitediff(f, x)
-    ϵ = cbrt(eps(typeof(x))) * max(one(typeof(x)), abs(x))
-    return (f(x + ϵ) - f(x - ϵ)) / (ϵ + ϵ)
-end
+const finitediff = central_fdm(5, 1)
 
 @testset "DiffRules" begin
 @testset "check rules" begin
@@ -20,22 +18,25 @@ for (M, f, arity) in DiffRules.diffrules(; filter_modules=nothing)
     if arity == 1
         @test DiffRules.hasdiffrule(M, f, 1)
         deriv = DiffRules.diffrule(M, f, :goo)
-        modifier = if f in (:asec, :acsc, :asecd, :acscd, :acosh, :acoth)
-            1.0
-        elseif f === :log1mexp
-            -1.0
-        elseif f === :log2mexp
-            -0.5
-        else
-            0.0
-        end
         @eval begin
             let
-                goo = rand() + $modifier
-                @test isapprox($deriv, finitediff($M.$f, goo), rtol=0.05)
+                goo = if $(f in (:asec, :acsc, :asecd, :acscd, :acosh, :acoth))
+                    # avoid singularities with finite differencing
+                    rand() + 1.5
+                elseif $(f === :log)
+                    # avoid singularities with finite differencing
+                    rand() + 0.5
+                elseif $(f === :log1mexp)
+                    rand() - 1.0
+                elseif $(f in (:log2mexp, :asin, :acos, :erfinv))
+                    rand() - 0.5
+                else
+                    rand()
+                end
+                @test $deriv ≈ finitediff($M.$f, goo) rtol=1e-9 atol=1e-9
                 # test for 2pi functions
-                if "mod2pi" == string($M.$f)
-                    goo = 4pi + $modifier
+                if $(f === :mod2pi)
+                    goo = 4 * pi
                     @test NaN === $deriv
                 end
             end
@@ -45,17 +46,31 @@ for (M, f, arity) in DiffRules.diffrules(; filter_modules=nothing)
         derivs = DiffRules.diffrule(M, f, :foo, :bar)
         @eval begin
             let
-                if "mod" == string($M.$f)
-                    foo, bar = rand() + 13, rand() + 5 # make sure x/y is not integer
+                foo, bar = if $(f === :mod)
+                    rand() + 13, rand() + 5 # make sure x/y is not integer
+                elseif $(f === :polygamma)
+                    rand(1:10), rand() # only supports integers as first arguments
+                elseif $(f === :bessely)
+                    # avoid singularities with finite differencing
+                    rand(), rand() + 0.5
+                elseif $(f === :log)
+                    # avoid singularities with finite differencing
+                    rand() + 1.5, rand()
+                elseif $(f in (:^, :pow))
+                    # avoid singularities with finite differencing
+                    rand() + 0.5, rand()
+                elseif $(f === :logbeta)
+                    # avoid singularities with finite differencing
+                    rand() + 0.5, rand() + 0.5
                 else
-                    foo, bar = rand(1:10), rand()
+                    rand(), rand()
                 end
                 dx, dy = $(derivs[1]), $(derivs[2])
-                if !(isnan(dx))
-                    @test isapprox(dx, finitediff(z -> $M.$f(z, bar), float(foo)), rtol=0.05)
+                if !isnan(dx)
+                    @test dx ≈ finitediff(z -> $M.$f(z, bar), float(foo)) rtol=1e-9 atol=1e-9
                 end
-                if !(isnan(dy))
-                    @test isapprox(dy, finitediff(z -> $M.$f(foo, z), bar), rtol=0.05)
+                if !isnan(dy)
+                    @test dy ≈ finitediff(z -> $M.$f(foo, z), bar) rtol=1e-9 atol=1e-9
                 end
             end
         end
@@ -89,7 +104,7 @@ for xtype in [:Float64, :BigFloat, :Int64]
                 x = $xtype(rand(1 : 10))
                 y = $mode
                 dx, dy = $(derivs[1]), $(derivs[2])
-                @test isapprox(dx, finitediff(z -> rem2pi(z, y), float(x)), rtol=0.05)
+                @test dx ≈ finitediff(z -> rem2pi(z, y), float(x)) rtol=1e-9 atol=1e-9
                 @test isnan(dy)
             end
         end
@@ -105,7 +120,7 @@ for xtype in [:Float64, :BigFloat]
                 x = rand($xtype)
                 y = $ytype(rand(1 : 10))
                 dx, dy = $(derivs[1]), $(derivs[2])
-                @test isapprox(dx, finitediff(z -> ldexp(z, y), x), rtol=0.05)
+                @test dx ≈ finitediff(z -> ldexp(z, y), x) rtol=1e-9 atol=1e-9
                 @test isnan(dy)
             end
         end