Trac #33733: allow to use flint for Stirling numbers

of both kind

URL: https://trac.sagemath.org/33733
Reported by: chapoton
Ticket author(s): Frédéric Chapoton
Reviewer(s): Marc Mezzarobba, Travis Scrimshaw

src/sage/combinat/combinat.py

@@ -364,9 +364,9 @@ def bell_number(n, algorithm='flint', **options) -> Integer:
 
     TESTS::
 
-        sage: all(bell_number(n) == bell_number(n,'dobinski') for n in range(200))
+        sage: all(bell_number(n) == bell_number(n,'dobinski') for n in range(100))
         True
-        sage: all(bell_number(n) == bell_number(n,'gap') for n in range(200))
+        sage: all(bell_number(n) == bell_number(n,'gap') for n in range(100))
         True
         sage: all(bell_number(n) == bell_number(n,'mpmath', prec=500) for n in range(200, 220))
         True
@@ -854,16 +854,23 @@ def lucas_number2(n, P, Q):
     return libgap.Lucas(P, Q, n)[1].sage()
 
 
-def stirling_number1(n, k) -> Integer:
+def stirling_number1(n, k, algorithm="gap") -> Integer:
     r"""
     Return the `n`-th Stirling number `S_1(n,k)` of the first kind.
 
     This is the number of permutations of `n` points with `k` cycles.
 
-    This wraps GAP's ``Stirling1``.
-
     See :wikipedia:`Stirling_numbers_of_the_first_kind`.
 
+    INPUT:
+
+    - ``n`` -- nonnegative machine-size integer
+    - ``k`` -- nonnegative machine-size integer
+    - ``algorithm``:
+
+      * ``"gap"`` (default) -- use GAP's ``Stirling1`` function
+      * ``"flint"`` -- use flint's ``arith_stirling_number_1u`` function
+
     EXAMPLES::
 
         sage: stirling_number1(3,2)
@@ -876,31 +883,48 @@ def stirling_number1(n, k) -> Integer:
         269325
 
     Indeed, `S_1(n,k) = S_1(n-1,k-1) + (n-1)S_1(n-1,k)`.
+
+    TESTS::
+
+        sage: stirling_number1(10,5, algorithm='flint')
+        269325
+
+        sage: s_sage = stirling_number1(50,3, algorithm="mutta")
+        Traceback (most recent call last):
+        ...
+        ValueError: unknown algorithm: mutta
     """
     n = ZZ(n)
     k = ZZ(k)
