diff --git a/src/StatsFuns.jl b/src/StatsFuns.jl
index 41a8a5d..3e51787 100644
--- a/src/StatsFuns.jl
+++ b/src/StatsFuns.jl
@@ -217,6 +217,16 @@ export
     tdistinvlogcdf,     # inverse-logcdf of student's t distribution
     tdistinvlogccdf,    # inverse-logccdf of student's t distribution
 
+    # distrs/srdist
+    srdistcdf,           # cdf of studentized range distribution
+    srdistccdf,          # ccdf of studentized range distribution
+    srdistlogcdf,        # logcdf of studentized range distribution
+    srdistlogccdf,       # logccdf of studentized range distribution
+    srdistinvcdf,        # inverse-cdf of studentized range distribution
+    srdistinvccdf,       # inverse-ccdf of studentized range distribution
+    srdistinvlogcdf,     # inverse-logcdf of studentized range distribution
+    srdistinvlogccdf,    # inverse-logccdf of studentized range distribution
+
     # misc
     logmvgamma,         # logarithm of multivariate gamma function
     lstirling_asym
@@ -244,5 +254,6 @@ include(joinpath("distrs", "norm.jl"))
 include(joinpath("distrs", "ntdist.jl"))
 include(joinpath("distrs", "pois.jl"))
 include(joinpath("distrs", "tdist.jl"))
+include(joinpath("distrs", "srdist.jl"))
 
 end # module
diff --git a/src/distrs/fdist.jl b/src/distrs/fdist.jl
index b0b9fdc..b46c796 100644
--- a/src/distrs/fdist.jl
+++ b/src/distrs/fdist.jl
@@ -13,7 +13,7 @@ import .RFunctions:
     fdistinvlogccdf
 
 # pdf for numbers with generic types
-fdistpdf(d1::Real, d2::Real, x::Number) = sqrt((d1 * x)^d1 * d2^d2 / (d1 * x + d2)^(d1 + d2)) / (x * beta(d1 / 2, d2 / 2))
+fdistpdf(ν1::Real, ν2::Real, x::Number) = sqrt((ν1 * x)^ν1 * ν2^ν2 / (ν1 * x + ν2)^(ν1 + ν2)) / (x * beta(ν1 / 2, ν2 / 2))
 
 # logpdf for numbers with generic types
-fdistlogpdf(d1::Real, d2::Real, x::Number) = (d1 * log(d1 * x) + d2 * log(d2) - (d1 + d2) * log(d1 * x + d2)) / 2 - log(x) - lbeta(d1 / 2, d2 / 2)
+fdistlogpdf(ν1::Real, ν2::Real, x::Number) = (ν1 * log(ν1 * x) + ν2 * log(ν2) - (ν1 + ν2) * log(ν1 * x + ν2)) / 2 - log(x) - lbeta(ν1 / 2, ν2 / 2)
diff --git a/src/distrs/srdist.jl b/src/distrs/srdist.jl
new file mode 100644
index 0000000..8d1e180
--- /dev/null
+++ b/src/distrs/srdist.jl
@@ -0,0 +1,11 @@
+# functions related to studentized range distribution
+
+import .RFunctions:
+    srdistcdf,
+    srdistccdf,
+    srdistlogcdf,
+    srdistlogccdf,
+    srdistinvcdf,
+    srdistinvccdf,
+    srdistinvlogcdf,
+    srdistinvlogccdf
diff --git a/src/rmath.jl b/src/rmath.jl
index bb85bdb..e2b4c77 100644
--- a/src/rmath.jl
+++ b/src/rmath.jl
@@ -34,29 +34,45 @@ function _import_rmath(rname::Symbol, jname::Symbol, pargs)
     invlogcdf = Symbol(jname, "invlogcdf")
     invlogccdf = Symbol(jname, "invlogccdf")
 
+    is_tukey = rname == :tukey
+
     rand = Symbol(jname, "rand")
     has_rand = true
-    if rname == :nbeta || rname == :nf || rname == :nt
+    if rname == :nbeta || rname == :nf || rname == :nt || is_tukey
         has_rand = false
     end
 
     # arguments & argument types
     _pts = fill(:Cdouble, length(pargs))
 
