-
Notifications
You must be signed in to change notification settings - Fork 59
/
factorials.jl
108 lines (93 loc) · 2.41 KB
/
factorials.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
#Factorials and elementary coefficients
export
derangement,
partialderangement,
factorial,
subfactorial,
doublefactorial,
hyperfactorial,
multifactorial,
gamma,
primorial,
multinomial
# TODO: This should really live in Base, otherwise it's type piracy
"""
factorial(n, k)
Compute ``n!/k!``.
"""
function Base.factorial(n::T, k::T) where T<:Integer
if k < 0 || n < 0 || k > n
throw(DomainError((n, k), "n and k must be nonnegative with k ≤ n"))
end
f = one(T)
while n > k
f = Base.checked_mul(f, n)
n -= 1
end
return f
end
Base.factorial(n::Integer, k::Integer) = factorial(promote(n, k)...)
"""
derangement(n)
Compute the number of permutations of `n` with no fixed points, also known as the
subfactorial. An alias `subfactorial` for this function is provided for convenience.
"""
function derangement(sn::Integer)
n = BigInt(sn)
return numerator(factorial(n) * sum([(-1)^k // factorial(k) for k = 0:n]))
end
const subfactorial = derangement
function doublefactorial(n::Integer)
if n < 0
throw(DomainError(n, "n must be nonnegative"))
end
z = Ref{BigInt}(0)
ccall((:__gmpz_2fac_ui, :libgmp), Cvoid, (Ref{BigInt}, UInt), z, UInt(n))
return z[]
end
"""
partialderangement(n, k)
Compute the number of permutations of `n` with exactly k fixed points.
"""
function partialderangement(n::Integer, k::Integer)
if n < 0
throw(DomainError(n))
end
if k < 0 || k > n
throw(DomainError(k))
end
a = BigInt(n)
b = BigInt(k)
return subfactorial(a - b) * binomial(a, b)
end
# Hyperfactorial
hyperfactorial(n::Integer) = n==0 ? BigInt(1) : prod(i->i^i, BigInt(1):n)
function multifactorial(n::Integer, m::Integer)
if n < 0
throw(DomainError(n, "n must be nonnegative"))
end
z = Ref{BigInt}(0)
ccall((:__gmpz_mfac_uiui, :libgmp), Cvoid, (Ref{BigInt}, UInt, UInt), z, UInt(n), UInt(m))
return z[]
end
function primorial(n::Integer)
if n < 0
throw(DomainError(n, "n must be nonnegative"))
end
z = Ref{BigInt}(0)
ccall((:__gmpz_primorial_ui, :libgmp), Cvoid, (Ref{BigInt}, UInt), z, UInt(n))
return z[]
end
"""
multinomial(k...)
Multinomial coefficient where `n = sum(k)`.
"""
function multinomial(k...)
s = 0
result = 1
@inbounds for i in k
s += i
result *= binomial(s, i)
end
result
end