-    from sage.libs.gap.libgap import libgap
-    return libgap.Stirling1(n, k).sage()
+    if algorithm == 'gap':
+        from sage.libs.gap.libgap import libgap
+        return libgap.Stirling1(n, k).sage()
+    if algorithm == 'flint':
+        import sage.libs.flint.arith
+        return sage.libs.flint.arith.stirling_number_1(n, k)
+    raise ValueError("unknown algorithm: %s" % algorithm)
 
 
 def stirling_number2(n, k, algorithm=None) -> Integer:
     r"""
-    Return the `n`-th Stirling number `S_2(n,k)` of the second
-    kind.
+    Return the `n`-th Stirling number `S_2(n,k)` of the second kind.
 
     This is the number of ways to partition a set of `n` elements into `k`
     pairwise disjoint nonempty subsets. The `n`-th Bell number is the
     sum of the `S_2(n,k)`'s, `k=0,...,n`.
 
+    See :wikipedia:`Stirling_numbers_of_the_second_kind`.
+
     INPUT:
 
-       *  ``n`` - nonnegative machine-size integer
-       *  ``k`` - nonnegative machine-size integer
-       * ``algorithm``:
+    -  ``n`` -- nonnegative machine-size integer
+    -  ``k`` -- nonnegative machine-size integer
+    - ``algorithm``:
 
-         * None (default) - use native implementation
-         * ``"maxima"`` - use Maxima's stirling2 function
-         * ``"gap"`` - use GAP's Stirling2 function
+      * ``None`` (default) -- use native implementation
+      * ``"flint"`` -- use flint's ``arith_stirling_number_2`` function
+      * ``"gap"`` -- use GAP's ``Stirling2`` function
+      * ``"maxima"`` -- use Maxima's ``stirling2`` function
 
     EXAMPLES:
 
@@ -922,10 +946,10 @@ def stirling_number2(n, k, algorithm=None) -> Integer:
 
     Stirling numbers satisfy `S_2(n,k) = S_2(n-1,k-1) + kS_2(n-1,k)`::
 
-         sage: 5*stirling_number2(9,5) + stirling_number2(9,4)
-         42525
-         sage: stirling_number2(10,5)
-         42525
+        sage: 5*stirling_number2(9,5) + stirling_number2(9,4)
+        42525
+        sage: stirling_number2(10,5)
+        42525
 
     TESTS::
 
@@ -961,58 +985,57 @@ def stirling_number2(n, k, algorithm=None) -> Integer:
         1900842429486
         sage: type(n)
         <class 'sage.rings.integer.Integer'>
-        sage: n = stirling_number2(20,11,algorithm='maxima')
+        sage: n = stirling_number2(20,11,algorithm='flint')
         sage: n
         1900842429486
         sage: type(n)
         <class 'sage.rings.integer.Integer'>
 
-     Sage's implementation splitting the computation of the Stirling
-     numbers of the second kind in two cases according to `n`, let us
-     check the result it gives agree with both maxima and gap.
-
-     For `n<200`::
+    Sage's implementation splitting the computation of the Stirling
+    numbers of the second kind in two cases according to `n`, let us
+    check the result it gives agree with both flint and gap.
 
-         sage: for n in Subsets(range(100,200), 5).random_element():
-         ....:     for k in Subsets(range(n), 5).random_element():
-         ....:         s_sage = stirling_number2(n,k)
-         ....:         s_maxima = stirling_number2(n,k, algorithm = "maxima")
-         ....:         s_gap = stirling_number2(n,k, algorithm = "gap")
-         ....:         if not (s_sage == s_maxima and s_sage == s_gap):
-         ....:             print("Error with n<200")
+    For `n<200`::
 
-     For `n\geq 200`::
+        sage: for n in Subsets(range(100,200), 5).random_element():
+        ....:     for k in Subsets(range(n), 5).random_element():
+        ....:         s_sage = stirling_number2(n,k)
+        ....:         s_flint = stirling_number2(n,k, algorithm = "flint")
+        ....:         s_gap = stirling_number2(n,k, algorithm = "gap")
+        ....:         if not (s_sage == s_flint and s_sage == s_gap):
+        ....:             print("Error with n<200")
 
-         sage: for n in Subsets(range(200,300), 5).random_element():
-         ....:     for k in Subsets(range(n), 5).random_element():
-         ....:         s_sage = stirling_number2(n,k)
-         ....:         s_maxima = stirling_number2(n,k, algorithm = "maxima")
-         ....:         s_gap = stirling_number2(n,k, algorithm = "gap")
-         ....:         if not (s_sage == s_maxima and s_sage == s_gap):
-         ....:             print("Error with n<200")
+    For `n\geq 200`::
 
+        sage: for n in Subsets(range(200,300), 5).random_element():
+        ....:     for k in Subsets(range(n), 5).random_element():
+        ....:         s_sage = stirling_number2(n,k)
+        ....:         s_flint = stirling_number2(n,k, algorithm = "flint")
+        ....:         s_gap = stirling_number2(n,k, algorithm = "gap")
+        ....:         if not (s_sage == s_flint and s_sage == s_gap):
+        ....:             print("Error with n<200")
 
-     TESTS:
+        sage: stirling_number2(20,3, algorithm="maxima")
+        580606446
 
-     Checking an exception is raised whenever a wrong value is given
-     for ``algorithm``::
-
-         sage: s_sage = stirling_number2(50,3, algorithm = "CloudReading")
-         Traceback (most recent call last):
-         ...
-         ValueError: unknown algorithm: CloudReading
+        sage: s_sage = stirling_number2(5,3, algorithm="namba")
+        Traceback (most recent call last):
+        ...
+        ValueError: unknown algorithm: namba
     """
     n = ZZ(n)
     k = ZZ(k)
     if algorithm is None:
         return _stirling_number2(n, k)