-    dtypes = Expr(:tuple, :Cdouble, _pts..., :Cint)
-    ptypes = Expr(:tuple, :Cdouble, _pts..., :Cint, :Cint)
-    qtypes = Expr(:tuple, :Cdouble, _pts..., :Cint, :Cint)
-    rtypes = Expr(:tuple, _pts...)
+    if is_tukey
+        dtypes = Expr(:tuple, :Cdouble, :Cdouble, _pts..., :Cint)
+        ptypes = Expr(:tuple, :Cdouble, :Cdouble, _pts..., :Cint, :Cint)
+        qtypes = Expr(:tuple, :Cdouble, :Cdouble, _pts..., :Cint, :Cint)
+        rtypes = Expr(:tuple, :Cdouble, _pts...)
+    else
+        dtypes = Expr(:tuple, :Cdouble, _pts..., :Cint)
+        ptypes = Expr(:tuple, :Cdouble, _pts..., :Cint, :Cint)
+        qtypes = Expr(:tuple, :Cdouble, _pts..., :Cint, :Cint)
+        rtypes = Expr(:tuple, _pts...)
+    end
 
     pdecls = [Expr(:(::), ps, :(Union{Float64,Int})) for ps in pargs] # [:(p1::Union{Float64, Int}), :(p2::Union{...}), ...]
 
-    # function implementation
+    if is_tukey
+        # ptukey and qtukey have an extra literal 1 argument
+        pargs = (1, pargs...)
+    end
+
+    # Function implementation
     quote
-        $pdf($(pdecls...), x::Union{Float64,Int}) =
-            ccall(($dfun, libRmath), Float64, $dtypes, x, $(pargs...), 0)
+        if $(!is_tukey)
+            $pdf($(pdecls...), x::Union{Float64,Int}) =
+                ccall(($dfun, libRmath), Float64, $dtypes, x, $(pargs...), 0)
 
-        $logpdf($(pdecls...), x::Union{Float64,Int}) =
-            ccall(($dfun, libRmath), Float64, $dtypes, x, $(pargs...), 1)
+            $logpdf($(pdecls...), x::Union{Float64,Int}) =
+                ccall(($dfun, libRmath), Float64, $dtypes, x, $(pargs...), 1)
+        end
 
         $cdf($(pdecls...), x::Union{Float64,Int}) =
             ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 1, 0)
@@ -99,13 +115,14 @@ end
 @import_rmath beta beta α β
 @import_rmath binom binom n p
 @import_rmath chisq chisq k
-@import_rmath f fdist d1 d2
+@import_rmath f fdist ν1 ν2
 @import_rmath gamma gamma α β
 @import_rmath hyper hyper ms mf n
 @import_rmath norm norm μ σ
 @import_rmath nbinom nbinom r p
 @import_rmath pois pois λ
 @import_rmath t tdist k
+@import_rmath tukey srdist k ν
 
 @import_rmath nbeta nbeta α β λ
 @import_rmath nchisq nchisq k λ
diff --git a/test/generic.jl b/test/generic.jl
index cd5f717..1b4bf0e 100644
--- a/test/generic.jl
+++ b/test/generic.jl
@@ -43,85 +43,95 @@ end
 
 ### Test cases
 
