diff --git a/src/distrs/binom.jl b/src/distrs/binom.jl
index 7d16c48..78c12d7 100644
--- a/src/distrs/binom.jl
+++ b/src/distrs/binom.jl
@@ -30,20 +30,158 @@ for l in ("", "log"), compl in (false, true)
     end
 end
 
-# Rmath implementations
+# Inverse CDF: find smallest k such that binomcdf(n, p, k) >= cprob.
+# Based on VBA critbinomial by Ian Smith.
+# PMF ratios: PMF(k+1)/PMF(k) = (n-k)*p/((k+1)*(1-p))
+#             PMF(k-1)/PMF(k) = k*(1-p)/((n-k+1)*p)
+function _critbinomial(n::Float64, p::Float64, cprob::Float64)
+    q = 1.0 - p
+
+    # Normal approximation for initial guess
+    σ = sqrt(n * p * q)
+    i = clamp(floor(n * p + norminvcdf(min(cprob, 1.0 - 1.0e-15)) * σ + 0.5), 0.0, n)
+
+    # Compute CDF at the guess
+    pr = Float64(binomcdf(n, p, i))
+
+    if pr >= cprob
+        # Search left: find smallest k with CDF(k) >= cprob
+        while i > 0
+            tpr = Float64(binompdf(n, p, i))
+            if pr - tpr < cprob
+                return i  # CDF(i-1) < cprob <= CDF(i)
+            end
+            pr -= tpr
+            i -= 1.0
+        end
+        return 0.0
+    else
+        # Search right
+        while i < n
+            i += 1.0
+            tpr = Float64(binompdf(n, p, i))
+            pr += tpr
+            if pr >= cprob
+                return i
+            end
+        end
+        return n
+    end
+end
+
+# Inverse CCDF: find smallest k such that binomccdf(n, p, k) <= cprob.
+# Based on VBA critcompbinomial by Ian Smith.
+function _critcompbinomial(n::Float64, p::Float64, cprob::Float64)
+    q = 1.0 - p
+
+    # Normal approximation
+    σ = sqrt(n * p * q)
+    i = clamp(floor(n * p - norminvcdf(min(cprob, 1.0 - 1.0e-15)) * σ + 0.5), 0.0, n)
+
+    # Compute CCDF at the guess
+    pr = Float64(binomccdf(n, p, i))
+
+    if pr > cprob
+        # Search right
+        while i < n
+            i += 1.0
+            tpr = Float64(binompdf(n, p, i))
+            pr -= tpr
+            if pr <= cprob
+                return i
+            end
+        end
+        return n
+    else
+        # Search left: find smallest k with CCDF(k) <= cprob
+        while i > 0
+            tpr = Float64(binompdf(n, p, i))
+            if pr + tpr > cprob
+                return i  # CCDF(i) <= cprob < CCDF(i-1)
+            end
+            pr += tpr
+            i -= 1.0
+        end
+        return 0.0
+    end
+end
+
+# Wrapper with edge cases and post-correction (VBA crit_binomial)
+function _binom_invcdf(n::Float64, p::Float64, q::Float64)
+    if q < 0 || q > 1 || n < 0 || p < 0 || p > 1 || isnan(q) || isnan(n) || isnan(p)
+        return NaN
+    elseif q == 0 || p == 0
+        return 0.0
+    elseif p == 1
+        return n
+    end
+
+    i = _critbinomial(n, p, q)
+
+    # Post-correction
+    pr = Float64(binomcdf(n, p, i))
+    if pr >= q
+        while i > 0
+            pr2 = Float64(binomcdf(n, p, i - 1.0))
+            if pr2 < q
+                return i
+            end
+            i -= 1.0
+            pr = pr2
+        end
+        return 0.0
+    else
+        return i + 1.0
+    end
+end
+
+# Wrapper with edge cases and post-correction (VBA comp_crit_binomial)
+function _binom_invccdf(n::Float64, p::Float64, q::Float64)
+    if q < 0 || q > 1 || n < 0 || p < 0 || p > 1 || isnan(q) || isnan(n) || isnan(p)
+        return NaN
+    elseif q == 1 || p == 0
+        return 0.0
+    elseif q == 0 || p == 1
+        return n
+    end
+
+    i = _critcompbinomial(n, p, q)
+
+    # Post-correction
+    pr = Float64(binomccdf(n, p, i))
+    if pr <= q
+        while i > 0
+            pr2 = Float64(binomccdf(n, p, i - 1.0))
+            if pr2 > q
+                return i
+            end
+            i -= 1.0
+            pr = pr2
+        end
+        return 0.0
+    else
+        return i + 1.0
+    end
+end
+
+# Public API
+
 function binominvcdf(n::Real, p::Real, q::Real)
     T = float(Base.promote_typeof(n, p, q))
