88 doublefactorial,
99 hyperfactorial,
1010 multifactorial,
11- gamma,
1211 primorial,
1312 multinomial
1413
@@ -29,6 +28,18 @@ function Base.factorial(n::T, k::T) where T<:Integer
2928 end
3029 return f
3130end
31+ function Base. factorial (n:: BigInt , k:: BigInt )
32+ if k < 0 || n < 0 || k > n
33+ throw (DomainError ((n, k), " n and k must be nonnegative with k ≤ n" ))
34+ end
35+ f = BigInt (1 )
36+ n = deepcopy (n) # avoid mutating input
37+ while n > k
38+ Base. GMP. MPZ. mul! (f, n)
39+ Base. GMP. MPZ. sub_ui! (n, 1 )
40+ end
41+ return f
42+ end
3243Base. factorial (n:: Integer , k:: Integer ) = factorial (promote (n, k)... )
3344
3445
@@ -38,9 +49,23 @@ Base.factorial(n::Integer, k::Integer) = factorial(promote(n, k)...)
3849Compute the number of permutations of `n` with no fixed points, also known as the
3950subfactorial. An alias `subfactorial` for this function is provided for convenience.
4051"""
41- function derangement (sn:: Integer )
42- n = BigInt (sn)
43- return numerator (factorial (n) * sum ([(- 1 )^ k // factorial (k) for k = 0 : n]))
52+ function derangement (n:: Integer )
53+ if n < 0
54+ throw (DomainError (n, " n must be nonnegative" ))
55+ elseif n <= 1
56+ return BigInt (1 - n)
57+ end
58+ d = BigInt (0 )
59+ for i in 2 : n
60+ # d = i * d + (iseven(i) ? 1 : -1)
61+ Base. GMP. MPZ. mul_ui! (d, i)
62+ if iseven (i)
63+ Base. GMP. MPZ. add_ui! (d, 1 )
64+ else
65+ Base. GMP. MPZ. sub_ui! (d, 1 )
66+ end
67+ end
68+ return d
4469end
4570const subfactorial = derangement
4671
95120"""
96121 multinomial(k...)
97122
98- Multinomial coefficient where `n = sum(k)`.
123+ Compute the multinomial coefficient
124+ ``\\ binom{n}{k_1,k_2,...,k_i} = \\ frac{n!}{k_1!k_2! \\ cdots k_i!}, n = \\ sum{k_i}``.
125+ Throws an `OverflowError` when the input is too large.
126+
127+ See Also: `binomial`.
128+
129+ # Examples
130+ ```jldoctest
131+ julia> # (x+y)^2 = x^2 + 2xy + y^2
132+
133+ julia> multinomial(2, 0)
134+ 1
135+
136+ julia> multinomial(1, 1)
137+ 2
138+
139+ julia> multinomial(0, 2)
140+ 1
141+
142+ julia> multinomial(10, 10, 10, 10)
143+ ERROR: OverflowError: 5550996791340 * 847660528 overflowed for type Int64
144+ Stacktrace:
145+ [...]
146+ ```
147+
148+ # External links
149+ - [Definitions](https://dlmf.nist.gov/26.4.2) on DLMF
150+ - [Multinomial theorem](https://en.wikipedia.org/wiki/Multinomial_theorem) on Wikipedia
99151"""
100152function multinomial (k... )
101153 s = 0
102154 result = 1
103155 @inbounds for i in k
104156 s += i
105- result *= binomial (s, i)
157+ bi = binomial (s, i)
158+ result = Base. Checked. checked_mul (result, bi)
106159 end
107160 result
108161end
0 commit comments