-genericcomp_tests("beta", [
-    ((1.0, 1.0), 0.01:0.01:0.99),
-    ((2.0, 3.0), 0.01:0.01:0.99),
-    ((10.0, 2.0), 0.01:0.01:0.99),
-])
-
-genericcomp_tests("binom", [
-    ((1, 0.5), 0.0:1.0),
-    ((1, 0.7), 0.0:1.0),
-    ((8, 0.6), 0.0:8.0),
-    ((20, 0.1), 0.0:20.0),
-    ((20, 0.9), 0.0:20.0),
-    ((20, 0.9), 0:20),
-])
-
-genericcomp_tests("chisq", [
-    ((1,), 0.0:0.1:8.0),
-    ((4,), 0.0:0.1:8.0),
-    ((9,), 0.0:0.1:8.0),
-])
-
-genericcomp_tests("fdist", [
-    ((1, 1), (0.0:0.1:5.0)),
-    ((2, 1), (0.0:0.1:5.0)),
-    ((5, 2), (0.0:0.1:5.0)),
-    ((10, 1), (0.0:0.1:5.0)),
-    ((10, 3), (0.0:0.1:5.0)),
-])
-
-genericcomp_tests("gamma", [
-    ((1.0, 1.0), (0.05:0.05:12.0)),
-    ((0.5, 1.0), (0.05:0.05:12.0)),
-    ((3.0, 1.0), (0.05:0.05:12.0)),
-    ((9.0, 1.0), (0.05:0.05:12.0)),
-    ((2.0, 3.0), (0.05:0.05:12.0)),
-])
-
-# genericcomp_tests("nbeta", [
-#     ((1.0, 1.0, 0.0), 0.01:0.01:0.99),
-#     ((2.0, 3.0, 0.0), 0.01:0.01:0.99),
-#     ((1.0, 1.0, 2.0), 0.01:0.01:0.99),
-#     ((3.0, 4.0, 2.0), 0.01:0.01:0.99),
-# ])
-
-# genericcomp_tests("nchisq", [
-#     ((2, 1), 0.0:0.2:8.0),
-#     ((2, 3), 0.0:0.2:8.0),
-#     ((4, 1), 0.0:0.2:8.0),
-#     ((4, 3), 0.0:0.2:8.0),
-# ])
-
-# genericcomp_tests("nfdist", [
-#     ((1.0, 1.0, 0.0), 0.1:0.1:10.0),
-#     ((1.0, 1.0, 2.0), 0.1:0.1:10.0),
-#     ((2.0, 3.0, 1.0), 0.1:0.1:10.0),
-# ])
-
-genericcomp_tests("norm", [
-    ((0.0, 1.0), -6.0:0.01:6.0),
-    ((2.0, 1.0), -3.0:0.01:7.0),
-    ((0.0, 0.5), -3.0:0.01:3.0),
-    ((0, 1), -3.0:3.0),
-    ((0, 1), -3.0:0.01:3.0)
-])
-
-# genericcomp_tests("ntdist", [
-#     ((0, 1), -4.0:0.1:10.0),
-#     ((0, 4), -4.0:0.1:10.0),
-#     ((2, 1), -4.0:0.1:10.0),
-#     ((2, 4), -4.0:0.1:10.0),
-# ])
-
-genericcomp_tests("pois", [
-    ((1.0,), 0:30),
-    ((10.0,), 0:42)
-])
-
-genericcomp_tests("tdist", [
-    ((1,), -5.0:0.1:5.0),
-    ((2,), -5.0:0.1:5.0),
-    ((5,), -5.0:0.1:5.0),
-])
+@testset "Generic" begin
+    genericcomp_tests("beta", [
+        ((1.0, 1.0), 0.01:0.01:0.99),
+        ((2.0, 3.0), 0.01:0.01:0.99),
+        ((10.0, 2.0), 0.01:0.01:0.99),
+    ])
+
+    genericcomp_tests("binom", [
+        ((1, 0.5), 0.0:1.0),
+        ((1, 0.7), 0.0:1.0),
+        ((8, 0.6), 0.0:8.0),
+        ((20, 0.1), 0.0:20.0),
+        ((20, 0.9), 0.0:20.0),
+        ((20, 0.9), 0:20),
+    ])
+
+    genericcomp_tests("chisq", [
+        ((1,), 0.0:0.1:8.0),
+        ((4,), 0.0:0.1:8.0),
+        ((9,), 0.0:0.1:8.0),
+    ])
+
+    genericcomp_tests("fdist", [
+        ((1, 1), (0.0:0.1:5.0)),
+        ((2, 1), (0.0:0.1:5.0)),
+        ((5, 2), (0.0:0.1:5.0)),
+        ((10, 1), (0.0:0.1:5.0)),
+        ((10, 3), (0.0:0.1:5.0)),
+    ])
+
+    genericcomp_tests("gamma", [
+        ((1.0, 1.0), (0.05:0.05:12.0)),
+        ((0.5, 1.0), (0.05:0.05:12.0)),
+        ((3.0, 1.0), (0.05:0.05:12.0)),
+        ((9.0, 1.0), (0.05:0.05:12.0)),
+        ((2.0, 3.0), (0.05:0.05:12.0)),
+    ])
+
+    # genericcomp_tests("nbeta", [
+    #     ((1.0, 1.0, 0.0), 0.01:0.01:0.99),
+    #     ((2.0, 3.0, 0.0), 0.01:0.01:0.99),
+    #     ((1.0, 1.0, 2.0), 0.01:0.01:0.99),
+    #     ((3.0, 4.0, 2.0), 0.01:0.01:0.99),
+    # ])
+
+    # genericcomp_tests("nchisq", [
+    #     ((2, 1), 0.0:0.2:8.0),
+    #     ((2, 3), 0.0:0.2:8.0),
+    #     ((4, 1), 0.0:0.2:8.0),
+    #     ((4, 3), 0.0:0.2:8.0),
+    # ])
+
+    # genericcomp_tests("nfdist", [
+    #     ((1.0, 1.0, 0.0), 0.1:0.1:10.0),
+    #     ((1.0, 1.0, 2.0), 0.1:0.1:10.0),
+    #     ((2.0, 3.0, 1.0), 0.1:0.1:10.0),
+    # ])
+
+    genericcomp_tests("norm", [
+        ((0.0, 1.0), -6.0:0.01:6.0),
+        ((2.0, 1.0), -3.0:0.01:7.0),
+        ((0.0, 0.5), -3.0:0.01:3.0),
+        ((0, 1), -3.0:3.0),
+        ((0, 1), -3.0:0.01:3.0)
+    ])
+
+    # genericcomp_tests("ntdist", [
+    #     ((0, 1), -4.0:0.1:10.0),
+    #     ((0, 4), -4.0:0.1:10.0),
+    #     ((2, 1), -4.0:0.1:10.0),
+    #     ((2, 4), -4.0:0.1:10.0),
+    # ])
+
+    genericcomp_tests("pois", [
+        ((1.0,), 0:30),
+        ((10.0,), 0:42)
+    ])
+
+    genericcomp_tests("tdist", [
+        ((1,), -5.0:0.1:5.0),
+        ((2,), -5.0:0.1:5.0),
+        ((5,), -5.0:0.1:5.0),
+    ])
+
+    #genericcomp_tests("srdist", [
+    #    ((1,2), (0.0:0.1:5.0)),
+    #    ((2,2), (0.0:0.1:5.0)),
+    #    ((5,3), (0.0:0.1:5.0)),
+    #    ((10,2), (0.0:0.1:5.0)),
+    #    ((10,5), (0.0:0.1:5.0))
+    #])
+end
diff --git a/test/rmath.jl b/test/rmath.jl
index 4926213..138cce8 100644
--- a/test/rmath.jl
+++ b/test/rmath.jl
@@ -15,8 +15,16 @@ get_statsfun(fname) = eval(Symbol(fname))
 get_rmathfun(fname) = eval(Meta.parse(string("RFunctions.", fname)))
 
 function rmathcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(eltype(X)))))