-    return convert(T, Rmath.qbinom(q, n, p, true, false))
+    return convert(T, _binom_invcdf(Float64(n), Float64(p), Float64(q)))
 end
+
 function binominvccdf(n::Real, p::Real, q::Real)
     T = float(Base.promote_typeof(n, p, q))
-    return convert(T, Rmath.qbinom(q, n, p, false, false))
+    return convert(T, _binom_invccdf(Float64(n), Float64(p), Float64(q)))
 end
+
 function binominvlogcdf(n::Real, p::Real, lq::Real)
     T = float(Base.promote_typeof(n, p, lq))
-    return convert(T, Rmath.qbinom(lq, n, p, true, true))
+    return convert(T, _binom_invcdf(Float64(n), Float64(p), exp(Float64(lq))))
 end
+
 function binominvlogccdf(n::Real, p::Real, lq::Real)
     T = float(Base.promote_typeof(n, p, lq))
-    return convert(T, Rmath.qbinom(lq, n, p, false, true))
+    return convert(T, _binom_invccdf(Float64(n), Float64(p), exp(Float64(lq))))
 end
diff --git a/src/distrs/nbinom.jl b/src/distrs/nbinom.jl
index 580e5db..343d6db 100644
--- a/src/distrs/nbinom.jl
+++ b/src/distrs/nbinom.jl
@@ -66,22 +66,172 @@ function nbinomlogccdf(r::Real, p::Real, k::Real)
     end
 end
 
