Fix stirlings2 bug

README.md

@@ -32,7 +32,7 @@ This library provides the following functions:
  - `multiexponents(m,n)`: returns the exponents in the multinomial expansion (x₁ + x₂ + ... + xₘ)ⁿ;
  - `primorial(n)`: returns the product of all positive prime numbers <= n; always returns a `BigInt`;
  - `stirlings1(n, k, signed=false)`: returns the `(n,k)`-th Stirling number of the first kind; the number is signed if `signed` is true; returns a `BigInt` only if given a `BigInt`.
- - `stirlings2(n, k)`: returns the `(n,k)`-th Stirling number of the second kind; returns a `BigInt` only if given a `BigInt`.
+ - `stirlings2(n, k)`: returns the `(n,k)`-th Stirling number of the second kind; may returns a `BigInt`, always when `3 <= k <= n-2`.
  - `nthperm(a, k)`: Compute the `k`th lexicographic permutation of the vector `a`.
  - `permutations(a)`: Generate all permutations of an indexable object `a` in lexicographic order.
 

src/numbers.jl

@@ -151,18 +151,23 @@ function stirlings1(n::Int, k::Int, signed::Bool=false)
     return (n - 1) * stirlings1(n - 1, k) + stirlings1(n - 1, k - 1)
 end
 
-function stirlings2(n::Int, k::Int)
-    if n < 0
-        throw(DomainError(n, "n must be nonnegative"))
-    elseif n == k == 0
-        return 1
-    elseif n == 0 || k == 0
-        return 0
-    elseif k == n - 1
-        return binomial(n, 2)
-    elseif k == 2
-        return 2^(n-1) - 1
+function stirlings2(n::Integer, k::Integer)
+    n >= 0 || throw(DomainError("n"))
+
+    if      n == k == 0;            one(n)
+    elseif  n == 0 || k == 0;       zero(n)
+    elseif  k == 1;                 one(n)
+    elseif  k == 2;                 2^(n-1)-1
+    elseif  k == n-1;               binomial(n, 2)
+    elseif  k == n;                 one(n)
+    else
+        # note: summing on itr leads to silent overflow without bigint
+        (bn, bk) = (big(n), big(k))
+        div(
+            sum(
+                ((-1)^(bk-i) * binomial(bk,i) * i^n
+                for i in 0:bk)),
+            factorial(bk))
     end
-
-    return k * stirlings2(n - 1, k) + stirlings2(n - 1, k - 1)
 end
+

src/partitions.jl

@@ -250,7 +250,7 @@ struct FixedSetPartitions{T<:AbstractVector}
     m::Int
 end
 
-Base.length(p::FixedSetPartitions) = nfixedsetpartitions(length(p.s),p.m)
+Base.length(p::FixedSetPartitions) = stirlings2(length(p.s),p.m)
 Base.eltype(p::FixedSetPartitions) = Vector{Vector{eltype(p.s)}}
 
 """
@@ -335,15 +335,6 @@ function nextfixedsetpartition(s::AbstractVector, m, a, b, n)
     return (part, (a,b,n))
 end
 
-function nfixedsetpartitions(n::Int, m::Int)
-    numpart = 0
-    for k = 0:m
-        numpart += (-1)^(m-k) * binomial(m, k) * (k^n)
-    end
-    numpart = div(numpart, factorial(m))
-    return numpart
-end
-
 # TODO: Base.DSP is no longer a thing in Julia 0.7
 #This function is still defined in Base because it is being used by Base.DSP
 #"""

test/numbers.jl

@@ -64,6 +64,7 @@
 @test stirlings2(6, 4) == 65
 @test stirlings2(6, 5) == 15
 @test stirlings2(6, 6) == 1
+@test stirlings2(25, 7) == 227832482998716310
 
 # bell
 @test bellnum.(0:10) == [