Repository: sympy/sympy
Pull request: Functions: Combinatorial-DP based implementation of Stirling Numbers
Fixes #11650

Recursive implementation was slower as it was getting executed in exponential time hence maximum recursion depth got reached, when both kinds of stirling numbers were called.

**[stirling(200, i, kind=1) for i in range(201)]
RuntimeError: maximum recursion depth exceeded while calling a Python object**

Similarly,

**[stirling(200, i, kind=2) for i in range(201)]
RuntimeError: maximum recursion depth exceeded while calling a Python object**

Now DP based implementation takes O(n*k) time and is also a space optimized solution and **we don't get the above Runtime Error**.

<!-- BEGIN RELEASE NOTES -->
* functions
  * Fixed maximum recursion errors and improved performance for `stirling`. 

<!-- END RELEASE NOTES -->


Problem statement:
Slow Recursive implementation of Stirling Numbers
Hi folks,

I noticed the other day that the current implementation for Stirling Numbers is recursive and thus it works very slowly. I was wondering if there is any idea behind it why it should be recursive.

If there is no specific reason then I would like to work on this to update the implementation to a dynamic one.

Hints from issue discussion:
Same question here. May I also ask that what is the reason for the recursive implementation? I wrote a dynamic programming version but the tests only cover small n and k so the there's no big difference in running time...

Why is this issue open? Is there anything which still needs to be addressed? 
Is Dynamic programming solution needs to be be implemented?
@oscarbenjamin 
I haven't looked at this at all myself. If you are interested then you can check whether the implementation is recursive or slow.
How to see whether it is slow or not?
Though I can see the code is indeed a recursive implementation and therefore will definitely be slow for large numbers.


There is nothing inherently slow about recursion although we will hit the recursion limit of we recurse too many times. Branched recursion can be very slow though - the classic example of this is computing numbers in the fibonacci sequence recursively using `fib(n) = fib(n-1) + fib(n-2)` where the number of calls to `fib` grows exponentially in `n`. That is an example where "dynamic programming" can make things more efficient. As I said though I haven't looked at the code for the Stirling numbers in SymPy so I don't know whether that is what happens.
Yes it is exceeding maximum recursion depth.

[stirling(200, i, kind=1) for i in range(201)]
RuntimeError: maximum recursion depth exceeded while calling a Python object
I have a DP space optimized solution which is working fine for numbers >=200 as well.
Recursive implementation takes exponential time, but DP implementation will take O(n*k) time.
Shall I go ahead and implement it?
If it is an improvement then go ahead
Yeah, it is an improvement, just please review it and merge it after the build.

Commit message: Add back some special formula

Patch to review:
diff --git a/sympy/functions/combinatorial/numbers.py b/sympy/functions/combinatorial/numbers.py
index 1693827d5e8a..bf2f8938379f 100644
--- a/sympy/functions/combinatorial/numbers.py
+++ b/sympy/functions/combinatorial/numbers.py
@@ -1715,22 +1715,25 @@ def _stirling1(n, k):
         return S.One
     if 0 in (n, k):
         return S.Zero
-    n1 = n - 1
 
     # some special values
     if n == k:
         return S.One
-    elif k == 1:
-        return factorial(n1)
-    elif k == n1:
+    elif k == n - 1:
         return binomial(n, 2)
     elif k == n - 2:
         return (3*n - 1)*binomial(n, 3)/4
     elif k == n - 3:
         return binomial(n, 2)*binomial(n, 4)
 
-    # general recurrence
-    return n1*_stirling1(n1, k) + _stirling1(n1, k - 1)
+    row = [1]+[0 for _ in range(k)]
+    for i in range(1, n+1):
+        new = [0]
+        for j in range(1, k+1):
+            stirling = (i-1) * row[j] + row[j-1]
+            new.append(stirling)
+        row = new
+    return Integer(row[k])
 
 
 @cacheit
@@ -1739,16 +1742,21 @@ def _stirling2(n, k):
         return S.One
     if 0 in (n, k):
         return S.Zero
-    n1 = n - 1
 
     # some special values
-    if k == n1:
+    if k == n - 1:
         return binomial(n, 2)
     elif k == 2:
-        return 2**n1 - 1
-
-    # general recurrence
-    return k*_stirling2(n1, k) + _stirling2(n1, k - 1)
+        return Integer(2**(n - 1) - 1)
+
+    row = [1]+[0 for _ in range(k)]
+    for i in range(1, n+1):
+        new = [0]
+        for j in range(1, k+1):
+            stirling = j * row[j] + row[j-1]
+            new.append(stirling)
+        row = new
+    return Integer(row[k])
 
 
 def stirling(n, k, d=None, kind=2, signed=False):
diff --git a/sympy/functions/combinatorial/tests/test_comb_numbers.py b/sympy/functions/combinatorial/tests/test_comb_numbers.py
index 0ef8d6405e5c..e5878e78ddd0 100644
--- a/sympy/functions/combinatorial/tests/test_comb_numbers.py
+++ b/sympy/functions/combinatorial/tests/test_comb_numbers.py
@@ -12,7 +12,7 @@
 
 from sympy.core.compatibility import range
 from sympy.core.expr import unchanged
-from sympy.core.numbers import GoldenRatio
+from sympy.core.numbers import GoldenRatio, Integer
 
 from sympy.utilities.pytest import XFAIL, raises
 