-# Rmath implementations
-# TODO: implement https://arxiv.org/abs/2001.03953
-# for inverting the incomplete beta function wrt the 2nd argument.
+# Inverse CDF: find smallest k such that nbinomcdf(r, p, k) >= cprob.
+# Based on VBA critnegbinom by Ian Smith.
+# PMF ratio: PMF(k+1)/PMF(k) = (k+r)/(k+1) * (1-p)
+function _critnbinom(r::Float64, p::Float64, cprob::Float64)
+    q = 1.0 - p
+
+    # Normal approximation: mean = r*q/p, var = r*q/p^2
+    μ = r * q / p
+    σ = sqrt(μ / p)
+    i = max(0.0, floor(μ + norminvcdf(min(cprob, 1.0 - 1.0e-15)) * σ + 0.5))
+
+    # Compute CDF at the guess
+    pr = Float64(nbinomcdf(r, p, i))
+
+    if pr >= cprob
+        # Search left
+        while i > 0
+            tpr = Float64(nbinompdf(r, p, i))
+            if pr - tpr < cprob
+                return i
+            end
+            pr -= tpr
+            i -= 1.0
+        end
+        return 0.0
+    else
+        # Search right
+        for _ in 1:10_000
+            i += 1.0
+            tpr = Float64(nbinompdf(r, p, i))
+            pr += tpr
+            if pr >= cprob
+                return i
+            end
+        end
+        return i
+    end
+end
+
+# Inverse CCDF: find smallest k such that nbinomccdf(r, p, k) <= cprob.
+# Based on VBA critcompnegbinom by Ian Smith.
+function _critcompnbinom(r::Float64, p::Float64, cprob::Float64)
+    q = 1.0 - p
+
+    # Normal approximation
+    μ = r * q / p
+    σ = sqrt(μ / p)
+    i = max(0.0, floor(μ - norminvcdf(min(cprob, 1.0 - 1.0e-15)) * σ + 0.5))
+
+    # Compute CCDF at the guess
+    pr = Float64(nbinomccdf(r, p, i))
+
+    if pr > cprob
+        # Search right
+        for _ in 1:10_000
+            i += 1.0
+            tpr = Float64(nbinompdf(r, p, i))
+            pr -= tpr
+            if pr <= cprob
+                return i
+            end
+        end
+        return i
+    else
+        # Search left
+        while i > 0
+            tpr = Float64(nbinompdf(r, p, i))
+            if pr + tpr > cprob
+                return i
+            end
+            pr += tpr
+            i -= 1.0
+        end
+        return 0.0
+    end
+end
+
+# Wrapper with edge cases
+function _nbinom_invcdf(r::Float64, p::Float64, q::Float64)
+    if q < 0 || q > 1 || r <= 0 || p < 0 || p > 1 || isnan(q) || isnan(r) || isnan(p)
+        return NaN
+    elseif q == 0 || p == 1
+        return 0.0
+    elseif q == 1
+        return Inf
+    elseif p == 0
+        return Inf
+    end
+
+    i = _critnbinom(r, p, q)
+
+    # Post-correction
+    pr = Float64(nbinomcdf(r, p, i))
+    if pr >= q
+        while i > 0
+            pr2 = Float64(nbinomcdf(r, p, i - 1.0))
+            if pr2 < q
+                return i
+            end
+            i -= 1.0
+        end
+        return 0.0
+    else
+        return i + 1.0
+    end
+end
+
+function _nbinom_invccdf(r::Float64, p::Float64, q::Float64)
+    if q < 0 || q > 1 || r <= 0 || p < 0 || p > 1 || isnan(q) || isnan(r) || isnan(p)
+        return NaN
+    elseif q == 0 || p == 0
+        return Inf
+    elseif q == 1
+        return 0.0
+    elseif p == 1
+        return 0.0
+    end
+
+    i = _critcompnbinom(r, p, q)
+
+    # Post-correction
+    pr = Float64(nbinomccdf(r, p, i))
+    if pr <= q
+        while i > 0
+            pr2 = Float64(nbinomccdf(r, p, i - 1.0))
+            if pr2 > q
+                return i
+            end
+            i -= 1.0
+        end
+        return 0.0
+    else
+        return i + 1.0
+    end
+end
+
+# Public API
+
 function nbinominvcdf(r::Real, p::Real, q::Real)
     T = float(Base.promote_typeof(r, p, q))
-    return convert(T, Rmath.qnbinom(q, r, p, true, false))
+    return convert(T, _nbinom_invcdf(Float64(r), Float64(p), Float64(q)))
 end
+
 function nbinominvccdf(r::Real, p::Real, q::Real)
     T = float(Base.promote_typeof(r, p, q))
-    return convert(T, Rmath.qnbinom(q, r, p, false, false))
+    return convert(T, _nbinom_invccdf(Float64(r), Float64(p), Float64(q)))
 end
+
 function nbinominvlogcdf(r::Real, p::Real, lq::Real)
     T = float(Base.promote_typeof(r, p, lq))
-    return convert(T, Rmath.qnbinom(lq, r, p, true, true))
+    _lq = Float64(lq)
+    result = if _lq > -1
+        _nbinom_invccdf(Float64(r), Float64(p), -expm1(_lq))
+    else
+        _nbinom_invcdf(Float64(r), Float64(p), exp(_lq))
+    end
+    return convert(T, result)
 end
+
 function nbinominvlogccdf(r::Real, p::Real, lq::Real)
     T = float(Base.promote_typeof(r, p, lq))
