diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 05f4788..b0d446f 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -22,6 +22,7 @@ jobs:
   test:
     name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ github.event_name }}
     runs-on: ${{ matrix.os }}
+    timeout-minutes: 30
     strategy:
       fail-fast: false
       matrix:
diff --git a/src/distrs/hyper.jl b/src/distrs/hyper.jl
index e7cf1c7..dbabedf 100644
--- a/src/distrs/hyper.jl
+++ b/src/distrs/hyper.jl
@@ -1,32 +1,86 @@
 # functions related to hyper-geometric distribution
 
-# Rmath implementations
-function hyperpdf(ms::Real, mf::Real, n::Real, x::Real)
-    T = float(Base.promote_typeof(ms, mf, n, x))
-    return convert(T, Rmath.dhyper(x, ms, mf, n, false))
+hyperpdf(ms::Real, mf::Real, n::Real, x::Real) = hyperpdf(promote(ms,mf,n,x)...)
+    
+function hyperpdf(ms::T, mf::T, n::T, x::T) where T<:Real
+    M, N = ms, ms+mf
+    fxnprev = fMNprev = one(T)
+    fxn = 1 - 2x/n
+    fMN = 1 - 2M/N
+    eiprev = ones(T, Int(n+1))
+    k = 0:n
+    ei = 1 .- k .*(2/n)
+    square(x) = x^2
+    res = fma(fxn, fMN/sum(square, ei), 1/(n+1))
+    ei = collect(ei)
+    for i in 1:n-1
+        fxnnext = ((2i + 1) * (n - 2x) * fxn - i * (n + i + 1) * fxnprev) / ((i + 1) * (n - i))
+        fxn, fxnprev = fxnnext, fxn
+        fMNnext = ((2i + 1) * (N - 2M) * fMN - i * (N + i + 1) * fMNprev) / ((i + 1) * (N - i))
+        fMN, fMNprev = fMNnext, fMN
+        @. eiprev = ((2i + 1) * (n - 2k) * ei - i * (n + i + 1) * eiprev) / ((i + 1) * (n - i))
+        ei, eiprev = eiprev, ei
+        correction =  fxn*fMN/sum(square, ei)
+        abs(correction) < res*eps(T) && break
+        res += correction
+    end
+    res
 end
-function hyperlogpdf(ms::Real, mf::Real, n::Real, x::Real)
-    T = float(Base.promote_typeof(ms, mf, n, x))
-    return convert(T, Rmath.dhyper(x, ms, mf, n, true))
+
+hyperlogpdf(ms::Real, mf::Real, n::Real, x::Real) = hyperlogpdf(promote(ms,mf,n,x)...)
+function hyperlogpdf(ms::T, mf::T, n::T, x::T) where T<:Real
+    return log1p(ms+mf) - log1p(ms) - log1p(mf) + 
+    (logbeta(n + 1, ms + mf - n + 1) - logbeta(x + 1, ms - x + 1) - 
+    logbeta(n - x + 1, mf - n + x + 1))
 end
 
-function hypercdf(ms::Real, mf::Real, n::Real, x::Real)
-    T = float(Base.promote_typeof(ms, mf, n, x))
-    return convert(T, Rmath.phyper(x, ms, mf, n, true, false))
+_hypercdf(ms::Real, mf::Real, n::Real, x::Real, invert::Bool) = _hypercdf(promote(ms, mf, n, x)..., invert::Bool)
+function _hypercdf(ms::T, mf::T, n::T, x::T, invert::Bool) where T<:Real
+    N = ms+mf
+    mode = fld((n + 1) * (ms + 1), N + 2)
+    local result
+    if x < mode
+        result = diff = hyperpdf(ms, mf, n, x)
+        lower_limit = max(0, n+ms-N)
+        while x != lower_limit && diff > (invert ? result : 1)*eps(T)
+            diff *= x * (N + x - n - ms) / ((1+n-x)*(1+ms -x))
+            result += diff
+            x -= 1
+        end
+    else
+        invert = !invert
+        x += 1
+        result = diff = hyperpdf(ms, mf, n, x)
+        upper_limit = min(ms, n)
+        while x != upper_limit && diff > (invert ? 1 : result)*eps(T)
+            diff *= (n-x) * (ms-x) / ((1+x) * (N+x+1-n-ms))
+            result += diff
+            x += 1
+        end
+    end
+    result, invert
 end
-function hyperccdf(ms::Real, mf::Real, n::Real, x::Real)
-    T = float(Base.promote_typeof(ms, mf, n, x))
-    return convert(T, Rmath.phyper(x, ms, mf, n, false, false))
+
+function hypercdf(ms, mf, n, x)
+    result, invert = _hypercdf(ms, mf, n, x, false)
+    invert ? 1-result : result
 end
-function hyperlogcdf(ms::Real, mf::Real, n::Real, x::Real)
-    T = float(Base.promote_typeof(ms, mf, n, x))
-    return convert(T, Rmath.phyper(x, ms, mf, n, true, true))
+
+function hyperccdf(ms, mf, n, x)
+    result, invert = _hypercdf(ms, mf, n, x, true)
+    invert ? 1-result : result
 end
-function hyperlogccdf(ms::Real, mf::Real, n::Real, x::Real)
-    T = float(Base.promote_typeof(ms, mf, n, x))
-    return convert(T, Rmath.phyper(x, ms, mf, n, false, true))
+
+function hyperlogcdf(ms, mf, n, x)
+    result, invert = _hypercdf(ms, mf, n, x, false)
+    invert ? log1p(-result) : log(result)
+end
+function hyperlogccdf(ms, mf, n, x)
+    result, invert = _hypercdf(ms, mf, n, x, true)
+    invert ? log1p(-result) : log(result)
 end
 
+# Rmath implementations
 function hyperinvcdf(ms::Real, mf::Real, n::Real, q::Real)
     T = float(Base.promote_typeof(ms, mf, n, q))
     return convert(T, Rmath.qhyper(q, ms, mf, n, true, false))
diff --git a/test/rmath.jl b/test/rmath.jl
index 8d5bd60..eb19422 100644
--- a/test/rmath.jl
+++ b/test/rmath.jl
@@ -8,6 +8,16 @@ function check_rmath(statsfun, rmathfun, params, a, isprob, rtol)
     v = @inferred(statsfun(params..., a))
     rv = @inferred(rmathfun(a, params...))
     @test v isa float(Base.promote_typeof(params..., a))
+    @test rv isa float(Base.promote_typeof(params..., a))
+    # vbig == bigfloat version if available
+    vbig = try
+        statsfun(params..., big(a))
+    catch
+        rv
+    end
+    # loose accuracy test to confirm shape is right
+    @test v ≈ rv nans = true
+    # tight accuracy test to confirm implimentation is correct numerically
     return if isprob
         @test v ≈ oftype(v, rv) rtol = rtol nans = true
     else
@@ -288,6 +298,7 @@ end
             ((2, 3, 4), 0.0f0:4.0f0),
             ((2, 3, 4), Float16(0):Float16(4)),
             ((2, 3, 4), (0 // 1):(4 // 1)),
+            ((40, 60, 40), 0.0:50.0),
         ]
     )