-    pdf     = string(basename, "pdf")
-    logpdf  = string(basename, "logpdf")
+    # tackle pdf specially
+    has_pdf = true
+    if basename == "srdist"
+        has_pdf = false
+    end
+
+    if has_pdf
+        pdf     = string(basename, "pdf")
+        logpdf  = string(basename, "logpdf")
+    end
     cdf     = string(basename, "cdf")
     ccdf    = string(basename, "ccdf")
     logcdf  = string(basename, "logcdf")
@@ -27,8 +35,10 @@ function rmathcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(elt
     invlogccdf = string(basename, "invlogccdf")
     rand       = string(basename, "rand")
 
-    stats_pdf     = get_statsfun(pdf)
-    stats_logpdf  = get_statsfun(logpdf)
+    if has_pdf
+        stats_pdf     = get_statsfun(pdf)
+        stats_logpdf  = get_statsfun(logpdf)
+    end
     stats_cdf     = get_statsfun(cdf)
     stats_ccdf    = get_statsfun(ccdf)
     stats_logcdf  = get_statsfun(logcdf)
@@ -38,8 +48,10 @@ function rmathcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(elt
     stats_invlogcdf  = get_statsfun(invlogcdf)
     stats_invlogccdf = get_statsfun(invlogccdf)
 
-    rmath_pdf     = get_rmathfun(pdf)
-    rmath_logpdf  = get_rmathfun(logpdf)
+    if has_pdf
+        rmath_pdf     = get_rmathfun(pdf)
+        rmath_logpdf  = get_rmathfun(logpdf)
+    end
     rmath_cdf     = get_rmathfun(cdf)
     rmath_ccdf    = get_rmathfun(ccdf)
     rmath_logcdf  = get_rmathfun(logcdf)
@@ -51,16 +63,18 @@ function rmathcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(elt
 
     # tackle rand specially
     has_rand = true
-    if basename == "nbeta" || basename == "nfdist" || basename == "ntdist"
+    if basename == "nbeta" || basename == "nfdist" || basename == "ntdist" || basename == "srdist"
         has_rand = false
     end
     rmath_rand = has_rand ? get_rmathfun(rand) : nothing
 
     for x in X
-        check_rmath(pdf, stats_pdf, rmath_pdf,
-            params, "x", x, true, rtol)
-        check_rmath(logpdf, stats_logpdf, rmath_logpdf,
-            params, "x", x, false, rtol)
+        if has_pdf
+            check_rmath(pdf, stats_pdf, rmath_pdf,
+                params, "x", x, true, rtol)
+            check_rmath(logpdf, stats_logpdf, rmath_logpdf,
+                params, "x", x, false, rtol)
+        end
         check_rmath(cdf, stats_cdf, rmath_cdf,
             params, "x", x, true, rtol)
         check_rmath(ccdf, stats_ccdf, rmath_ccdf,
@@ -100,112 +114,130 @@ end
 
 ### Test cases
 
-rmathcomp_tests("beta", [
-    ((1.0, 1.0), 0.01:0.01:0.99),
-    ((2.0, 3.0), 0.01:0.01:0.99),
-    ((10.0, 2.0), 0.01:0.01:0.99),
-    ((10, 2), 0.01:0.01:0.99),
-])
-
-rmathcomp_tests("binom", [
-    ((1, 0.5), 0.0:1.0),
-    ((1, 0.7), 0.0:1.0),
-    ((8, 0.6), 0.0:8.0),
-    ((20, 0.1), 0.0:20.0),
-    ((20, 0.9), 0.0:20.0),
-    ((20, 0.9), 0:20),
-])
-
-rmathcomp_tests("chisq", [
-    ((1,), 0.0:0.1:8.0),
-    ((4,), 0.0:0.1:8.0),
-    ((9,), 0.0:0.1:8.0),
-    ((9,), 0:8),
-])
-
-rmathcomp_tests("fdist", [
-    ((1, 1), (0.0:0.1:5.0)),
-    ((2, 1), (0.0:0.1:5.0)),
-    ((5, 2), (0.0:0.1:5.0)),
-    ((10, 1), (0.0:0.1:5.0)),
-    ((10, 3), (0.0:0.1:5.0)),
-    ((10, 3), (0:5)),
-])
-
-rmathcomp_tests("gamma", [
-    ((1.0, 1.0), (0.05:0.05:12.0)),
-    ((0.5, 1.0), (0.05:0.05:12.0)),
-    ((3.0, 1.0), (0.05:0.05:12.0)),
-    ((9.0, 1.0), (0.05:0.05:12.0)),
-    ((2.0, 3.0), (0.05:0.05:12.0)),
-    ((2, 3), (1:12)),
-])
-
-rmathcomp_tests("hyper", [
-    ((2, 3, 4), 0.0:4.0),
-    ((2, 3, 4), 0:4)
-])
-
-rmathcomp_tests("nbeta", [
-    ((1.0, 1.0, 0.0), 0.01:0.01:0.99),
-    ((2.0, 3.0, 0.0), 0.01:0.01:0.99),
-    ((1.0, 1.0, 2.0), 0.01:0.01:0.99),
-    ((3.0, 4.0, 2.0), 0.01:0.01:0.99),
-    ((3, 4, 2), 0.01:0.01:0.99),
-])
-
-rmathcomp_tests("nbinom", [
-    ((1, 0.5), 0.0:20.0),
-    ((3, 0.5), 0.0:20.0),
-    ((3, 0.2), 0.0:20.0),
-    ((3, 0.8), 0.0:20.0),
-    ((3, 0.8), 0:20),
-])
-
-rmathcomp_tests("nchisq", [
-    ((2, 1), 0.0:0.2:8.0),
-    ((2, 3), 0.0:0.2:8.0),
-    ((4, 1), 0.0:0.2:8.0),
-    ((4, 3), 0.0:0.2:8.0),
-    ((4, 3), 0:8),
-])
-
-rmathcomp_tests("nfdist", [
-    ((1.0, 1.0, 0.0), 0.1:0.1:10.0),
-    ((1.0, 1.0, 2.0), 0.1:0.1:10.0),
-    ((2.0, 3.0, 1.0), 0.1:0.1:10.0),
-    ((2, 3, 1), 1:10),
-])
-
-rmathcomp_tests("norm", [
-    ((0.0, 1.0), -6.0:0.01:6.0),
-    ((2.0, 1.0), -3.0:0.01:7.0),
-    ((0.0, 0.5), -3.0:0.01:3.0),
-    ((0, 1), -3.0:3.0),
-    ((0, 1), -3.0:0.01:3.0),
-    ((0, 1), -3:3),
-    ((0.0, 0.0), -6.0:0.1:6.0),
-])
-
-rmathcomp_tests("ntdist", [
-    ((0, 1), -4.0:0.1:10.0),
-    ((0, 4), -4.0:0.1:10.0),
-    ((2, 1), -4.0:0.1:10.0),
-    ((2, 4), -4.0:0.1:10.0),
-    ((2, 4), -4:10),
-])
-
-rmathcomp_tests("pois", [
-    ((0.5,), 0.0:20.0),
-    ((1.0,), 0.0:20.0),
-    ((2.0,), 0.0:20.0),
-    ((10.0,), 0.0:20.0),
-    ((10,), 0:20),
-])
-
-rmathcomp_tests("tdist", [
-    ((1,), -5.0:0.1:5.0),
-    ((2,), -5.0:0.1:5.0),
-    ((5,), -5.0:0.1:5.0),
-    ((5,), -5:5),
-])
+@testset "RMath" begin
+    rmathcomp_tests("beta", [
+        ((1.0, 1.0), 0.01:0.01:0.99),
+        ((2.0, 3.0), 0.01:0.01:0.99),
+        ((10.0, 2.0), 0.01:0.01:0.99),
+        ((10, 2), 0.01:0.01:0.99),
+    ])
+
+    rmathcomp_tests("binom", [
+        ((1, 0.5), 0.0:1.0),
+        ((1, 0.7), 0.0:1.0),
+        ((8, 0.6), 0.0:8.0),
+        ((20, 0.1), 0.0:20.0),
+        ((20, 0.9), 0.0:20.0),
+        ((20, 0.9), 0:20),
+    ])
+
+    rmathcomp_tests("chisq", [
+        ((1,), 0.0:0.1:8.0),
+        ((4,), 0.0:0.1:8.0),
+        ((9,), 0.0:0.1:8.0),
+        ((9,), 0:8),
+    ])
+
+    rmathcomp_tests("fdist", [
+        ((1, 1), (0.0:0.1:5.0)),
+        ((2, 1), (0.0:0.1:5.0)),
+        ((5, 2), (0.0:0.1:5.0)),
+        ((10, 1), (0.0:0.1:5.0)),
+        ((10, 3), (0.0:0.1:5.0)),
+        ((10, 3), (0:5)),
+    ])
+
+    rmathcomp_tests("gamma", [
+        ((1.0, 1.0), (0.05:0.05:12.0)),
+        ((0.5, 1.0), (0.05:0.05:12.0)),
+        ((3.0, 1.0), (0.05:0.05:12.0)),
+        ((9.0, 1.0), (0.05:0.05:12.0)),
+        ((2.0, 3.0), (0.05:0.05:12.0)),
+        ((2, 3), (1:12)),
+    ])
+
+    rmathcomp_tests("hyper", [
+        ((2, 3, 4), 0.0:4.0),
+        ((2, 3, 4), 0:4)
+    ])
+
+    rmathcomp_tests("nbeta", [
+        ((1.0, 1.0, 0.0), 0.01:0.01:0.99),
+        ((2.0, 3.0, 0.0), 0.01:0.01:0.99),
+        ((1.0, 1.0, 2.0), 0.01:0.01:0.99),
+        ((3.0, 4.0, 2.0), 0.01:0.01:0.99),
+        ((3, 4, 2), 0.01:0.01:0.99),
+    ])
+
+    rmathcomp_tests("nbinom", [
+        ((1, 0.5), 0.0:20.0),
+        ((3, 0.5), 0.0:20.0),
+        ((3, 0.2), 0.0:20.0),
+        ((3, 0.8), 0.0:20.0),
+        ((3, 0.8), 0:20),
+    ])
+
+    rmathcomp_tests("nchisq", [
+        ((2, 1), 0.0:0.2:8.0),
+        ((2, 3), 0.0:0.2:8.0),
+        ((4, 1), 0.0:0.2:8.0),
+        ((4, 3), 0.0:0.2:8.0),
+        ((4, 3), 0:8),
+    ])
+
+    rmathcomp_tests("nfdist", [
+        ((1.0, 1.0, 0.0), 0.1:0.1:10.0),
+        ((1.0, 1.0, 2.0), 0.1:0.1:10.0),
+        ((2.0, 3.0, 1.0), 0.1:0.1:10.0),
+        ((2, 3, 1), 1:10),
+    ])
+
+    rmathcomp_tests("norm", [
+        ((0.0, 1.0), -6.0:0.01:6.0),
+        ((2.0, 1.0), -3.0:0.01:7.0),
+        ((0.0, 0.5), -3.0:0.01:3.0),
+        ((0, 1), -3.0:3.0),
+        ((0, 1), -3.0:0.01:3.0),
+        ((0, 1), -3:3),
+        ((0.0, 0.0), -6.0:0.1:6.0),
+    ])
+
+    rmathcomp_tests("ntdist", [
+        ((0, 1), -4.0:0.1:10.0),
+        ((0, 4), -4.0:0.1:10.0),
+        ((2, 1), -4.0:0.1:10.0),
+        ((2, 4), -4.0:0.1:10.0),
+        ((2, 4), -4:10),
+    ])
+
+    rmathcomp_tests("pois", [
+        ((0.5,), 0.0:20.0),
+        ((1.0,), 0.0:20.0),
+        ((2.0,), 0.0:20.0),
+        ((10.0,), 0.0:20.0),
+        ((10,), 0:20),
+    ])
+
+    rmathcomp_tests("tdist", [
+        ((1,), -5.0:0.1:5.0),
+        ((2,), -5.0:0.1:5.0),
+        ((5,), -5.0:0.1:5.0),
+        ((5,), -5:5),
+    ])
+
+    rmathcomp_tests("srdist", [
+        ((1,2), (0.0:0.2:5.0)),
+        ((2,2), (0.0:0.2:5.0)),
+        ((5,3), (0.0:0.2:5.0)),
+        ((10,2), (0.0:0.2:5.0)),
+        ((10,5), (0.0:0.2:5.0))
+    ])
+
+    # Note: Convergence fails in srdist with cdf values below 0.16 with df = 10, k = 5.
+    # Reduced df or k allows convergence. This test documents this behavior.
+    x = 0.15
+    q = srdistcdf(10, 5, 0.15)
+    rx = srdistinvcdf(10, 5, q)
+    rtol = 100eps(1.0)
+    @test_broken x ≈ rx atol=rtol rtol=rtol nans=true
+end