-    return convert(T, Rmath.qnbinom(lq, r, p, false, true))
+    _lq = Float64(lq)
+    result = if _lq > -1
+        _nbinom_invcdf(Float64(r), Float64(p), -expm1(_lq))
+    else
+        _nbinom_invccdf(Float64(r), Float64(p), exp(_lq))
+    end
+    return convert(T, result)
 end
diff --git a/src/distrs/pois.jl b/src/distrs/pois.jl
index 1c9b565..d65c187 100644
--- a/src/distrs/pois.jl
+++ b/src/distrs/pois.jl
@@ -18,20 +18,163 @@ poislogcdf(λ::Real, x::Real) = gammalogccdf(max(0, floor(x + 1)), 1, λ)
 
 poislogccdf(λ::Real, x::Real) = gammalogcdf(max(0, floor(x + 1)), 1, λ)
 
-# Rmath implementations
+# Inverse CDF: find smallest k such that poiscdf(λ, k) >= cprob.
+# Based on VBA critpoiss by Ian Smith.
+# PMF ratio: PMF(k+1)/PMF(k) = λ/(k+1), PMF(k-1)/PMF(k) = k/λ
+function _critpoiss(λ::Float64, cprob::Float64)
+    # Normal approximation for initial guess
+    σ = sqrt(λ)
+    i = max(0.0, floor(λ + norminvcdf(min(cprob, 1.0 - 1.0e-15)) * σ + 0.5))
+
+    # Compute CDF at the guess
+    pr = Float64(poiscdf(λ, i))
+
+    if pr >= cprob
+        # Search left: find smallest k with CDF(k) >= cprob
+        while i > 0
+            tpr = Float64(poispdf(λ, i))
+            if pr - tpr < cprob
+                return i
+            end
+            pr -= tpr
+            i -= 1.0
+        end
+        return 0.0
+    else
+        # Search right (with guard against infinite loop)
+        for _ in 1:10_000
+            i += 1.0
+            tpr = Float64(poispdf(λ, i))
+            pr += tpr
+            if pr >= cprob
+                return i
+            end
+        end
+        return i
+    end
+end
+
+# Inverse CCDF: find smallest k such that poisccdf(λ, k) <= cprob.
+# Based on VBA critcomppoiss by Ian Smith.
+function _critcomppoiss(λ::Float64, cprob::Float64)
+    # Normal approximation
+    σ = sqrt(λ)
+    i = max(0.0, floor(λ - norminvcdf(min(cprob, 1.0 - 1.0e-15)) * σ + 0.5))
+
+    # Compute CCDF at the guess
+    pr = Float64(poisccdf(λ, i))
+
+    if pr > cprob
+        # Search right (with guard)
+        for _ in 1:10_000
+            i += 1.0
+            tpr = Float64(poispdf(λ, i))
+            pr -= tpr
+            if pr <= cprob
+                return i
+            end
+        end
+        return i
+    else
+        # Search left
+        while i > 0
+            tpr = Float64(poispdf(λ, i))
+            if pr + tpr > cprob
+                return i
+            end
+            pr += tpr
+            i -= 1.0
+        end
+        return 0.0
+    end
+end
+
+# Wrapper with edge cases
+function _pois_invcdf(λ::Float64, q::Float64)
+    if q < 0 || q > 1 || λ < 0 || isnan(q) || isnan(λ)
+        return NaN
+    elseif q == 0 || λ == 0
+        return 0.0
+    elseif q == 1
+        return Inf
+    end
+
+    i = _critpoiss(λ, q)
+
+    # Post-correction
+    pr = Float64(poiscdf(λ, i))
+    if pr >= q
+        while i > 0
+            pr2 = Float64(poiscdf(λ, i - 1.0))
+            if pr2 < q
+                return i
+            end
+            i -= 1.0
+        end
+        return 0.0
+    else
+        return i + 1.0
+    end
+end
+
+function _pois_invccdf(λ::Float64, q::Float64)
+    if q < 0 || q > 1 || λ < 0 || isnan(q) || isnan(λ)
+        return NaN
+    elseif q == 0
+        return Inf
+    elseif λ == 0
+        return 0.0
+    end
+
+    i = _critcomppoiss(λ, q)
+
+    # Post-correction
+    pr = Float64(poisccdf(λ, i))
+    if pr <= q
+        while i > 0
+            pr2 = Float64(poisccdf(λ, i - 1.0))
+            if pr2 > q
+                return i
+            end
+            i -= 1.0
+        end
+        return 0.0
+    else
+        return i + 1.0
+    end
+end
+
+# Public API
+
 function poisinvcdf(λ::Real, q::Real)
     T = float(Base.promote_typeof(λ, q))
