Check overflow in multinomial

src/factorials.jl

@@ -95,14 +95,42 @@ end
 """
     multinomial(k...)
 
-Multinomial coefficient where `n = sum(k)`.
+Compute the multinomial coefficient
+``\\binom{n}{k_1,k_2,...,k_i} = \\frac{n!}{k_1!k_2! \\cdots k_i!}, n = \\sum{k_i}``.
+Throws an `OverflowError` when the input is too large.
+
+See Also: `binomial`.
+
+# Examples
+```jldoctest
+julia> # (x+y)^2 = x^2 + 2xy + y^2
+
+julia> multinomial(2, 0)
+1
+
+julia> multinomial(1, 1)
+2
+
+julia> multinomial(0, 2)
+1
+
+julia> multinomial(10, 10, 10, 10)
+ERROR: OverflowError: 5550996791340 * 847660528 overflowed for type Int64
+Stacktrace:
+[...]
+```
+
+# External links
+- [Definitions](https://dlmf.nist.gov/26.4.2) on DLMF
+- [Multinomial theorem](https://en.wikipedia.org/wiki/Multinomial_theorem) on Wikipedia
 """
 function multinomial(k...)
     s = 0
     result = 1
     @inbounds for i in k
         s += i
-        result *= binomial(s, i)
+        bi = binomial(s, i)
+        result = Base.Checked.checked_mul(result, bi)
     end
     result
 end

test/factorials.jl

@@ -31,8 +31,22 @@
     @test multifactorial(40, 2) == doublefactorial(40)
     @test_throws DomainError multifactorial(-1, 1)
 
-    # multinomial
-    @test multinomial(1, 4, 4, 2) == 34650
+    @testset "multinomial" begin
+        # > For k=0,1, the multinomial coefficient is defined to be 1
+        # https://dlmf.nist.gov/26.4#i.p1
+        @test multinomial() == 1
+        @test multinomial(0) == 1
+
+        @test multinomial(1, 4, 4, 2) == 34650
+        # wolfram:  Multinomial[10, 10, 10, 5]
+        @test multinomial(10, 10, 10, 5) == 1_802_031_190_366_286_880
+
+        # checked_mul overflowed for type Int64
+        @test_throws OverflowError multinomial(10, 10, 10, 6)
+        @test_throws OverflowError multinomial(10, 10, 10, 10)
+        # binomial(200, 100) overflows
+        @test_throws OverflowError multinomial(100, 100)  
+    end
 
     # primorial
     @test primorial(17) == 510510