@@ -473,7 +473,8 @@ def test_nC_nP_nT():
         multiset_permutations, multiset_combinations, multiset_partitions,
         partitions, subsets, permutations)
     from sympy.functions.combinatorial.numbers import (
-        nP, nC, nT, stirling, _multiset_histogram, _AOP_product)
+        nP, nC, nT, stirling, _stirling1, _stirling2, _multiset_histogram,
+        _AOP_product)
     from sympy.combinatorics.permutations import Permutation
     from sympy.core.numbers import oo
     from random import choice
@@ -596,6 +597,10 @@ def test_nC_nP_nT():
     assert stirling(3, 4, kind=1) == stirling(3, 4, kind=1) == 0
     raises(ValueError, lambda: stirling(-2, 2))
 
+    # Assertion that the return type is SymPy Integer.
+    assert isinstance(_stirling1(6, 3), Integer)
+    assert isinstance(_stirling2(6, 3), Integer)
+
     def delta(p):
         if len(p) == 1:
             return oo


Code review comments:
diff_hunk: @@ -1715,22 +1715,25 @@ def _stirling1(n, k):
         return S.One
     if 0 in (n, k):
         return S.Zero
-    n1 = n - 1
 
     # some special values
     if n == k:
         return S.One
-    elif k == 1:
-        return factorial(n1)
-    elif k == n1:
+    elif k == n - 1:
         return binomial(n, 2)
     elif k == n - 2:
         return (3*n - 1)*binomial(n, 3)/4
     elif k == n - 3:
         return binomial(n, 2)*binomial(n, 4)
 
-    # general recurrence
-    return n1*_stirling1(n1, k) + _stirling1(n1, k - 1)
+    row = [1]+[0 for _ in range(k)]
+    for i in range(1, n+1):
+        new = [0]
+        for j in range(1, k+1):
+            stirling = (i-1) * row[j] + row[j-1]
+            new.append(stirling)
path: sympy/functions/combinatorial/numbers.py
text: @user1:
Is it possible to reuse `row` (instead of appending to `new`) if `j` would run from `k` down to `1`?

@author:
I don't think so, as far as I can understand, you are suggesting a Top-Down approach, where we will try to reuse row, but I don't know whether the same relation will hold or not.

@user1:
I suggest `row[j] = (i-1) * row[j] + row[j-1]`.

@author:
> I suggest `row[j] = (i-1) * row[j] + row[j-1]`.

Some tests are failing after doing this change. 

@user1:
>  if j would run from k down to 1

Together with this change?

@author:
> > if j would run from k down to 1
> 
> Together with this change?

for j in range(k+1, 1) ??

@user1:
`for j in range(k, 0, -1)` should do that.

@author:
Yes they are failing. Answer is not coming right, and it won't because we want to generate from bottom to up, so j must run down from 1 o k

@author:
Please review the pull request @user2 @user3 @user1 
If any changes are still required, please suggest, otherwise I think the request is ready to be merged.

@user1:
The tests succeed on my system with this diff:
```
--- a/sympy/functions/combinatorial/numbers.py
+++ b/sympy/functions/combinatorial/numbers.py
@@ -1726,13 +1726,10 @@ def _stirling1(n, k):
     elif k == n - 3:
         return binomial(n, 2)*binomial(n, 4)
 
-    row = [1]+[0 for _ in range(k)]
-    for i in range(1, n+1):
-        new = [0]
-        for j in range(1, k+1):
-            stirling = (i-1) * row[j] + row[j-1]
-            new.append(stirling)
-        row = new
+    row = [0, 1] + [0]*(k-1)  # for n = 1 
+    for i in range(2, n+1):
+        for j in range(k, 0, -1):
+            row[j] = (i-1) * row[j] + row[j-1]
     return Integer(row[k])
```
The initial value `row[0] = 1` is only needed for computing the case `i = 1`. It is necessary to replace `row[0]` by 0 afterwards as is done here: `new = [0]`. This could be set explicitly after the `j`-loop. However, I would prefer starting from the case `n = 1` where `row[0] = 0` is already set. (The case `n = 0` has been handled above.)

@author:
That's true, yeah these are small optimizations which should have been made. Extremely sorry for the inconvenience.

@user1:
Another optimization is possible. For small values of `i`, `row` will mostly consist of zeros. It is not necessary to reassign them for `j > i`. So it seems that the `j` loop could be `for j in range(min(k, i), 0, -1)`.

@author:
Any further improvements after this one? @user1 @user3 @user2 
If there are, please suggest. Otherwise, I will commit all the changes once and for all, so that the request gets ready to be merged.

diff_hunk: @@ -1715,22 +1715,25 @@ def _stirling1(n, k):
         return S.One
     if 0 in (n, k):
         return S.Zero
-    n1 = n - 1
 
     # some special values
     if n == k:
         return S.One
-    elif k == 1:
-        return factorial(n1)
-    elif k == n1:
+    elif k == n - 1:
         return binomial(n, 2)
     elif k == n - 2:
         return (3*n - 1)*binomial(n, 3)/4
     elif k == n - 3:
         return binomial(n, 2)*binomial(n, 4)
 
-    # general recurrence
-    return n1*_stirling1(n1, k) + _stirling1(n1, k - 1)
+    row = [1]+[0 for _ in range(k)]
path: sympy/functions/combinatorial/numbers.py
text: @user1:
I think that `[0]*k` could be used.

@author:
for row declaration? Would it make a significant difference?