-    return convert(T, Rmath.qpois(q, λ, true, false))
+    return convert(T, _pois_invcdf(Float64(λ), Float64(q)))
 end
+
 function poisinvccdf(λ::Real, q::Real)
     T = float(Base.promote_typeof(λ, q))
-    return convert(T, Rmath.qpois(q, λ, false, false))
+    return convert(T, _pois_invccdf(Float64(λ), Float64(q)))
 end
+
 function poisinvlogcdf(λ::Real, lq::Real)
     T = float(Base.promote_typeof(λ, lq))
-    return convert(T, Rmath.qpois(lq, λ, true, true))
+    _lq = Float64(lq)
+    # Use -expm1(lq) = 1-exp(lq) via invccdf for accuracy when lq ≈ 0
+    result = if _lq > -1
+        _pois_invccdf(Float64(λ), -expm1(_lq))
+    else
+        _pois_invcdf(Float64(λ), exp(_lq))
+    end
+    return convert(T, result)
 end
+
 function poisinvlogccdf(λ::Real, lq::Real)
     T = float(Base.promote_typeof(λ, lq))
-    return convert(T, Rmath.qpois(lq, λ, false, true))
+    _lq = Float64(lq)
+    result = if _lq > -1
+        _pois_invcdf(Float64(λ), -expm1(_lq))
+    else
+        _pois_invccdf(Float64(λ), exp(_lq))
+    end
+    return convert(T, result)
 end
diff --git a/test/rmath.jl b/test/rmath.jl
index 8d5bd60..9f8087e 100644
--- a/test/rmath.jl
+++ b/test/rmath.jl
@@ -219,8 +219,8 @@ end
         end
     end
 
