From 55725ed6a0fde396137d340600aba63c0b6c6cda Mon Sep 17 00:00:00 2001
From: Moritz Schauer <moritzschauer@web.de>
Date: Fri, 11 Nov 2022 09:53:07 +0100
Subject: [PATCH] Reviving #84

> The documentation for stirlings1 and stirlings2 indicates that they can use BigInts if provided with such as argument. However the existing functions specify n and k arguments to be Int. This patch allows support for BigInt integers without breaking type consistency of return types, and also provides for some memoization, since otherwise the functions, since they are recursive, can be slow when handling very large integers.

(credit https://github.com/wherrera10)
---
 src/numbers.jl | 34 +++++++++++++++++++++++-----------
 1 file changed, 23 insertions(+), 11 deletions(-)

diff --git a/src/numbers.jl b/src/numbers.jl
index 4039f9f..60e989a 100644
--- a/src/numbers.jl
+++ b/src/numbers.jl
@@ -68,7 +68,7 @@ end
 """
     narayana(n,k)
 
-Compute the Narayana number `N(n,k)` given by ``\\frac{1}{n}\\binom{n}{k}\\binom{n}{k-1}``
+Compute the Narayana number `N(n,k)`` given by ``\\frac{1}{n}\\binom{n}{k}\\binom{n}{k-1}``
 Wikipedia : https://en.wikipedia.org/wiki/Narayana_number
 """
 function narayana(bn::Integer,bk::Integer)
@@ -125,19 +125,23 @@ function lucasnum(n::Integer)
     return z[]
 end
 
-function stirlings1(n::Int, k::Int, signed::Bool=false)
+stirlings1cache = Dict()
+
+function stirlings1(n::Integer, k::Integer, signed::Bool=false)
     if signed == true
         return (-1)^(n - k) * stirlings1(n, k)
     end
 
-    if n < 0
+    if haskey(stirlings1cache, Pair(n, k))
+        return stirlings1cache[Pair(n, k)]
+    elseif n < 0
         throw(DomainError(n, "n must be nonnegative"))
     elseif n == k == 0
-        return 1
+        return one(n)
     elseif n == 0 || k == 0
-        return 0
+        return zero(n)
     elseif n == k
-        return 1
+        return one(n)
     elseif k == 1
         return factorial(n-1)
     elseif k == n - 1
@@ -148,21 +152,29 @@ function stirlings1(n::Int, k::Int, signed::Bool=false)
         return binomial(n, 2) * binomial(n, 4)
     end
 
-    return (n - 1) * stirlings1(n - 1, k) + stirlings1(n - 1, k - 1)
+    ret = (n - 1) * stirlings1(n - 1, k) + stirlings1(n - 1, k - 1)
+    stirlings1cache[Pair(n, k)] = ret
+    return ret
 end
 
-function stirlings2(n::Int, k::Int)
+stirlings2cache = Dict()
+
+function stirlings2(n::Integer, k::Integer)
     if n < 0
         throw(DomainError(n, "n must be nonnegative"))
+    elseif haskey(stirlings2cache, Pair(n, k))
+        return stirlings2cache[Pair(n, k)]
     elseif n == k == 0
-        return 1
+        return one(n)
     elseif n == 0 || k == 0
-        return 0
+        return zero(n)
     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)
+    ret = k * stirlings2(n - 1, k) + stirlings2(n - 1, k - 1)
+    stirlings2cache[Pair(n, k)] = ret
+    ret
 end