Merge pull request #16 from ymer/master

New implementation for stirlings

Author: andreasnoack

README.md

@@ -25,7 +25,8 @@ This library provides the following functions:
  - `multifactorial(n)`: returns the m-multifactorial n(!^m); always returns a `BigInt`;
  - `multinomial(k...)`: receives a tuple of `k_1, ..., k_n` and calculates the multinomial coefficient `(n k)`, where `n = sum(k)`; returns a `BigInt` only if given a `BigInt`;
  - `primorial(n)`: returns the product of all positive prime numbers <= n; always returns a `BigInt`;
- - `stirlings1(n, k)`: the signed `(n,k)`-th Stirling number of the first kind; returns a `BigInt` only if given 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`.
 
 Young diagrams
 --------------

src/numbers.jl

@@ -7,8 +7,11 @@ export bellnum,
     lassallenum,
     legendresymbol,
     lucasnum,
-    stirlings1
+    stirlings1,
+    stirlings2
 
+import Base: factorial, binomial
+	
 "Returns the n-th Bell number"
 function bellnum(bn::Integer)
     if bn < 0
@@ -85,9 +88,44 @@ function lucasnum(n::Integer)
     return z
 end
 
-"Returns s(n, k), the signed Stirling number of first kind"
-function stirlings1(n::Integer, k::Integer)
-    p = poly(0:(n-1))
-    p[n - k + 1]
+function stirlings1(n::Int, k::Int, signed::Bool=false)
+    if signed == true
+        return (-1)^(n - k) * stirlings1(n, k)
+    end
+    
+    if n < 0
+        throw(DomainError())
+    elseif n == k == 0
+        return 1
+    elseif n == 0 || k == 0
+        return 0
+    elseif n == k
+        return 1
+    elseif k == 1
+        return factorial(n-1)
+    elseif k == n - 1
+        return binomial(n, 2)
+    elseif k == n - 2
+        return div((3 * n - 1) * binomial(n, 3), 4)
+    elseif k == n - 3
+        return binomial(n, 2) * binomial(n, 4)
+    end
+    
+    return (n - 1) * stirlings1(n - 1, k) + stirlings1(n - 1, k - 1)
 end
 
+function stirlings2(n::Int, k::Int)
+    if n < 0
+        throw(DomainError())
+    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
+    end
+    
+    return k * stirlings2(n - 1, k) + stirlings2(n - 1, k - 1)
+end

test/numbers.jl

@@ -23,11 +23,42 @@ using Base.Test
 @test_throws DomainError lucasnum(-1)
 
 # stirlings1
-@test sum([abs(stirlings1(8, i)) for i = 0:8]) == factorial(8)
+@test_throws DomainError stirlings1(-1, 1)
+@test typeof(stirlings1(6, 2)) == Int
+@test stirlings1(0, 0) == 1
+@test stirlings1(1, 1) == 1
+@test stirlings1(2, 6) == 0
+@test stirlings1(6, 0) == 0
+@test stirlings1(6, 0, true) == 0
+@test stirlings1(6, 1) == 120
+@test stirlings1(6, 1, true) == -120
+@test stirlings1(6, 2) == 274
+@test stirlings1(6, 2, true) == 274
+@test stirlings1(6, 3) == 225
+@test stirlings1(6, 3, true) == -225
+@test stirlings1(6, 4) == 85
+@test stirlings1(6, 4, true) == 85
+@test stirlings1(6, 5) == 15
+@test stirlings1(6, 5, true) == -15
+@test stirlings1(6, 6) == 1
+@test stirlings1(6, 6, true) == 1
+@test sum([abs(stirlings1(8, i, true)) for i = 0:8]) == factorial(8)
+
+# stirlings2
+@test_throws DomainError stirlings2(-1, 1)
+@test typeof(stirlings2(6, 2)) == Int
+@test stirlings2(0, 0) == 1
+@test stirlings2(1, 1) == 1
+@test stirlings2(2, 6) == 0
+@test stirlings2(6, 0) == 0
+@test stirlings2(6, 1) == 1
+@test stirlings2(6, 2) == 31
+@test stirlings2(6, 3) == 90
+@test stirlings2(6, 4) == 65
+@test stirlings2(6, 5) == 15
+@test stirlings2(6, 6) == 1
 
 # bell
 @test bellnum(7) == 877
 @test bellnum(42) == parse(BigInt,"35742549198872617291353508656626642567")
 @test_throws DomainError bellnum(-1)
-
-