-    rmathcomp_tests(
-        "binom", [
+    @testset "binom" begin
+        configs = [
             ((1, 0.5), 0.0:1.0),
             ((1, 0.7), 0.0:1.0),
             ((8, 0.6), 0.0:8.0),
@@ -234,7 +234,52 @@ end
             ((10, 1 // 2), (0 // 1):(10 // 1)),
             ((1, 0.5), [-Inf, Inf]),
         ]
-    )
+        @testset "params: $params" for (params, data) in configs
+            (n, p) = params
+            T = float(Base.promote_typeof(n, p, first(data)))
+            rtol = _default_rtol(T)
+            @testset "pdf with x=$x" for x in data
+                check_rmath(binompdf, (a, n, p) -> Rmath.dbinom(a, n, p, false), params, x, true, rtol)
+            end
+            @testset "logpdf with x=$x" for x in data
+                check_rmath(binomlogpdf, (a, n, p) -> Rmath.dbinom(a, n, p, true), params, x, false, rtol)
+            end
+            @testset "cdf with x=$x" for x in data
+                check_rmath(binomcdf, (a, n, p) -> Rmath.pbinom(a, n, p, true, false), params, x, true, rtol)
+            end
+            @testset "ccdf with x=$x" for x in data
+                check_rmath(binomccdf, (a, n, p) -> Rmath.pbinom(a, n, p, false, false), params, x, true, rtol)
+            end
+            @testset "logcdf with x=$x" for x in data
+                check_rmath(binomlogcdf, (a, n, p) -> Rmath.pbinom(a, n, p, true, true), params, x, false, rtol)
+            end
+            @testset "logccdf with x=$x" for x in data
+                check_rmath(binomlogccdf, (a, n, p) -> Rmath.pbinom(a, n, p, false, true), params, x, false, rtol)
+            end
+            if T === Float64
+                @testset "invcdf round-trip with x=$x" for x in data
+                    q = binomcdf(n, p, x)
+                    (q == 0 || q == 1) && continue
+                    @test binominvcdf(n, p, q) ≈ x rtol = rtol nans = true
+                end
+                @testset "invccdf round-trip with x=$x" for x in data
+                    q = binomccdf(n, p, x)
+                    (q == 0 || q == 1) && continue
+                    @test binominvccdf(n, p, q) ≈ x rtol = rtol nans = true
+                end
+                @testset "invlogcdf round-trip with x=$x" for x in data
+                    lq = binomlogcdf(n, p, x)
+                    (lq == 0 || lq == -Inf) && continue
+                    @test binominvlogcdf(n, p, lq) ≈ x atol = 1 rtol = rtol nans = true
+                end
+                @testset "invlogccdf round-trip with x=$x" for x in data
+                    lq = binomlogccdf(n, p, x)
+                    (lq == 0 || lq == -Inf) && continue
+                    @test binominvlogccdf(n, p, lq) ≈ x atol = 1 rtol = rtol nans = true
+                end
+            end
+        end
+    end
 
     rmathcomp_tests(
         "chisq", [
@@ -307,8 +352,8 @@ end
         ]
     )
 
-    rmathcomp_tests(
-        "nbinom", [
+    @testset "nbinom" begin
+        configs = [
             ((1, 0.5), 0.0:20.0),
             ((3, 0.5), 0.0:20.0),
             ((3, 0.2), 0.0:20.0),
@@ -321,7 +366,60 @@ end
             ((1, 0.5), [-Inf, Inf, NaN]),
             ((1, 0.5f0), [-Inf32, Inf32, NaN32, -2, 1.1]),
         ]
-    )
+        @testset "params: $params" for (params, data) in configs
+            (r, p) = params
+            T = float(Base.promote_typeof(r, p, first(data)))
+            rtol = _default_rtol(T)
+            @testset "pdf with x=$x" for x in data
+                check_rmath(nbinompdf, (a, r, p) -> Rmath.dnbinom(a, r, p, false), params, x, true, rtol)
+            end
+            @testset "logpdf with x=$x" for x in data
+                check_rmath(nbinomlogpdf, (a, r, p) -> Rmath.dnbinom(a, r, p, true), params, x, false, rtol)
+            end
+            @testset "cdf with x=$x" for x in data
+                check_rmath(nbinomcdf, (a, r, p) -> Rmath.pnbinom(a, r, p, true, false), params, x, true, rtol)
+            end
+            @testset "ccdf with x=$x" for x in data
+                check_rmath(nbinomccdf, (a, r, p) -> Rmath.pnbinom(a, r, p, false, false), params, x, true, rtol)
+            end
+            @testset "logcdf with x=$x" for x in data
+                check_rmath(nbinomlogcdf, (a, r, p) -> Rmath.pnbinom(a, r, p, true, true), params, x, false, rtol)
+            end
+            @testset "logccdf with x=$x" for x in data
+                check_rmath(nbinomlogccdf, (a, r, p) -> Rmath.pnbinom(a, r, p, false, true), params, x, false, rtol)
+            end
+            if T === Float64
+                @testset "invcdf round-trip with x=$x" for x in data
+                    (isnan(x) || isinf(x)) && continue
+                    xf = floor(Float64(x))
+                    q = nbinomcdf(r, p, x)
+                    (q == 0 || q == 1 || xf < 0) && continue
+                    @test Float64(nbinominvcdf(r, p, q)) == xf
+                end
+                @testset "invccdf round-trip with x=$x" for x in data
+                    (isnan(x) || isinf(x)) && continue
+                    xf = floor(Float64(x))
+                    q = nbinomccdf(r, p, x)
+                    (q == 0 || q == 1 || xf < 0) && continue
+                    @test Float64(nbinominvccdf(r, p, q)) == xf
+                end
+                @testset "invlogcdf round-trip with x=$x" for x in data
+                    (isnan(x) || isinf(x)) && continue
+                    xf = floor(Float64(x))
+                    lq = nbinomlogcdf(r, p, x)
+                    (!isfinite(lq) || lq == 0 || xf < 0) && continue
+                    @test abs(Float64(nbinominvlogcdf(r, p, lq)) - xf) <= 1
+                end
+                @testset "invlogccdf round-trip with x=$x" for x in data
+                    (isnan(x) || isinf(x)) && continue
+                    xf = floor(Float64(x))
+                    lq = nbinomlogccdf(r, p, x)
+                    (!isfinite(lq) || lq == 0 || xf < 0) && continue
+                    @test abs(Float64(nbinominvlogccdf(r, p, lq)) - xf) <= 1
+                end
+            end
+        end
+    end
 
     rmathcomp_tests(
         "nchisq", [
@@ -386,8 +484,8 @@ end
         ]
     )
 
-    rmathcomp_tests(
-        "pois", [
+    @testset "pois" begin
+        configs = [
             ((0.5,), 0.0:20.0),
             ((1.0,), 0.0:20.0),
             ((2.0,), 0.0:20.0),
@@ -399,7 +497,52 @@ end
             ((0.5f0,), Float16(0):Float16(20)),
             ((1 // 2,), (0 // 1):(20 // 1)),
         ]
-    )
+        @testset "params: $params" for (params, data) in configs
+            (λ,) = params
+            T = float(Base.promote_typeof(λ, first(data)))
+            rtol = _default_rtol(T)
+            @testset "pdf with x=$x" for x in data
+                check_rmath(poispdf, (a, λ) -> Rmath.dpois(a, λ, false), params, x, true, rtol)
+            end
+            @testset "logpdf with x=$x" for x in data
+                check_rmath(poislogpdf, (a, λ) -> Rmath.dpois(a, λ, true), params, x, false, rtol)
+            end
+            @testset "cdf with x=$x" for x in data
+                check_rmath(poiscdf, (a, λ) -> Rmath.ppois(a, λ, true, false), params, x, true, rtol)
+            end
+            @testset "ccdf with x=$x" for x in data
+                check_rmath(poisccdf, (a, λ) -> Rmath.ppois(a, λ, false, false), params, x, true, rtol)
+            end
+            @testset "logcdf with x=$x" for x in data
+                check_rmath(poislogcdf, (a, λ) -> Rmath.ppois(a, λ, true, true), params, x, false, rtol)
+            end
+            @testset "logccdf with x=$x" for x in data
+                check_rmath(poislogccdf, (a, λ) -> Rmath.ppois(a, λ, false, true), params, x, false, rtol)
+            end
+            if T === Float64
+                @testset "invcdf round-trip with x=$x" for x in data
+                    q = poiscdf(λ, x)
+                    (q == 0 || q == 1) && continue
+                    @test poisinvcdf(λ, q) ≈ x rtol = rtol nans = true
+                end
+                @testset "invccdf round-trip with x=$x" for x in data
+                    q = poisccdf(λ, x)
+                    (q == 0 || q == 1) && continue
+                    @test poisinvccdf(λ, q) ≈ x rtol = rtol nans = true
+                end
+                @testset "invlogcdf round-trip with x=$x" for x in data
+                    lq = poislogcdf(λ, x)
+                    (lq == 0 || lq == -Inf) && continue
+                    @test poisinvlogcdf(λ, lq) ≈ x atol = 1 rtol = rtol nans = true
+                end
+                @testset "invlogccdf round-trip with x=$x" for x in data
+                    lq = poislogccdf(λ, x)
+                    (lq == 0 || lq == -Inf) && continue
+                    @test poisinvlogccdf(λ, lq) ≈ x atol = 1 rtol = rtol nans = true
+                end
+            end
+        end
+    end
 
     rmathcomp_tests(
         "tdist", [
@@ -483,13 +626,36 @@ end
             ((1, 1), [-10, -6, 2, 24]),
         ]
     )
-    rmathcomp_tests(
-        "binom", [
+    @testset "binom (edge cases)" begin
+        configs = [
             ((5, 0.5), [-8, -2.3, 1.2, 5.4, 11.9]),
             ((5, 1 // 2), [-8, -23 // 10, 6 // 5, 27 // 5, 119 // 10]),
             ((5, 1 // 2), [-8, -2, 6, 12]),
         ]
-    )
+        @testset "params: $params" for (params, data) in configs
+            (n, p) = params
+            T = float(Base.promote_typeof(n, p, first(data)))
+            rtol = _default_rtol(T)
+            @testset "pdf with x=$x" for x in data
+                check_rmath(binompdf, (a, n, p) -> Rmath.dbinom(a, n, p, false), params, x, true, rtol)
+            end
+            @testset "logpdf with x=$x" for x in data
+                check_rmath(binomlogpdf, (a, n, p) -> Rmath.dbinom(a, n, p, true), params, x, false, rtol)
+            end
+            @testset "cdf with x=$x" for x in data
+                check_rmath(binomcdf, (a, n, p) -> Rmath.pbinom(a, n, p, true, false), params, x, true, rtol)
+            end
+            @testset "ccdf with x=$x" for x in data
+                check_rmath(binomccdf, (a, n, p) -> Rmath.pbinom(a, n, p, false, false), params, x, true, rtol)
+            end
+            @testset "logcdf with x=$x" for x in data
+                check_rmath(binomlogcdf, (a, n, p) -> Rmath.pbinom(a, n, p, true, true), params, x, false, rtol)
+            end
+            @testset "logccdf with x=$x" for x in data
+                check_rmath(binomlogccdf, (a, n, p) -> Rmath.pbinom(a, n, p, false, true), params, x, false, rtol)
+            end
+        end
+    end
     rmathcomp_tests(
         "fdist", [
             ((1.0, 1.0), [-10.0, -6.3]),
@@ -504,11 +670,34 @@ end
             ((1, 1), [-10, -6]),
         ]
     )
-    rmathcomp_tests(
-        "pois", [
+    @testset "pois (edge cases)" begin
+        configs = [
             ((0.5,), [-10, -2.5, 1.3, 8.7]),
             ((1 // 2,), [-10, -5 // 2, 13 // 10, 87 // 10]),
             ((1,), [-10, -3]),
         ]
-    )
+        @testset "params: $params" for (params, data) in configs
+            (λ,) = params
+            T = float(Base.promote_typeof(λ, first(data)))
+            rtol = _default_rtol(T)
+            @testset "pdf with x=$x" for x in data
+                check_rmath(poispdf, (a, λ) -> Rmath.dpois(a, λ, false), params, x, true, rtol)
+            end
+            @testset "logpdf with x=$x" for x in data
+                check_rmath(poislogpdf, (a, λ) -> Rmath.dpois(a, λ, true), params, x, false, rtol)
+            end
+            @testset "cdf with x=$x" for x in data
+                check_rmath(poiscdf, (a, λ) -> Rmath.ppois(a, λ, true, false), params, x, true, rtol)
+            end
+            @testset "ccdf with x=$x" for x in data
+                check_rmath(poisccdf, (a, λ) -> Rmath.ppois(a, λ, false, false), params, x, true, rtol)
+            end
+            @testset "logcdf with x=$x" for x in data
+                check_rmath(poislogcdf, (a, λ) -> Rmath.ppois(a, λ, true, true), params, x, false, rtol)
+            end
+            @testset "logccdf with x=$x" for x in data
+                check_rmath(poislogccdf, (a, λ) -> Rmath.ppois(a, λ, false, true), params, x, false, rtol)
+            end
+        end
+    end
 end