diff --git a/docs/src/functions_list.md b/docs/src/functions_list.md
index dbb81548..7afaf2a8 100644
--- a/docs/src/functions_list.md
+++ b/docs/src/functions_list.md
@@ -63,4 +63,5 @@ SpecialFunctions.beta
 SpecialFunctions.logbeta
 SpecialFunctions.logabsbeta
 SpecialFunctions.logabsbinomial
+SpecialFunctions.logmultinomial
 ```
diff --git a/docs/src/functions_overview.md b/docs/src/functions_overview.md
index cb85eb7d..146f0b77 100644
--- a/docs/src/functions_overview.md
+++ b/docs/src/functions_overview.md
@@ -17,6 +17,7 @@ Here the *Special Functions* are listed according to the structure of [NIST Digi
 | [`logabsgamma(x)`](@ref SpecialFunctions.logabsgamma)                                           | accurate `log(abs(gamma(x)))` for large `x`                                                                                  |
 | [`lgamma(x)`](@ref SpecialFunctions.lgamma)                                           | accurate `log(gamma(x))` for large `x`                                                                                                                          |
 | [`logfactorial(x)`](@ref SpecialFunctions.logfactorial)                                            | accurate `log(factorial(x))` for large `x`; same as `lgamma(x+1)` for `x > 1`, zero otherwise                                                                   |
+| [`logmultinomial(k...)`](@ref SpecialFunctions.logmultinomial)                                            | accurate `log(multinomial(k...))` for large multisets                                                                   |
 | [`beta(x,y)`](@ref SpecialFunctions.beta)                                           | [beta function](https://en.wikipedia.org/wiki/Beta_function) at `x,y`                                                                                           |
 | [`logbeta(x,y)`](@ref SpecialFunctions.logbeta)                                          | accurate `log(beta(x,y))` for large `x` or `y`      |
 | [`logabsbeta(x,y)`](@ref SpecialFunctions.logabsbeta)                                          | accurate `log(abs(beta(x,y)))` for large `x` or `y`     |
diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl
index 814e2a46..d718cc3c 100644
--- a/src/SpecialFunctions.jl
+++ b/src/SpecialFunctions.jl
@@ -57,7 +57,8 @@ export
     zeta,
     sinint,
     cosint,
-    lbinomial
+    lbinomial,
+    logmultinomial
 
 include("bessel.jl")
 include("erf.jl")
diff --git a/src/gamma.jl b/src/gamma.jl
index 8273bf31..1c8c8cbe 100644
--- a/src/gamma.jl
+++ b/src/gamma.jl
@@ -566,7 +566,7 @@ function polygamma(m::Integer, z::Number)
     polygamma(m, x)
 end
 
-export gamma, loggamma, logabsgamma, beta, logbeta, logabsbeta, logfactorial, logabsbinomial
+export gamma, loggamma, logabsgamma, beta, logbeta, logabsbeta, logfactorial, logabsbinomial, logmultinomial
 
 ## from base/special/gamma.jl
 
@@ -922,3 +922,25 @@ function logabsbinomial(n::T, k::T) where {T<:Integer}
     end
 end
 logabsbinomial(n::Integer, k::Integer) = logabsbinomial(promote(n, k)...)
+
+#=
+logmultinomial computes the log of the multinomial coefficient.
+From Wolfram Mathworld (https://mathworld.wolfram.com/MultinomialCoefficient.html):
+
+The multinomial coefficients
+
+$$(n_1, n_2, \ldots, n_k)! = \frac{(n_1+n_2+...+n_k)!}{n_1! n_2! \ldots n_k!}$$
+
+are the terms in the multinomial series expansion. In other words, the number
+of distinct permutations in a multiset of $k$ distinct elements of multiplicity
+$n_i$ $(1 \le i \le k)$ is $(n_1, \ldots, n_k)!$.
+=#
+function logmultinomial(multiset...)
+    numerator = 0
+    denominator = 0.0
+    @inbounds for multiplicity in multiset
+        numerator += multiplicity
+        denominator += loggamma(multiplicity + 1)
+    end
+    return loggamma(numerator + 1) - denominator
+end