From 85fc819b020e9f373c070c17de14591fe1b1b5dc Mon Sep 17 00:00:00 2001 From: Olivier MATHIEU Date: Tue, 14 Apr 2020 03:13:47 +0200 Subject: [PATCH 1/6] Update numbers.jl --- src/numbers.jl | 32 +++++++++++++++++++------------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/src/numbers.jl b/src/numbers.jl index 4465c77..91e4ac5 100644 --- a/src/numbers.jl +++ b/src/numbers.jl @@ -151,18 +151,24 @@ 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")) + k in 0:n || throw(DomainError("k")) + + if n == k == 0; one(n) + elseif n == 0 || k == 0; zero(n) + elseif k == 1; 1 + elseif k == 2; 2^(n-1)-1 + elseif k == n-1; binomial(n, 2) + elseif k == n; 1 + 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 + From 59c1bcdba646fcbf3318327b5a2765d0c0fb6a49 Mon Sep 17 00:00:00 2001 From: Olivier MATHIEU Date: Tue, 14 Apr 2020 03:14:48 +0200 Subject: [PATCH 2/6] Update partitions.jl --- src/partitions.jl | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/src/partitions.jl b/src/partitions.jl index 9ba2969..2f6e17e 100644 --- a/src/partitions.jl +++ b/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 #""" From 1f5b291c01c82cf9aaea72e4ec52cd5d035d07e1 Mon Sep 17 00:00:00 2001 From: Olivier MATHIEU Date: Tue, 14 Apr 2020 03:19:40 +0200 Subject: [PATCH 3/6] Update numbers.jl --- src/numbers.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/numbers.jl b/src/numbers.jl index 91e4ac5..7961250 100644 --- a/src/numbers.jl +++ b/src/numbers.jl @@ -153,14 +153,13 @@ end function stirlings2(n::Integer, k::Integer) n >= 0 || throw(DomainError("n")) - k in 0:n || throw(DomainError("k")) if n == k == 0; one(n) elseif n == 0 || k == 0; zero(n) - elseif k == 1; 1 + elseif k == 1; zero(n) elseif k == 2; 2^(n-1)-1 elseif k == n-1; binomial(n, 2) - elseif k == n; 1 + elseif k == n; zero(n) else # note: summing on itr leads to silent overflow without bigint (bn, bk) = (big(n), big(k)) From 07494def64d7f7e9592ba943c630a1c0eb9c76fc Mon Sep 17 00:00:00 2001 From: Olivier MATHIEU Date: Tue, 14 Apr 2020 03:21:27 +0200 Subject: [PATCH 4/6] Update numbers.jl --- src/numbers.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/numbers.jl b/src/numbers.jl index 7961250..7bb0033 100644 --- a/src/numbers.jl +++ b/src/numbers.jl @@ -156,10 +156,10 @@ function stirlings2(n::Integer, k::Integer) if n == k == 0; one(n) elseif n == 0 || k == 0; zero(n) - elseif k == 1; zero(n) + elseif k == 1; one(n) elseif k == 2; 2^(n-1)-1 elseif k == n-1; binomial(n, 2) - elseif k == n; zero(n) + elseif k == n; one(n) else # note: summing on itr leads to silent overflow without bigint (bn, bk) = (big(n), big(k)) From 917d4c14d9ae090f3d86157cb1cb1a4dedb6efb1 Mon Sep 17 00:00:00 2001 From: Olivier MATHIEU Date: Tue, 14 Apr 2020 03:26:52 +0200 Subject: [PATCH 5/6] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index e8b97d6..0d549ce 100644 --- a/README.md +++ b/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. From d4d32d930fb8855bf4923704e9b2ea64caff7cb3 Mon Sep 17 00:00:00 2001 From: Olivier MATHIEU Date: Tue, 14 Apr 2020 03:29:57 +0200 Subject: [PATCH 6/6] Update numbers.jl --- test/numbers.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/numbers.jl b/test/numbers.jl index b729d36..ebabe09 100644 --- a/test/numbers.jl +++ b/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) == [