diff --git a/src/sage/combinat/combinat.py b/src/sage/combinat/combinat.py index 42fed39d896..48475d117dc 100644 --- a/src/sage/combinat/combinat.py +++ b/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) - sage: n = stirling_number2(20,11,algorithm='maxima') + sage: n = stirling_number2(20,11,algorithm='flint') sage: n 1900842429486 sage: type(n) - 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): diff --git a/src/sage/libs/flint/arith.pxd b/src/sage/libs/flint/arith.pxd index d4c63fcaf13..c8e1fb35566 100644 --- a/src/sage/libs/flint/arith.pxd +++ b/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) diff --git a/src/sage/libs/flint/arith.pyx b/src/sage/libs/flint/arith.pyx index 4a378ce9862..f24446b4ae4 100644 --- a/src/sage/libs/flint/arith.pyx +++ b/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`.