-    elif algorithm == 'gap':
+    if algorithm == 'gap':
         from sage.libs.gap.libgap import libgap
         return libgap.Stirling2(n, k).sage()
-    elif algorithm == 'maxima':
+    if algorithm == 'flint':
+        import sage.libs.flint.arith
+        return sage.libs.flint.arith.stirling_number_2(n, k)
+    if algorithm == 'maxima':
         return ZZ(maxima.stirling2(n, k))  # type:ignore
-    else:
-        raise ValueError("unknown algorithm: %s" % algorithm)
+    raise ValueError("unknown algorithm: %s" % algorithm)
 
 
 def polygonal_number(s, n):
@@ -1404,8 +1427,6 @@ def __bool__(self) -> bool:
         """
         return bool(self._list)
 
-
-
     def __len__(self) -> Integer:
         """
         EXAMPLES::
@@ -2442,7 +2463,6 @@ def __iter__(self):
 
 ##############################################################################
 from sage.sets.image_set import ImageSubobject
-from sage.categories.map import is_Map
 
 
 class MapCombinatorialClass(ImageSubobject, CombinatorialClass):

src/sage/libs/flint/arith.pxd

@@ -1,13 +1,15 @@
 # distutils: libraries = flint
 # distutils: depends = flint/arith.h
 
-from sage.libs.flint.types cimport fmpz_t, fmpq_t, ulong
+from sage.libs.flint.types cimport fmpz_t, fmpq_t, ulong, slong
 
 # flint/arith.h
 cdef extern from "flint_wrap.h":
     void arith_bell_number(fmpz_t b, ulong n)
     void arith_bernoulli_number(fmpq_t x, ulong n)
-    void arith_euler_number ( fmpz_t res , ulong n )
+    void arith_euler_number(fmpz_t res, ulong n)
+    void arith_stirling_number_1u(fmpz_t res, slong n, slong k)
+    void arith_stirling_number_2(fmpz_t res, slong n, slong k)
     void arith_number_of_partitions(fmpz_t x, ulong n)
     void arith_dedekind_sum(fmpq_t, fmpz_t, fmpz_t)
     void arith_harmonic_number(fmpq_t, unsigned long n)

src/sage/libs/flint/arith.pyx

@@ -116,6 +116,56 @@ def euler_number(unsigned long n):
     return ans
 
 
+def stirling_number_1(long n, long k):
+    """
+    Return the unsigned Stirling number of the first kind.
+
+    EXAMPLES::
+
+        sage: from sage.libs.flint.arith import stirling_number_1
+        sage: [stirling_number_1(8,i) for i in range(9)]
+        [0, 5040, 13068, 13132, 6769, 1960, 322, 28, 1]
+    """
+    cdef fmpz_t ans_fmpz
+    cdef Integer ans = Integer(0)
+
+    fmpz_init(ans_fmpz)
+
+    if n > 1000:
+        sig_on()
+    arith_stirling_number_1u(ans_fmpz, n, k)
+    fmpz_get_mpz(ans.value, ans_fmpz)
+    fmpz_clear(ans_fmpz)
+    if n > 1000:
+        sig_off()
+    return ans
+
+
+def stirling_number_2(long n, long k):
+    """
+    Return the Stirling number of the second kind.
+
+    EXAMPLES::
+
+        sage: from sage.libs.flint.arith import stirling_number_2
+        sage: [stirling_number_2(8,i) for i in range(9)]
+        [0, 1, 127, 966, 1701, 1050, 266, 28, 1]
+    """
+    cdef fmpz_t ans_fmpz
+    cdef Integer ans = Integer(0)
+
+    fmpz_init(ans_fmpz)
+
+    if n > 1000:
+        sig_on()
+    arith_stirling_number_2(ans_fmpz, n, k)
+    fmpz_get_mpz(ans.value, ans_fmpz)
+    fmpz_clear(ans_fmpz)
+    if n > 1000:
+        sig_off()
+    return ans
+
+
 def number_of_partitions(unsigned long n):
     """
     Return the number of partitions of the integer `n`.