Add Lah numbers

Author: kwankyu

src/sage/combinat/combinat.py

@@ -49,15 +49,6 @@
    and is represented by a sorted list of length k containing elements
    from S.
 
-.. WARNING::
-
-   The following function is deprecated and will soon be removed.
-
-    - Permutations of a multiset, :func:`permutations`,
-      :func:`permutations_iterator`, :func:`number_of_permutations`. A
-      permutation is a list that contains exactly the same elements but possibly
-      in different order.
-
 **Related functions:**
 
 -  Bernoulli polynomials, :func:`bernoulli_polynomial`
@@ -95,8 +86,7 @@
 The ``sage.groups.perm_gps.permgroup_elements``
 contains the following combinatorial functions:
 
-
--  matrix method of PermutationGroupElement yielding the
+-  matrix method of :class:`PermutationGroupElement` yielding the
    permutation matrix of the group element.
 
 .. TODO::
@@ -114,7 +104,7 @@
 
 AUTHORS:
 
-- David Joyner (2006-07): initial implementation.
+- David Joyner (2006-07): initial implementation
 
 - William Stein (2006-07): editing of docs and code; many
   optimizations, refinements, and bug fixes in corner cases
@@ -132,8 +122,7 @@
 - Punarbasu Purkayastha (2012-12): deprecate arrangements, combinations,
   combinations_iterator, and clean up very old deprecated methods.
 
-Functions and classes
----------------------
+- Kwankyu Lee (2025-01): added Lah numbers
 """
 
 # ****************************************************************************
@@ -505,7 +494,7 @@ def narayana_number(n: Integer, k: Integer) -> Integer:
 
     - ``n`` -- integer
 
-    - ``k`` -- integer between ``0`` and ``n - 1``
+    - ``k`` -- integer between `0` and `n - 1`
 
     OUTPUT: integer
 
@@ -890,6 +879,10 @@ def stirling_number1(n, k, algorithm='gap') -> Integer:
     """
     n = ZZ(n)
     k = ZZ(k)
+    if k > n:
+        return ZZ.zero()
+    if k == 0:
+        return ZZ.zero() if n else ZZ.one()
     if algorithm == 'gap':
         from sage.libs.gap.libgap import libgap
         return libgap.Stirling1(n, k).sage()
@@ -1016,6 +1009,10 @@ def stirling_number2(n, k, algorithm=None) -> Integer:
     """
     n = ZZ(n)
     k = ZZ(k)
+    if k > n:
+        return ZZ.zero()
+    if k == 0:
+        return ZZ.zero() if n else ZZ.one()
     if algorithm is None:
         return _stirling_number2(n, k)
     if algorithm == 'gap':
@@ -1029,6 +1026,44 @@ def stirling_number2(n, k, algorithm=None) -> Integer:
     raise ValueError("unknown algorithm: %s" % algorithm)
 
 
+def lah_number(n, k) -> Integer:
+    r"""
+    Return the Lah number `L(n,k)`
+
+    This is the number of ways to partition a set of `n` elements into `k`
+    pairwise disjoint nonempty linearly-ordered subsets.
+
+    This is also called the Stirling number of the third kind.
+
+    See :wikipedia:`Lah_number`.
+
+    INPUT:
+
+    - ``n`` -- nonnegative integer
+    - ``k`` -- nonnegative integer
+
+    EXAMPLES:
+
+        sage: lah_number(50, 30)
+        3242322638238907670866645288893161825894400000
+
+    We verify a well-known identity::
+
+        sage: S1 = stirling_number1; S2 = stirling_number2
+        sage: all(lah_number(n, k) == sum(S1(n, j) * S2(j, k) for j in [k..n])
+        ....:     for n in range(10) for k in range(10))
+        True
+    """
+    n = ZZ(n)
+    k = ZZ(k)
+    if k > n:
+        return ZZ.zero()
+    if k == 0:
+        return ZZ.zero() if n else ZZ.one()
+    a = n.factorial() // k.factorial()
+    return a**2 * k // (n * (n - k).factorial())
+
+
 def polygonal_number(s, n):
     r"""
     Return the `n`-th `s`-gonal number.