Coverage for /usr/lib/python3/dist-packages/sympy/functions/combinatorial/numbers.py: 18%
1011 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
2This module implements some special functions that commonly appear in
3combinatorial contexts (e.g. in power series); in particular,
4sequences of rational numbers such as Bernoulli and Fibonacci numbers.
6Factorials, binomial coefficients and related functions are located in
7the separate 'factorials' module.
8"""
9from math import prod
10from collections import defaultdict
11from typing import Tuple as tTuple
13from sympy.core import S, Symbol, Add, Dummy
14from sympy.core.cache import cacheit
15from sympy.core.expr import Expr
16from sympy.core.function import ArgumentIndexError, Function, expand_mul
17from sympy.core.logic import fuzzy_not
18from sympy.core.mul import Mul
19from sympy.core.numbers import E, I, pi, oo, Rational, Integer
20from sympy.core.relational import Eq, is_le, is_gt
21from sympy.external.gmpy import SYMPY_INTS
22from sympy.functions.combinatorial.factorials import (binomial,
23 factorial, subfactorial)
24from sympy.functions.elementary.exponential import log
25from sympy.functions.elementary.piecewise import Piecewise
26from sympy.ntheory.primetest import isprime, is_square
27from sympy.polys.appellseqs import bernoulli_poly, euler_poly, genocchi_poly
28from sympy.utilities.enumerative import MultisetPartitionTraverser
29from sympy.utilities.exceptions import sympy_deprecation_warning
30from sympy.utilities.iterables import multiset, multiset_derangements, iterable
31from sympy.utilities.memoization import recurrence_memo
32from sympy.utilities.misc import as_int
34from mpmath import mp, workprec
35from mpmath.libmp import ifib as _ifib
38def _product(a, b):
39 return prod(range(a, b + 1))
42# Dummy symbol used for computing polynomial sequences
43_sym = Symbol('x')
46#----------------------------------------------------------------------------#
47# #
48# Carmichael numbers #
49# #
50#----------------------------------------------------------------------------#
52def _divides(p, n):
53 return n % p == 0
55class carmichael(Function):
56 r"""
57 Carmichael Numbers:
59 Certain cryptographic algorithms make use of big prime numbers.
60 However, checking whether a big number is prime is not so easy.
61 Randomized prime number checking tests exist that offer a high degree of
62 confidence of accurate determination at low cost, such as the Fermat test.
64 Let 'a' be a random number between $2$ and $n - 1$, where $n$ is the
65 number whose primality we are testing. Then, $n$ is probably prime if it
66 satisfies the modular arithmetic congruence relation:
68 .. math :: a^{n-1} = 1 \pmod{n}
70 (where mod refers to the modulo operation)
72 If a number passes the Fermat test several times, then it is prime with a
73 high probability.
75 Unfortunately, certain composite numbers (non-primes) still pass the Fermat
76 test with every number smaller than themselves.
77 These numbers are called Carmichael numbers.
79 A Carmichael number will pass a Fermat primality test to every base $b$
80 relatively prime to the number, even though it is not actually prime.
81 This makes tests based on Fermat's Little Theorem less effective than
82 strong probable prime tests such as the Baillie-PSW primality test and
83 the Miller-Rabin primality test.
85 Examples
86 ========
88 >>> from sympy import carmichael
89 >>> carmichael.find_first_n_carmichaels(5)
90 [561, 1105, 1729, 2465, 2821]
91 >>> carmichael.find_carmichael_numbers_in_range(0, 562)
92 [561]
93 >>> carmichael.find_carmichael_numbers_in_range(0,1000)
94 [561]
95 >>> carmichael.find_carmichael_numbers_in_range(0,2000)
96 [561, 1105, 1729]
98 References
99 ==========
101 .. [1] https://en.wikipedia.org/wiki/Carmichael_number
102 .. [2] https://en.wikipedia.org/wiki/Fermat_primality_test
103 .. [3] https://www.jstor.org/stable/23248683?seq=1#metadata_info_tab_contents
104 """
106 @staticmethod
107 def is_perfect_square(n):
108 sympy_deprecation_warning(
109 """
110is_perfect_square is just a wrapper around sympy.ntheory.primetest.is_square
111so use that directly instead.
112 """,
113 deprecated_since_version="1.11",
114 active_deprecations_target='deprecated-carmichael-static-methods',
115 )
116 return is_square(n)
118 @staticmethod
119 def divides(p, n):
120 sympy_deprecation_warning(
121 """
122 divides can be replaced by directly testing n % p == 0.
123 """,
124 deprecated_since_version="1.11",
125 active_deprecations_target='deprecated-carmichael-static-methods',
126 )
127 return n % p == 0
129 @staticmethod
130 def is_prime(n):
131 sympy_deprecation_warning(
132 """
133is_prime is just a wrapper around sympy.ntheory.primetest.isprime so use that
134directly instead.
135 """,
136 deprecated_since_version="1.11",
137 active_deprecations_target='deprecated-carmichael-static-methods',
138 )
139 return isprime(n)
141 @staticmethod
142 def is_carmichael(n):
143 if n >= 0:
144 if (n == 1) or isprime(n) or (n % 2 == 0):
145 return False
147 divisors = [1, n]
149 # get divisors
150 divisors.extend([i for i in range(3, n // 2 + 1, 2) if n % i == 0])
152 for i in divisors:
153 if is_square(i) and i != 1:
154 return False
155 if isprime(i):
156 if not _divides(i - 1, n - 1):
157 return False
159 return True
161 else:
162 raise ValueError('The provided number must be greater than or equal to 0')
164 @staticmethod
165 def find_carmichael_numbers_in_range(x, y):
166 if 0 <= x <= y:
167 if x % 2 == 0:
168 return [i for i in range(x + 1, y, 2) if carmichael.is_carmichael(i)]
169 else:
170 return [i for i in range(x, y, 2) if carmichael.is_carmichael(i)]
172 else:
173 raise ValueError('The provided range is not valid. x and y must be non-negative integers and x <= y')
175 @staticmethod
176 def find_first_n_carmichaels(n):
177 i = 1
178 carmichaels = []
180 while len(carmichaels) < n:
181 if carmichael.is_carmichael(i):
182 carmichaels.append(i)
183 i += 2
185 return carmichaels
188#----------------------------------------------------------------------------#
189# #
190# Fibonacci numbers #
191# #
192#----------------------------------------------------------------------------#
195class fibonacci(Function):
196 r"""
197 Fibonacci numbers / Fibonacci polynomials
199 The Fibonacci numbers are the integer sequence defined by the
200 initial terms `F_0 = 0`, `F_1 = 1` and the two-term recurrence
201 relation `F_n = F_{n-1} + F_{n-2}`. This definition
202 extended to arbitrary real and complex arguments using
203 the formula
205 .. math :: F_z = \frac{\phi^z - \cos(\pi z) \phi^{-z}}{\sqrt 5}
207 The Fibonacci polynomials are defined by `F_1(x) = 1`,
208 `F_2(x) = x`, and `F_n(x) = x*F_{n-1}(x) + F_{n-2}(x)` for `n > 2`.
209 For all positive integers `n`, `F_n(1) = F_n`.
211 * ``fibonacci(n)`` gives the `n^{th}` Fibonacci number, `F_n`
212 * ``fibonacci(n, x)`` gives the `n^{th}` Fibonacci polynomial in `x`, `F_n(x)`
214 Examples
215 ========
217 >>> from sympy import fibonacci, Symbol
219 >>> [fibonacci(x) for x in range(11)]
220 [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
221 >>> fibonacci(5, Symbol('t'))
222 t**4 + 3*t**2 + 1
224 See Also
225 ========
227 bell, bernoulli, catalan, euler, harmonic, lucas, genocchi, partition, tribonacci
229 References
230 ==========
232 .. [1] https://en.wikipedia.org/wiki/Fibonacci_number
233 .. [2] https://mathworld.wolfram.com/FibonacciNumber.html
235 """
237 @staticmethod
238 def _fib(n):
239 return _ifib(n)
241 @staticmethod
242 @recurrence_memo([None, S.One, _sym])
243 def _fibpoly(n, prev):
244 return (prev[-2] + _sym*prev[-1]).expand()
246 @classmethod
247 def eval(cls, n, sym=None):
248 if n is S.Infinity:
249 return S.Infinity
251 if n.is_Integer:
252 if sym is None:
253 n = int(n)
254 if n < 0:
255 return S.NegativeOne**(n + 1) * fibonacci(-n)
256 else:
257 return Integer(cls._fib(n))
258 else:
259 if n < 1:
260 raise ValueError("Fibonacci polynomials are defined "
261 "only for positive integer indices.")
262 return cls._fibpoly(n).subs(_sym, sym)
264 def _eval_rewrite_as_sqrt(self, n, **kwargs):
265 from sympy.functions.elementary.miscellaneous import sqrt
266 return 2**(-n)*sqrt(5)*((1 + sqrt(5))**n - (-sqrt(5) + 1)**n) / 5
268 def _eval_rewrite_as_GoldenRatio(self,n, **kwargs):
269 return (S.GoldenRatio**n - 1/(-S.GoldenRatio)**n)/(2*S.GoldenRatio-1)
272#----------------------------------------------------------------------------#
273# #
274# Lucas numbers #
275# #
276#----------------------------------------------------------------------------#
279class lucas(Function):
280 """
281 Lucas numbers
283 Lucas numbers satisfy a recurrence relation similar to that of
284 the Fibonacci sequence, in which each term is the sum of the
285 preceding two. They are generated by choosing the initial
286 values `L_0 = 2` and `L_1 = 1`.
288 * ``lucas(n)`` gives the `n^{th}` Lucas number
290 Examples
291 ========
293 >>> from sympy import lucas
295 >>> [lucas(x) for x in range(11)]
296 [2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123]
298 See Also
299 ========
301 bell, bernoulli, catalan, euler, fibonacci, harmonic, genocchi, partition, tribonacci
303 References
304 ==========
306 .. [1] https://en.wikipedia.org/wiki/Lucas_number
307 .. [2] https://mathworld.wolfram.com/LucasNumber.html
309 """
311 @classmethod
312 def eval(cls, n):
313 if n is S.Infinity:
314 return S.Infinity
316 if n.is_Integer:
317 return fibonacci(n + 1) + fibonacci(n - 1)
319 def _eval_rewrite_as_sqrt(self, n, **kwargs):
320 from sympy.functions.elementary.miscellaneous import sqrt
321 return 2**(-n)*((1 + sqrt(5))**n + (-sqrt(5) + 1)**n)
324#----------------------------------------------------------------------------#
325# #
326# Tribonacci numbers #
327# #
328#----------------------------------------------------------------------------#
331class tribonacci(Function):
332 r"""
333 Tribonacci numbers / Tribonacci polynomials
335 The Tribonacci numbers are the integer sequence defined by the
336 initial terms `T_0 = 0`, `T_1 = 1`, `T_2 = 1` and the three-term
337 recurrence relation `T_n = T_{n-1} + T_{n-2} + T_{n-3}`.
339 The Tribonacci polynomials are defined by `T_0(x) = 0`, `T_1(x) = 1`,
340 `T_2(x) = x^2`, and `T_n(x) = x^2 T_{n-1}(x) + x T_{n-2}(x) + T_{n-3}(x)`
341 for `n > 2`. For all positive integers `n`, `T_n(1) = T_n`.
343 * ``tribonacci(n)`` gives the `n^{th}` Tribonacci number, `T_n`
344 * ``tribonacci(n, x)`` gives the `n^{th}` Tribonacci polynomial in `x`, `T_n(x)`
346 Examples
347 ========
349 >>> from sympy import tribonacci, Symbol
351 >>> [tribonacci(x) for x in range(11)]
352 [0, 1, 1, 2, 4, 7, 13, 24, 44, 81, 149]
353 >>> tribonacci(5, Symbol('t'))
354 t**8 + 3*t**5 + 3*t**2
356 See Also
357 ========
359 bell, bernoulli, catalan, euler, fibonacci, harmonic, lucas, genocchi, partition
361 References
362 ==========
364 .. [1] https://en.wikipedia.org/wiki/Generalizations_of_Fibonacci_numbers#Tribonacci_numbers
365 .. [2] https://mathworld.wolfram.com/TribonacciNumber.html
366 .. [3] https://oeis.org/A000073
368 """
370 @staticmethod
371 @recurrence_memo([S.Zero, S.One, S.One])
372 def _trib(n, prev):
373 return (prev[-3] + prev[-2] + prev[-1])
375 @staticmethod
376 @recurrence_memo([S.Zero, S.One, _sym**2])
377 def _tribpoly(n, prev):
378 return (prev[-3] + _sym*prev[-2] + _sym**2*prev[-1]).expand()
380 @classmethod
381 def eval(cls, n, sym=None):
382 if n is S.Infinity:
383 return S.Infinity
385 if n.is_Integer:
386 n = int(n)
387 if n < 0:
388 raise ValueError("Tribonacci polynomials are defined "
389 "only for non-negative integer indices.")
390 if sym is None:
391 return Integer(cls._trib(n))
392 else:
393 return cls._tribpoly(n).subs(_sym, sym)
395 def _eval_rewrite_as_sqrt(self, n, **kwargs):
396 from sympy.functions.elementary.miscellaneous import cbrt, sqrt
397 w = (-1 + S.ImaginaryUnit * sqrt(3)) / 2
398 a = (1 + cbrt(19 + 3*sqrt(33)) + cbrt(19 - 3*sqrt(33))) / 3
399 b = (1 + w*cbrt(19 + 3*sqrt(33)) + w**2*cbrt(19 - 3*sqrt(33))) / 3
400 c = (1 + w**2*cbrt(19 + 3*sqrt(33)) + w*cbrt(19 - 3*sqrt(33))) / 3
401 Tn = (a**(n + 1)/((a - b)*(a - c))
402 + b**(n + 1)/((b - a)*(b - c))
403 + c**(n + 1)/((c - a)*(c - b)))
404 return Tn
406 def _eval_rewrite_as_TribonacciConstant(self, n, **kwargs):
407 from sympy.functions.elementary.integers import floor
408 from sympy.functions.elementary.miscellaneous import cbrt, sqrt
409 b = cbrt(586 + 102*sqrt(33))
410 Tn = 3 * b * S.TribonacciConstant**n / (b**2 - 2*b + 4)
411 return floor(Tn + S.Half)
414#----------------------------------------------------------------------------#
415# #
416# Bernoulli numbers #
417# #
418#----------------------------------------------------------------------------#
421class bernoulli(Function):
422 r"""
423 Bernoulli numbers / Bernoulli polynomials / Bernoulli function
425 The Bernoulli numbers are a sequence of rational numbers
426 defined by `B_0 = 1` and the recursive relation (`n > 0`):
428 .. math :: n+1 = \sum_{k=0}^n \binom{n+1}{k} B_k
430 They are also commonly defined by their exponential generating
431 function, which is `\frac{x}{1 - e^{-x}}`. For odd indices > 1,
432 the Bernoulli numbers are zero.
434 The Bernoulli polynomials satisfy the analogous formula:
436 .. math :: B_n(x) = \sum_{k=0}^n (-1)^k \binom{n}{k} B_k x^{n-k}
438 Bernoulli numbers and Bernoulli polynomials are related as
439 `B_n(1) = B_n`.
441 The generalized Bernoulli function `\operatorname{B}(s, a)`
442 is defined for any complex `s` and `a`, except where `a` is a
443 nonpositive integer and `s` is not a nonnegative integer. It is
444 an entire function of `s` for fixed `a`, related to the Hurwitz
445 zeta function by
447 .. math:: \operatorname{B}(s, a) = \begin{cases}
448 -s \zeta(1-s, a) & s \ne 0 \\ 1 & s = 0 \end{cases}
450 When `s` is a nonnegative integer this function reduces to the
451 Bernoulli polynomials: `\operatorname{B}(n, x) = B_n(x)`. When
452 `a` is omitted it is assumed to be 1, yielding the (ordinary)
453 Bernoulli function which interpolates the Bernoulli numbers and is
454 related to the Riemann zeta function.
456 We compute Bernoulli numbers using Ramanujan's formula:
458 .. math :: B_n = \frac{A(n) - S(n)}{\binom{n+3}{n}}
460 where:
462 .. math :: A(n) = \begin{cases} \frac{n+3}{3} &
463 n \equiv 0\ \text{or}\ 2 \pmod{6} \\
464 -\frac{n+3}{6} & n \equiv 4 \pmod{6} \end{cases}
466 and:
468 .. math :: S(n) = \sum_{k=1}^{[n/6]} \binom{n+3}{n-6k} B_{n-6k}
470 This formula is similar to the sum given in the definition, but
471 cuts `\frac{2}{3}` of the terms. For Bernoulli polynomials, we use
472 Appell sequences.
474 For `n` a nonnegative integer and `s`, `a`, `x` arbitrary complex numbers,
476 * ``bernoulli(n)`` gives the nth Bernoulli number, `B_n`
477 * ``bernoulli(s)`` gives the Bernoulli function `\operatorname{B}(s)`
478 * ``bernoulli(n, x)`` gives the nth Bernoulli polynomial in `x`, `B_n(x)`
479 * ``bernoulli(s, a)`` gives the generalized Bernoulli function
480 `\operatorname{B}(s, a)`
482 .. versionchanged:: 1.12
483 ``bernoulli(1)`` gives `+\frac{1}{2}` instead of `-\frac{1}{2}`.
484 This choice of value confers several theoretical advantages [5]_,
485 including the extension to complex parameters described above
486 which this function now implements. The previous behavior, defined
487 only for nonnegative integers `n`, can be obtained with
488 ``(-1)**n*bernoulli(n)``.
490 Examples
491 ========
493 >>> from sympy import bernoulli
494 >>> from sympy.abc import x
495 >>> [bernoulli(n) for n in range(11)]
496 [1, 1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0, 5/66]
497 >>> bernoulli(1000001)
498 0
499 >>> bernoulli(3, x)
500 x**3 - 3*x**2/2 + x/2
502 See Also
503 ========
505 andre, bell, catalan, euler, fibonacci, harmonic, lucas, genocchi,
506 partition, tribonacci, sympy.polys.appellseqs.bernoulli_poly
508 References
509 ==========
511 .. [1] https://en.wikipedia.org/wiki/Bernoulli_number
512 .. [2] https://en.wikipedia.org/wiki/Bernoulli_polynomial
513 .. [3] https://mathworld.wolfram.com/BernoulliNumber.html
514 .. [4] https://mathworld.wolfram.com/BernoulliPolynomial.html
515 .. [5] Peter Luschny, "The Bernoulli Manifesto",
516 https://luschny.de/math/zeta/The-Bernoulli-Manifesto.html
517 .. [6] Peter Luschny, "An introduction to the Bernoulli function",
518 https://arxiv.org/abs/2009.06743
520 """
522 args: tTuple[Integer]
524 # Calculates B_n for positive even n
525 @staticmethod
526 def _calc_bernoulli(n):
527 s = 0
528 a = int(binomial(n + 3, n - 6))
529 for j in range(1, n//6 + 1):
530 s += a * bernoulli(n - 6*j)
531 # Avoid computing each binomial coefficient from scratch
532 a *= _product(n - 6 - 6*j + 1, n - 6*j)
533 a //= _product(6*j + 4, 6*j + 9)
534 if n % 6 == 4:
535 s = -Rational(n + 3, 6) - s
536 else:
537 s = Rational(n + 3, 3) - s
538 return s / binomial(n + 3, n)
540 # We implement a specialized memoization scheme to handle each
541 # case modulo 6 separately
542 _cache = {0: S.One, 2: Rational(1, 6), 4: Rational(-1, 30)}
543 _highest = {0: 0, 2: 2, 4: 4}
545 @classmethod
546 def eval(cls, n, x=None):
547 if x is S.One:
548 return cls(n)
549 elif n.is_zero:
550 return S.One
551 elif n.is_integer is False or n.is_nonnegative is False:
552 if x is not None and x.is_Integer and x.is_nonpositive:
553 return S.NaN
554 return
555 # Bernoulli numbers
556 elif x is None:
557 if n is S.One:
558 return S.Half
559 elif n.is_odd and (n-1).is_positive:
560 return S.Zero
561 elif n.is_Number:
562 n = int(n)
563 # Use mpmath for enormous Bernoulli numbers
564 if n > 500:
565 p, q = mp.bernfrac(n)
566 return Rational(int(p), int(q))
567 case = n % 6
568 highest_cached = cls._highest[case]
569 if n <= highest_cached:
570 return cls._cache[n]
571 # To avoid excessive recursion when, say, bernoulli(1000) is
572 # requested, calculate and cache the entire sequence ... B_988,
573 # B_994, B_1000 in increasing order
574 for i in range(highest_cached + 6, n + 6, 6):
575 b = cls._calc_bernoulli(i)
576 cls._cache[i] = b
577 cls._highest[case] = i
578 return b
579 # Bernoulli polynomials
580 elif n.is_Number:
581 return bernoulli_poly(n, x)
583 def _eval_rewrite_as_zeta(self, n, x=1, **kwargs):
584 from sympy.functions.special.zeta_functions import zeta
585 return Piecewise((1, Eq(n, 0)), (-n * zeta(1-n, x), True))
587 def _eval_evalf(self, prec):
588 if not all(x.is_number for x in self.args):
589 return
590 n = self.args[0]._to_mpmath(prec)
591 x = (self.args[1] if len(self.args) > 1 else S.One)._to_mpmath(prec)
592 with workprec(prec):
593 if n == 0:
594 res = mp.mpf(1)
595 elif n == 1:
596 res = x - mp.mpf(0.5)
597 elif mp.isint(n) and n >= 0:
598 res = mp.bernoulli(n) if x == 1 else mp.bernpoly(n, x)
599 else:
600 res = -n * mp.zeta(1-n, x)
601 return Expr._from_mpmath(res, prec)
604#----------------------------------------------------------------------------#
605# #
606# Bell numbers #
607# #
608#----------------------------------------------------------------------------#
611class bell(Function):
612 r"""
613 Bell numbers / Bell polynomials
615 The Bell numbers satisfy `B_0 = 1` and
617 .. math:: B_n = \sum_{k=0}^{n-1} \binom{n-1}{k} B_k.
619 They are also given by:
621 .. math:: B_n = \frac{1}{e} \sum_{k=0}^{\infty} \frac{k^n}{k!}.
623 The Bell polynomials are given by `B_0(x) = 1` and
625 .. math:: B_n(x) = x \sum_{k=1}^{n-1} \binom{n-1}{k-1} B_{k-1}(x).
627 The second kind of Bell polynomials (are sometimes called "partial" Bell
628 polynomials or incomplete Bell polynomials) are defined as
630 .. math:: B_{n,k}(x_1, x_2,\dotsc x_{n-k+1}) =
631 \sum_{j_1+j_2+j_2+\dotsb=k \atop j_1+2j_2+3j_2+\dotsb=n}
632 \frac{n!}{j_1!j_2!\dotsb j_{n-k+1}!}
633 \left(\frac{x_1}{1!} \right)^{j_1}
634 \left(\frac{x_2}{2!} \right)^{j_2} \dotsb
635 \left(\frac{x_{n-k+1}}{(n-k+1)!} \right) ^{j_{n-k+1}}.
637 * ``bell(n)`` gives the `n^{th}` Bell number, `B_n`.
638 * ``bell(n, x)`` gives the `n^{th}` Bell polynomial, `B_n(x)`.
639 * ``bell(n, k, (x1, x2, ...))`` gives Bell polynomials of the second kind,
640 `B_{n,k}(x_1, x_2, \dotsc, x_{n-k+1})`.
642 Notes
643 =====
645 Not to be confused with Bernoulli numbers and Bernoulli polynomials,
646 which use the same notation.
648 Examples
649 ========
651 >>> from sympy import bell, Symbol, symbols
653 >>> [bell(n) for n in range(11)]
654 [1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975]
655 >>> bell(30)
656 846749014511809332450147
657 >>> bell(4, Symbol('t'))
658 t**4 + 6*t**3 + 7*t**2 + t
659 >>> bell(6, 2, symbols('x:6')[1:])
660 6*x1*x5 + 15*x2*x4 + 10*x3**2
662 See Also
663 ========
665 bernoulli, catalan, euler, fibonacci, harmonic, lucas, genocchi, partition, tribonacci
667 References
668 ==========
670 .. [1] https://en.wikipedia.org/wiki/Bell_number
671 .. [2] https://mathworld.wolfram.com/BellNumber.html
672 .. [3] https://mathworld.wolfram.com/BellPolynomial.html
674 """
676 @staticmethod
677 @recurrence_memo([1, 1])
678 def _bell(n, prev):
679 s = 1
680 a = 1
681 for k in range(1, n):
682 a = a * (n - k) // k
683 s += a * prev[k]
684 return s
686 @staticmethod
687 @recurrence_memo([S.One, _sym])
688 def _bell_poly(n, prev):
689 s = 1
690 a = 1
691 for k in range(2, n + 1):
692 a = a * (n - k + 1) // (k - 1)
693 s += a * prev[k - 1]
694 return expand_mul(_sym * s)
696 @staticmethod
697 def _bell_incomplete_poly(n, k, symbols):
698 r"""
699 The second kind of Bell polynomials (incomplete Bell polynomials).
701 Calculated by recurrence formula:
703 .. math:: B_{n,k}(x_1, x_2, \dotsc, x_{n-k+1}) =
704 \sum_{m=1}^{n-k+1}
705 \x_m \binom{n-1}{m-1} B_{n-m,k-1}(x_1, x_2, \dotsc, x_{n-m-k})
707 where
708 `B_{0,0} = 1;`
709 `B_{n,0} = 0; for n \ge 1`
710 `B_{0,k} = 0; for k \ge 1`
712 """
713 if (n == 0) and (k == 0):
714 return S.One
715 elif (n == 0) or (k == 0):
716 return S.Zero
717 s = S.Zero
718 a = S.One
719 for m in range(1, n - k + 2):
720 s += a * bell._bell_incomplete_poly(
721 n - m, k - 1, symbols) * symbols[m - 1]
722 a = a * (n - m) / m
723 return expand_mul(s)
725 @classmethod
726 def eval(cls, n, k_sym=None, symbols=None):
727 if n is S.Infinity:
728 if k_sym is None:
729 return S.Infinity
730 else:
731 raise ValueError("Bell polynomial is not defined")
733 if n.is_negative or n.is_integer is False:
734 raise ValueError("a non-negative integer expected")
736 if n.is_Integer and n.is_nonnegative:
737 if k_sym is None:
738 return Integer(cls._bell(int(n)))
739 elif symbols is None:
740 return cls._bell_poly(int(n)).subs(_sym, k_sym)
741 else:
742 r = cls._bell_incomplete_poly(int(n), int(k_sym), symbols)
743 return r
745 def _eval_rewrite_as_Sum(self, n, k_sym=None, symbols=None, **kwargs):
746 from sympy.concrete.summations import Sum
747 if (k_sym is not None) or (symbols is not None):
748 return self
750 # Dobinski's formula
751 if not n.is_nonnegative:
752 return self
753 k = Dummy('k', integer=True, nonnegative=True)
754 return 1 / E * Sum(k**n / factorial(k), (k, 0, S.Infinity))
757#----------------------------------------------------------------------------#
758# #
759# Harmonic numbers #
760# #
761#----------------------------------------------------------------------------#
764class harmonic(Function):
765 r"""
766 Harmonic numbers
768 The nth harmonic number is given by `\operatorname{H}_{n} =
769 1 + \frac{1}{2} + \frac{1}{3} + \ldots + \frac{1}{n}`.
771 More generally:
773 .. math:: \operatorname{H}_{n,m} = \sum_{k=1}^{n} \frac{1}{k^m}
775 As `n \rightarrow \infty`, `\operatorname{H}_{n,m} \rightarrow \zeta(m)`,
776 the Riemann zeta function.
778 * ``harmonic(n)`` gives the nth harmonic number, `\operatorname{H}_n`
780 * ``harmonic(n, m)`` gives the nth generalized harmonic number
781 of order `m`, `\operatorname{H}_{n,m}`, where
782 ``harmonic(n) == harmonic(n, 1)``
784 This function can be extended to complex `n` and `m` where `n` is not a
785 negative integer or `m` is a nonpositive integer as
787 .. math:: \operatorname{H}_{n,m} = \begin{cases} \zeta(m) - \zeta(m, n+1)
788 & m \ne 1 \\ \psi(n+1) + \gamma & m = 1 \end{cases}
790 Examples
791 ========
793 >>> from sympy import harmonic, oo
795 >>> [harmonic(n) for n in range(6)]
796 [0, 1, 3/2, 11/6, 25/12, 137/60]
797 >>> [harmonic(n, 2) for n in range(6)]
798 [0, 1, 5/4, 49/36, 205/144, 5269/3600]
799 >>> harmonic(oo, 2)
800 pi**2/6
802 >>> from sympy import Symbol, Sum
803 >>> n = Symbol("n")
805 >>> harmonic(n).rewrite(Sum)
806 Sum(1/_k, (_k, 1, n))
808 We can evaluate harmonic numbers for all integral and positive
809 rational arguments:
811 >>> from sympy import S, expand_func, simplify
812 >>> harmonic(8)
813 761/280
814 >>> harmonic(11)
815 83711/27720
817 >>> H = harmonic(1/S(3))
818 >>> H
819 harmonic(1/3)
820 >>> He = expand_func(H)
821 >>> He
822 -log(6) - sqrt(3)*pi/6 + 2*Sum(log(sin(_k*pi/3))*cos(2*_k*pi/3), (_k, 1, 1))
823 + 3*Sum(1/(3*_k + 1), (_k, 0, 0))
824 >>> He.doit()
825 -log(6) - sqrt(3)*pi/6 - log(sqrt(3)/2) + 3
826 >>> H = harmonic(25/S(7))
827 >>> He = simplify(expand_func(H).doit())
828 >>> He
829 log(sin(2*pi/7)**(2*cos(16*pi/7))/(14*sin(pi/7)**(2*cos(pi/7))*cos(pi/14)**(2*sin(pi/14)))) + pi*tan(pi/14)/2 + 30247/9900
830 >>> He.n(40)
831 1.983697455232980674869851942390639915940
832 >>> harmonic(25/S(7)).n(40)
833 1.983697455232980674869851942390639915940
835 We can rewrite harmonic numbers in terms of polygamma functions:
837 >>> from sympy import digamma, polygamma
838 >>> m = Symbol("m", integer=True, positive=True)
840 >>> harmonic(n).rewrite(digamma)
841 polygamma(0, n + 1) + EulerGamma
843 >>> harmonic(n).rewrite(polygamma)
844 polygamma(0, n + 1) + EulerGamma
846 >>> harmonic(n,3).rewrite(polygamma)
847 polygamma(2, n + 1)/2 + zeta(3)
849 >>> simplify(harmonic(n,m).rewrite(polygamma))
850 Piecewise((polygamma(0, n + 1) + EulerGamma, Eq(m, 1)),
851 (-(-1)**m*polygamma(m - 1, n + 1)/factorial(m - 1) + zeta(m), True))
853 Integer offsets in the argument can be pulled out:
855 >>> from sympy import expand_func
857 >>> expand_func(harmonic(n+4))
858 harmonic(n) + 1/(n + 4) + 1/(n + 3) + 1/(n + 2) + 1/(n + 1)
860 >>> expand_func(harmonic(n-4))
861 harmonic(n) - 1/(n - 1) - 1/(n - 2) - 1/(n - 3) - 1/n
863 Some limits can be computed as well:
865 >>> from sympy import limit, oo
867 >>> limit(harmonic(n), n, oo)
868 oo
870 >>> limit(harmonic(n, 2), n, oo)
871 pi**2/6
873 >>> limit(harmonic(n, 3), n, oo)
874 zeta(3)
876 For `m > 1`, `H_{n,m}` tends to `\zeta(m)` in the limit of infinite `n`:
878 >>> m = Symbol("m", positive=True)
879 >>> limit(harmonic(n, m+1), n, oo)
880 zeta(m + 1)
882 See Also
883 ========
885 bell, bernoulli, catalan, euler, fibonacci, lucas, genocchi, partition, tribonacci
887 References
888 ==========
890 .. [1] https://en.wikipedia.org/wiki/Harmonic_number
891 .. [2] https://functions.wolfram.com/GammaBetaErf/HarmonicNumber/
892 .. [3] https://functions.wolfram.com/GammaBetaErf/HarmonicNumber2/
894 """
896 @classmethod
897 def eval(cls, n, m=None):
898 from sympy.functions.special.zeta_functions import zeta
899 if m is S.One:
900 return cls(n)
901 if m is None:
902 m = S.One
903 if n.is_zero:
904 return S.Zero
905 elif m.is_zero:
906 return n
907 elif n is S.Infinity:
908 if m.is_negative:
909 return S.NaN
910 elif is_le(m, S.One):
911 return S.Infinity
912 elif is_gt(m, S.One):
913 return zeta(m)
914 elif m.is_Integer and m.is_nonpositive:
915 return (bernoulli(1-m, n+1) - bernoulli(1-m)) / (1-m)
916 elif n.is_Integer:
917 if n.is_negative and (m.is_integer is False or m.is_nonpositive is False):
918 return S.ComplexInfinity if m is S.One else S.NaN
919 if n.is_nonnegative:
920 return Add(*(k**(-m) for k in range(1, int(n)+1)))
922 def _eval_rewrite_as_polygamma(self, n, m=S.One, **kwargs):
923 from sympy.functions.special.gamma_functions import gamma, polygamma
924 if m.is_integer and m.is_positive:
925 return Piecewise((polygamma(0, n+1) + S.EulerGamma, Eq(m, 1)),
926 (S.NegativeOne**m * (polygamma(m-1, 1) - polygamma(m-1, n+1)) /
927 gamma(m), True))
929 def _eval_rewrite_as_digamma(self, n, m=1, **kwargs):
930 from sympy.functions.special.gamma_functions import polygamma
931 return self.rewrite(polygamma)
933 def _eval_rewrite_as_trigamma(self, n, m=1, **kwargs):
934 from sympy.functions.special.gamma_functions import polygamma
935 return self.rewrite(polygamma)
937 def _eval_rewrite_as_Sum(self, n, m=None, **kwargs):
938 from sympy.concrete.summations import Sum
939 k = Dummy("k", integer=True)
940 if m is None:
941 m = S.One
942 return Sum(k**(-m), (k, 1, n))
944 def _eval_rewrite_as_zeta(self, n, m=S.One, **kwargs):
945 from sympy.functions.special.zeta_functions import zeta
946 from sympy.functions.special.gamma_functions import digamma
947 return Piecewise((digamma(n + 1) + S.EulerGamma, Eq(m, 1)),
948 (zeta(m) - zeta(m, n+1), True))
950 def _eval_expand_func(self, **hints):
951 from sympy.concrete.summations import Sum
952 n = self.args[0]
953 m = self.args[1] if len(self.args) == 2 else 1
955 if m == S.One:
956 if n.is_Add:
957 off = n.args[0]
958 nnew = n - off
959 if off.is_Integer and off.is_positive:
960 result = [S.One/(nnew + i) for i in range(off, 0, -1)] + [harmonic(nnew)]
961 return Add(*result)
962 elif off.is_Integer and off.is_negative:
963 result = [-S.One/(nnew + i) for i in range(0, off, -1)] + [harmonic(nnew)]
964 return Add(*result)
966 if n.is_Rational:
967 # Expansions for harmonic numbers at general rational arguments (u + p/q)
968 # Split n as u + p/q with p < q
969 p, q = n.as_numer_denom()
970 u = p // q
971 p = p - u * q
972 if u.is_nonnegative and p.is_positive and q.is_positive and p < q:
973 from sympy.functions.elementary.exponential import log
974 from sympy.functions.elementary.integers import floor
975 from sympy.functions.elementary.trigonometric import sin, cos, cot
976 k = Dummy("k")
977 t1 = q * Sum(1 / (q * k + p), (k, 0, u))
978 t2 = 2 * Sum(cos((2 * pi * p * k) / S(q)) *
979 log(sin((pi * k) / S(q))),
980 (k, 1, floor((q - 1) / S(2))))
981 t3 = (pi / 2) * cot((pi * p) / q) + log(2 * q)
982 return t1 + t2 - t3
984 return self
986 def _eval_rewrite_as_tractable(self, n, m=1, limitvar=None, **kwargs):
987 from sympy.functions.special.zeta_functions import zeta
988 from sympy.functions.special.gamma_functions import polygamma
989 pg = self.rewrite(polygamma)
990 if not isinstance(pg, harmonic):
991 return pg.rewrite("tractable", deep=True)
992 arg = m - S.One
993 if arg.is_nonzero:
994 return (zeta(m) - zeta(m, n+1)).rewrite("tractable", deep=True)
996 def _eval_evalf(self, prec):
997 if not all(x.is_number for x in self.args):
998 return
999 n = self.args[0]._to_mpmath(prec)
1000 m = (self.args[1] if len(self.args) > 1 else S.One)._to_mpmath(prec)
1001 if mp.isint(n) and n < 0:
1002 return S.NaN
1003 with workprec(prec):
1004 if m == 1:
1005 res = mp.harmonic(n)
1006 else:
1007 res = mp.zeta(m) - mp.zeta(m, n+1)
1008 return Expr._from_mpmath(res, prec)
1010 def fdiff(self, argindex=1):
1011 from sympy.functions.special.zeta_functions import zeta
1012 if len(self.args) == 2:
1013 n, m = self.args
1014 else:
1015 n, m = self.args + (1,)
1016 if argindex == 1:
1017 return m * zeta(m+1, n+1)
1018 else:
1019 raise ArgumentIndexError
1022#----------------------------------------------------------------------------#
1023# #
1024# Euler numbers #
1025# #
1026#----------------------------------------------------------------------------#
1029class euler(Function):
1030 r"""
1031 Euler numbers / Euler polynomials / Euler function
1033 The Euler numbers are given by:
1035 .. math:: E_{2n} = I \sum_{k=1}^{2n+1} \sum_{j=0}^k \binom{k}{j}
1036 \frac{(-1)^j (k-2j)^{2n+1}}{2^k I^k k}
1038 .. math:: E_{2n+1} = 0
1040 Euler numbers and Euler polynomials are related by
1042 .. math:: E_n = 2^n E_n\left(\frac{1}{2}\right).
1044 We compute symbolic Euler polynomials using Appell sequences,
1045 but numerical evaluation of the Euler polynomial is computed
1046 more efficiently (and more accurately) using the mpmath library.
1048 The Euler polynomials are special cases of the generalized Euler function,
1049 related to the Genocchi function as
1051 .. math:: \operatorname{E}(s, a) = -\frac{\operatorname{G}(s+1, a)}{s+1}
1053 with the limit of `\psi\left(\frac{a+1}{2}\right) - \psi\left(\frac{a}{2}\right)`
1054 being taken when `s = -1`. The (ordinary) Euler function interpolating
1055 the Euler numbers is then obtained as
1056 `\operatorname{E}(s) = 2^s \operatorname{E}\left(s, \frac{1}{2}\right)`.
1058 * ``euler(n)`` gives the nth Euler number `E_n`.
1059 * ``euler(s)`` gives the Euler function `\operatorname{E}(s)`.
1060 * ``euler(n, x)`` gives the nth Euler polynomial `E_n(x)`.
1061 * ``euler(s, a)`` gives the generalized Euler function `\operatorname{E}(s, a)`.
1063 Examples
1064 ========
1066 >>> from sympy import euler, Symbol, S
1067 >>> [euler(n) for n in range(10)]
1068 [1, 0, -1, 0, 5, 0, -61, 0, 1385, 0]
1069 >>> [2**n*euler(n,1) for n in range(10)]
1070 [1, 1, 0, -2, 0, 16, 0, -272, 0, 7936]
1071 >>> n = Symbol("n")
1072 >>> euler(n + 2*n)
1073 euler(3*n)
1075 >>> x = Symbol("x")
1076 >>> euler(n, x)
1077 euler(n, x)
1079 >>> euler(0, x)
1080 1
1081 >>> euler(1, x)
1082 x - 1/2
1083 >>> euler(2, x)
1084 x**2 - x
1085 >>> euler(3, x)
1086 x**3 - 3*x**2/2 + 1/4
1087 >>> euler(4, x)
1088 x**4 - 2*x**3 + x
1090 >>> euler(12, S.Half)
1091 2702765/4096
1092 >>> euler(12)
1093 2702765
1095 See Also
1096 ========
1098 andre, bell, bernoulli, catalan, fibonacci, harmonic, lucas, genocchi,
1099 partition, tribonacci, sympy.polys.appellseqs.euler_poly
1101 References
1102 ==========
1104 .. [1] https://en.wikipedia.org/wiki/Euler_numbers
1105 .. [2] https://mathworld.wolfram.com/EulerNumber.html
1106 .. [3] https://en.wikipedia.org/wiki/Alternating_permutation
1107 .. [4] https://mathworld.wolfram.com/AlternatingPermutation.html
1109 """
1111 @classmethod
1112 def eval(cls, n, x=None):
1113 if n.is_zero:
1114 return S.One
1115 elif n is S.NegativeOne:
1116 if x is None:
1117 return S.Pi/2
1118 from sympy.functions.special.gamma_functions import digamma
1119 return digamma((x+1)/2) - digamma(x/2)
1120 elif n.is_integer is False or n.is_nonnegative is False:
1121 return
1122 # Euler numbers
1123 elif x is None:
1124 if n.is_odd and n.is_positive:
1125 return S.Zero
1126 elif n.is_Number:
1127 from mpmath import mp
1128 n = n._to_mpmath(mp.prec)
1129 res = mp.eulernum(n, exact=True)
1130 return Integer(res)
1131 # Euler polynomials
1132 elif n.is_Number:
1133 return euler_poly(n, x)
1135 def _eval_rewrite_as_Sum(self, n, x=None, **kwargs):
1136 from sympy.concrete.summations import Sum
1137 if x is None and n.is_even:
1138 k = Dummy("k", integer=True)
1139 j = Dummy("j", integer=True)
1140 n = n / 2
1141 Em = (S.ImaginaryUnit * Sum(Sum(binomial(k, j) * (S.NegativeOne**j *
1142 (k - 2*j)**(2*n + 1)) /
1143 (2**k*S.ImaginaryUnit**k * k), (j, 0, k)), (k, 1, 2*n + 1)))
1144 return Em
1145 if x:
1146 k = Dummy("k", integer=True)
1147 return Sum(binomial(n, k)*euler(k)/2**k*(x - S.Half)**(n - k), (k, 0, n))
1149 def _eval_rewrite_as_genocchi(self, n, x=None, **kwargs):
1150 if x is None:
1151 return Piecewise((S.Pi/2, Eq(n, -1)),
1152 (-2**n * genocchi(n+1, S.Half) / (n+1), True))
1153 from sympy.functions.special.gamma_functions import digamma
1154 return Piecewise((digamma((x+1)/2) - digamma(x/2), Eq(n, -1)),
1155 (-genocchi(n+1, x) / (n+1), True))
1157 def _eval_evalf(self, prec):
1158 if not all(i.is_number for i in self.args):
1159 return
1160 from mpmath import mp
1161 m, x = (self.args[0], None) if len(self.args) == 1 else self.args
1162 m = m._to_mpmath(prec)
1163 if x is not None:
1164 x = x._to_mpmath(prec)
1165 with workprec(prec):
1166 if mp.isint(m) and m >= 0:
1167 res = mp.eulernum(m) if x is None else mp.eulerpoly(m, x)
1168 else:
1169 if m == -1:
1170 res = mp.pi if x is None else mp.digamma((x+1)/2) - mp.digamma(x/2)
1171 else:
1172 y = 0.5 if x is None else x
1173 res = 2 * (mp.zeta(-m, y) - 2**(m+1) * mp.zeta(-m, (y+1)/2))
1174 if x is None:
1175 res *= 2**m
1176 return Expr._from_mpmath(res, prec)
1179#----------------------------------------------------------------------------#
1180# #
1181# Catalan numbers #
1182# #
1183#----------------------------------------------------------------------------#
1186class catalan(Function):
1187 r"""
1188 Catalan numbers
1190 The `n^{th}` catalan number is given by:
1192 .. math :: C_n = \frac{1}{n+1} \binom{2n}{n}
1194 * ``catalan(n)`` gives the `n^{th}` Catalan number, `C_n`
1196 Examples
1197 ========
1199 >>> from sympy import (Symbol, binomial, gamma, hyper,
1200 ... catalan, diff, combsimp, Rational, I)
1202 >>> [catalan(i) for i in range(1,10)]
1203 [1, 2, 5, 14, 42, 132, 429, 1430, 4862]
1205 >>> n = Symbol("n", integer=True)
1207 >>> catalan(n)
1208 catalan(n)
1210 Catalan numbers can be transformed into several other, identical
1211 expressions involving other mathematical functions
1213 >>> catalan(n).rewrite(binomial)
1214 binomial(2*n, n)/(n + 1)
1216 >>> catalan(n).rewrite(gamma)
1217 4**n*gamma(n + 1/2)/(sqrt(pi)*gamma(n + 2))
1219 >>> catalan(n).rewrite(hyper)
1220 hyper((1 - n, -n), (2,), 1)
1222 For some non-integer values of n we can get closed form
1223 expressions by rewriting in terms of gamma functions:
1225 >>> catalan(Rational(1, 2)).rewrite(gamma)
1226 8/(3*pi)
1228 We can differentiate the Catalan numbers C(n) interpreted as a
1229 continuous real function in n:
1231 >>> diff(catalan(n), n)
1232 (polygamma(0, n + 1/2) - polygamma(0, n + 2) + log(4))*catalan(n)
1234 As a more advanced example consider the following ratio
1235 between consecutive numbers:
1237 >>> combsimp((catalan(n + 1)/catalan(n)).rewrite(binomial))
1238 2*(2*n + 1)/(n + 2)
1240 The Catalan numbers can be generalized to complex numbers:
1242 >>> catalan(I).rewrite(gamma)
1243 4**I*gamma(1/2 + I)/(sqrt(pi)*gamma(2 + I))
1245 and evaluated with arbitrary precision:
1247 >>> catalan(I).evalf(20)
1248 0.39764993382373624267 - 0.020884341620842555705*I
1250 See Also
1251 ========
1253 andre, bell, bernoulli, euler, fibonacci, harmonic, lucas, genocchi,
1254 partition, tribonacci, sympy.functions.combinatorial.factorials.binomial
1256 References
1257 ==========
1259 .. [1] https://en.wikipedia.org/wiki/Catalan_number
1260 .. [2] https://mathworld.wolfram.com/CatalanNumber.html
1261 .. [3] https://functions.wolfram.com/GammaBetaErf/CatalanNumber/
1262 .. [4] http://geometer.org/mathcircles/catalan.pdf
1264 """
1266 @classmethod
1267 def eval(cls, n):
1268 from sympy.functions.special.gamma_functions import gamma
1269 if (n.is_Integer and n.is_nonnegative) or \
1270 (n.is_noninteger and n.is_negative):
1271 return 4**n*gamma(n + S.Half)/(gamma(S.Half)*gamma(n + 2))
1273 if (n.is_integer and n.is_negative):
1274 if (n + 1).is_negative:
1275 return S.Zero
1276 if (n + 1).is_zero:
1277 return Rational(-1, 2)
1279 def fdiff(self, argindex=1):
1280 from sympy.functions.elementary.exponential import log
1281 from sympy.functions.special.gamma_functions import polygamma
1282 n = self.args[0]
1283 return catalan(n)*(polygamma(0, n + S.Half) - polygamma(0, n + 2) + log(4))
1285 def _eval_rewrite_as_binomial(self, n, **kwargs):
1286 return binomial(2*n, n)/(n + 1)
1288 def _eval_rewrite_as_factorial(self, n, **kwargs):
1289 return factorial(2*n) / (factorial(n+1) * factorial(n))
1291 def _eval_rewrite_as_gamma(self, n, piecewise=True, **kwargs):
1292 from sympy.functions.special.gamma_functions import gamma
1293 # The gamma function allows to generalize Catalan numbers to complex n
1294 return 4**n*gamma(n + S.Half)/(gamma(S.Half)*gamma(n + 2))
1296 def _eval_rewrite_as_hyper(self, n, **kwargs):
1297 from sympy.functions.special.hyper import hyper
1298 return hyper([1 - n, -n], [2], 1)
1300 def _eval_rewrite_as_Product(self, n, **kwargs):
1301 from sympy.concrete.products import Product
1302 if not (n.is_integer and n.is_nonnegative):
1303 return self
1304 k = Dummy('k', integer=True, positive=True)
1305 return Product((n + k) / k, (k, 2, n))
1307 def _eval_is_integer(self):
1308 if self.args[0].is_integer and self.args[0].is_nonnegative:
1309 return True
1311 def _eval_is_positive(self):
1312 if self.args[0].is_nonnegative:
1313 return True
1315 def _eval_is_composite(self):
1316 if self.args[0].is_integer and (self.args[0] - 3).is_positive:
1317 return True
1319 def _eval_evalf(self, prec):
1320 from sympy.functions.special.gamma_functions import gamma
1321 if self.args[0].is_number:
1322 return self.rewrite(gamma)._eval_evalf(prec)
1325#----------------------------------------------------------------------------#
1326# #
1327# Genocchi numbers #
1328# #
1329#----------------------------------------------------------------------------#
1332class genocchi(Function):
1333 r"""
1334 Genocchi numbers / Genocchi polynomials / Genocchi function
1336 The Genocchi numbers are a sequence of integers `G_n` that satisfy the
1337 relation:
1339 .. math:: \frac{-2t}{1 + e^{-t}} = \sum_{n=0}^\infty \frac{G_n t^n}{n!}
1341 They are related to the Bernoulli numbers by
1343 .. math:: G_n = 2 (1 - 2^n) B_n
1345 and generalize like the Bernoulli numbers to the Genocchi polynomials and
1346 function as
1348 .. math:: \operatorname{G}(s, a) = 2 \left(\operatorname{B}(s, a) -
1349 2^s \operatorname{B}\left(s, \frac{a+1}{2}\right)\right)
1351 .. versionchanged:: 1.12
1352 ``genocchi(1)`` gives `-1` instead of `1`.
1354 Examples
1355 ========
1357 >>> from sympy import genocchi, Symbol
1358 >>> [genocchi(n) for n in range(9)]
1359 [0, -1, -1, 0, 1, 0, -3, 0, 17]
1360 >>> n = Symbol('n', integer=True, positive=True)
1361 >>> genocchi(2*n + 1)
1362 0
1363 >>> x = Symbol('x')
1364 >>> genocchi(4, x)
1365 -4*x**3 + 6*x**2 - 1
1367 See Also
1368 ========
1370 bell, bernoulli, catalan, euler, fibonacci, harmonic, lucas, partition, tribonacci
1371 sympy.polys.appellseqs.genocchi_poly
1373 References
1374 ==========
1376 .. [1] https://en.wikipedia.org/wiki/Genocchi_number
1377 .. [2] https://mathworld.wolfram.com/GenocchiNumber.html
1378 .. [3] Peter Luschny, "An introduction to the Bernoulli function",
1379 https://arxiv.org/abs/2009.06743
1381 """
1383 @classmethod
1384 def eval(cls, n, x=None):
1385 if x is S.One:
1386 return cls(n)
1387 elif n.is_integer is False or n.is_nonnegative is False:
1388 return
1389 # Genocchi numbers
1390 elif x is None:
1391 if n.is_odd and (n-1).is_positive:
1392 return S.Zero
1393 elif n.is_Number:
1394 return 2 * (1-S(2)**n) * bernoulli(n)
1395 # Genocchi polynomials
1396 elif n.is_Number:
1397 return genocchi_poly(n, x)
1399 def _eval_rewrite_as_bernoulli(self, n, x=1, **kwargs):
1400 if x == 1 and n.is_integer and n.is_nonnegative:
1401 return 2 * (1-S(2)**n) * bernoulli(n)
1402 return 2 * (bernoulli(n, x) - 2**n * bernoulli(n, (x+1) / 2))
1404 def _eval_rewrite_as_dirichlet_eta(self, n, x=1, **kwargs):
1405 from sympy.functions.special.zeta_functions import dirichlet_eta
1406 return -2*n * dirichlet_eta(1-n, x)
1408 def _eval_is_integer(self):
1409 if len(self.args) > 1 and self.args[1] != 1:
1410 return
1411 n = self.args[0]
1412 if n.is_integer and n.is_nonnegative:
1413 return True
1415 def _eval_is_negative(self):
1416 if len(self.args) > 1 and self.args[1] != 1:
1417 return
1418 n = self.args[0]
1419 if n.is_integer and n.is_nonnegative:
1420 if n.is_odd:
1421 return fuzzy_not((n-1).is_positive)
1422 return (n/2).is_odd
1424 def _eval_is_positive(self):
1425 if len(self.args) > 1 and self.args[1] != 1:
1426 return
1427 n = self.args[0]
1428 if n.is_integer and n.is_nonnegative:
1429 if n.is_zero or n.is_odd:
1430 return False
1431 return (n/2).is_even
1433 def _eval_is_even(self):
1434 if len(self.args) > 1 and self.args[1] != 1:
1435 return
1436 n = self.args[0]
1437 if n.is_integer and n.is_nonnegative:
1438 if n.is_even:
1439 return n.is_zero
1440 return (n-1).is_positive
1442 def _eval_is_odd(self):
1443 if len(self.args) > 1 and self.args[1] != 1:
1444 return
1445 n = self.args[0]
1446 if n.is_integer and n.is_nonnegative:
1447 if n.is_even:
1448 return fuzzy_not(n.is_zero)
1449 return fuzzy_not((n-1).is_positive)
1451 def _eval_is_prime(self):
1452 if len(self.args) > 1 and self.args[1] != 1:
1453 return
1454 n = self.args[0]
1455 # only G_6 = -3 and G_8 = 17 are prime,
1456 # but SymPy does not consider negatives as prime
1457 # so only n=8 is tested
1458 return (n-8).is_zero
1460 def _eval_evalf(self, prec):
1461 if all(i.is_number for i in self.args):
1462 return self.rewrite(bernoulli)._eval_evalf(prec)
1465#----------------------------------------------------------------------------#
1466# #
1467# Andre numbers #
1468# #
1469#----------------------------------------------------------------------------#
1472class andre(Function):
1473 r"""
1474 Andre numbers / Andre function
1476 The Andre number `\mathcal{A}_n` is Luschny's name for half the number of
1477 *alternating permutations* on `n` elements, where a permutation is alternating
1478 if adjacent elements alternately compare "greater" and "smaller" going from
1479 left to right. For example, `2 < 3 > 1 < 4` is an alternating permutation.
1481 This sequence is A000111 in the OEIS, which assigns the names *up/down numbers*
1482 and *Euler zigzag numbers*. It satisfies a recurrence relation similar to that
1483 for the Catalan numbers, with `\mathcal{A}_0 = 1` and
1485 .. math:: 2 \mathcal{A}_{n+1} = \sum_{k=0}^n \binom{n}{k} \mathcal{A}_k \mathcal{A}_{n-k}
1487 The Bernoulli and Euler numbers are signed transformations of the odd- and
1488 even-indexed elements of this sequence respectively:
1490 .. math :: \operatorname{B}_{2k} = \frac{2k \mathcal{A}_{2k-1}}{(-4)^k - (-16)^k}
1492 .. math :: \operatorname{E}_{2k} = (-1)^k \mathcal{A}_{2k}
1494 Like the Bernoulli and Euler numbers, the Andre numbers are interpolated by the
1495 entire Andre function:
1497 .. math :: \mathcal{A}(s) = (-i)^{s+1} \operatorname{Li}_{-s}(i) +
1498 i^{s+1} \operatorname{Li}_{-s}(-i) = \\ \frac{2 \Gamma(s+1)}{(2\pi)^{s+1}}
1499 (\zeta(s+1, 1/4) - \zeta(s+1, 3/4) \cos{\pi s})
1501 Examples
1502 ========
1504 >>> from sympy import andre, euler, bernoulli
1505 >>> [andre(n) for n in range(11)]
1506 [1, 1, 1, 2, 5, 16, 61, 272, 1385, 7936, 50521]
1507 >>> [(-1)**k * andre(2*k) for k in range(7)]
1508 [1, -1, 5, -61, 1385, -50521, 2702765]
1509 >>> [euler(2*k) for k in range(7)]
1510 [1, -1, 5, -61, 1385, -50521, 2702765]
1511 >>> [andre(2*k-1) * (2*k) / ((-4)**k - (-16)**k) for k in range(1, 8)]
1512 [1/6, -1/30, 1/42, -1/30, 5/66, -691/2730, 7/6]
1513 >>> [bernoulli(2*k) for k in range(1, 8)]
1514 [1/6, -1/30, 1/42, -1/30, 5/66, -691/2730, 7/6]
1516 See Also
1517 ========
1519 bernoulli, catalan, euler, sympy.polys.appellseqs.andre_poly
1521 References
1522 ==========
1524 .. [1] https://en.wikipedia.org/wiki/Alternating_permutation
1525 .. [2] https://mathworld.wolfram.com/EulerZigzagNumber.html
1526 .. [3] Peter Luschny, "An introduction to the Bernoulli function",
1527 https://arxiv.org/abs/2009.06743
1528 """
1530 @classmethod
1531 def eval(cls, n):
1532 if n is S.NaN:
1533 return S.NaN
1534 elif n is S.Infinity:
1535 return S.Infinity
1536 if n.is_zero:
1537 return S.One
1538 elif n == -1:
1539 return -log(2)
1540 elif n == -2:
1541 return -2*S.Catalan
1542 elif n.is_Integer:
1543 if n.is_nonnegative and n.is_even:
1544 return abs(euler(n))
1545 elif n.is_odd:
1546 from sympy.functions.special.zeta_functions import zeta
1547 m = -n-1
1548 return I**m * Rational(1-2**m, 4**m) * zeta(-n)
1550 def _eval_rewrite_as_zeta(self, s, **kwargs):
1551 from sympy.functions.elementary.trigonometric import cos
1552 from sympy.functions.special.gamma_functions import gamma
1553 from sympy.functions.special.zeta_functions import zeta
1554 return 2 * gamma(s+1) / (2*pi)**(s+1) * \
1555 (zeta(s+1, S.One/4) - cos(pi*s) * zeta(s+1, S(3)/4))
1557 def _eval_rewrite_as_polylog(self, s, **kwargs):
1558 from sympy.functions.special.zeta_functions import polylog
1559 return (-I)**(s+1) * polylog(-s, I) + I**(s+1) * polylog(-s, -I)
1561 def _eval_is_integer(self):
1562 n = self.args[0]
1563 if n.is_integer and n.is_nonnegative:
1564 return True
1566 def _eval_is_positive(self):
1567 if self.args[0].is_nonnegative:
1568 return True
1570 def _eval_evalf(self, prec):
1571 if not self.args[0].is_number:
1572 return
1573 s = self.args[0]._to_mpmath(prec+12)
1574 with workprec(prec+12):
1575 sp, cp = mp.sinpi(s/2), mp.cospi(s/2)
1576 res = 2*mp.dirichlet(-s, (-sp, cp, sp, -cp))
1577 return Expr._from_mpmath(res, prec)
1580#----------------------------------------------------------------------------#
1581# #
1582# Partition numbers #
1583# #
1584#----------------------------------------------------------------------------#
1586_npartition = [1, 1]
1587class partition(Function):
1588 r"""
1589 Partition numbers
1591 The Partition numbers are a sequence of integers `p_n` that represent the
1592 number of distinct ways of representing `n` as a sum of natural numbers
1593 (with order irrelevant). The generating function for `p_n` is given by:
1595 .. math:: \sum_{n=0}^\infty p_n x^n = \prod_{k=1}^\infty (1 - x^k)^{-1}
1597 Examples
1598 ========
1600 >>> from sympy import partition, Symbol
1601 >>> [partition(n) for n in range(9)]
1602 [1, 1, 2, 3, 5, 7, 11, 15, 22]
1603 >>> n = Symbol('n', integer=True, negative=True)
1604 >>> partition(n)
1605 0
1607 See Also
1608 ========
1610 bell, bernoulli, catalan, euler, fibonacci, harmonic, lucas, genocchi, tribonacci
1612 References
1613 ==========
1615 .. [1] https://en.wikipedia.org/wiki/Partition_(number_theory%29
1616 .. [2] https://en.wikipedia.org/wiki/Pentagonal_number_theorem
1618 """
1620 @staticmethod
1621 def _partition(n):
1622 L = len(_npartition)
1623 if n < L:
1624 return _npartition[n]
1625 # lengthen cache
1626 for _n in range(L, n + 1):
1627 v, p, i = 0, 0, 0
1628 while 1:
1629 s = 0
1630 p += 3*i + 1 # p = pentagonal number: 1, 5, 12, ...
1631 if _n >= p:
1632 s += _npartition[_n - p]
1633 i += 1
1634 gp = p + i # gp = generalized pentagonal: 2, 7, 15, ...
1635 if _n >= gp:
1636 s += _npartition[_n - gp]
1637 if s == 0:
1638 break
1639 else:
1640 v += s if i%2 == 1 else -s
1641 _npartition.append(v)
1642 return v
1644 @classmethod
1645 def eval(cls, n):
1646 is_int = n.is_integer
1647 if is_int == False:
1648 raise ValueError("Partition numbers are defined only for "
1649 "integers")
1650 elif is_int:
1651 if n.is_negative:
1652 return S.Zero
1654 if n.is_zero or (n - 1).is_zero:
1655 return S.One
1657 if n.is_Integer:
1658 return Integer(cls._partition(n))
1661 def _eval_is_integer(self):
1662 if self.args[0].is_integer:
1663 return True
1665 def _eval_is_negative(self):
1666 if self.args[0].is_integer:
1667 return False
1669 def _eval_is_positive(self):
1670 n = self.args[0]
1671 if n.is_nonnegative and n.is_integer:
1672 return True
1675#######################################################################
1676###
1677### Functions for enumerating partitions, permutations and combinations
1678###
1679#######################################################################
1682class _MultisetHistogram(tuple):
1683 pass
1686_N = -1
1687_ITEMS = -2
1688_M = slice(None, _ITEMS)
1691def _multiset_histogram(n):
1692 """Return tuple used in permutation and combination counting. Input
1693 is a dictionary giving items with counts as values or a sequence of
1694 items (which need not be sorted).
1696 The data is stored in a class deriving from tuple so it is easily
1697 recognized and so it can be converted easily to a list.
1698 """
1699 if isinstance(n, dict): # item: count
1700 if not all(isinstance(v, int) and v >= 0 for v in n.values()):
1701 raise ValueError
1702 tot = sum(n.values())
1703 items = sum(1 for k in n if n[k] > 0)
1704 return _MultisetHistogram([n[k] for k in n if n[k] > 0] + [items, tot])
1705 else:
1706 n = list(n)
1707 s = set(n)
1708 lens = len(s)
1709 lenn = len(n)
1710 if lens == lenn:
1711 n = [1]*lenn + [lenn, lenn]
1712 return _MultisetHistogram(n)
1713 m = dict(zip(s, range(lens)))
1714 d = dict(zip(range(lens), (0,)*lens))
1715 for i in n:
1716 d[m[i]] += 1
1717 return _multiset_histogram(d)
1720def nP(n, k=None, replacement=False):
1721 """Return the number of permutations of ``n`` items taken ``k`` at a time.
1723 Possible values for ``n``:
1725 integer - set of length ``n``
1727 sequence - converted to a multiset internally
1729 multiset - {element: multiplicity}
1731 If ``k`` is None then the total of all permutations of length 0
1732 through the number of items represented by ``n`` will be returned.
1734 If ``replacement`` is True then a given item can appear more than once
1735 in the ``k`` items. (For example, for 'ab' permutations of 2 would
1736 include 'aa', 'ab', 'ba' and 'bb'.) The multiplicity of elements in
1737 ``n`` is ignored when ``replacement`` is True but the total number
1738 of elements is considered since no element can appear more times than
1739 the number of elements in ``n``.
1741 Examples
1742 ========
1744 >>> from sympy.functions.combinatorial.numbers import nP
1745 >>> from sympy.utilities.iterables import multiset_permutations, multiset
1746 >>> nP(3, 2)
1747 6
1748 >>> nP('abc', 2) == nP(multiset('abc'), 2) == 6
1749 True
1750 >>> nP('aab', 2)
1751 3
1752 >>> nP([1, 2, 2], 2)
1753 3
1754 >>> [nP(3, i) for i in range(4)]
1755 [1, 3, 6, 6]
1756 >>> nP(3) == sum(_)
1757 True
1759 When ``replacement`` is True, each item can have multiplicity
1760 equal to the length represented by ``n``:
1762 >>> nP('aabc', replacement=True)
1763 121
1764 >>> [len(list(multiset_permutations('aaaabbbbcccc', i))) for i in range(5)]
1765 [1, 3, 9, 27, 81]
1766 >>> sum(_)
1767 121
1769 See Also
1770 ========
1771 sympy.utilities.iterables.multiset_permutations
1773 References
1774 ==========
1776 .. [1] https://en.wikipedia.org/wiki/Permutation
1778 """
1779 try:
1780 n = as_int(n)
1781 except ValueError:
1782 return Integer(_nP(_multiset_histogram(n), k, replacement))
1783 return Integer(_nP(n, k, replacement))
1786@cacheit
1787def _nP(n, k=None, replacement=False):
1789 if k == 0:
1790 return 1
1791 if isinstance(n, SYMPY_INTS): # n different items
1792 # assert n >= 0
1793 if k is None:
1794 return sum(_nP(n, i, replacement) for i in range(n + 1))
1795 elif replacement:
1796 return n**k
1797 elif k > n:
1798 return 0
1799 elif k == n:
1800 return factorial(k)
1801 elif k == 1:
1802 return n
1803 else:
1804 # assert k >= 0
1805 return _product(n - k + 1, n)
1806 elif isinstance(n, _MultisetHistogram):
1807 if k is None:
1808 return sum(_nP(n, i, replacement) for i in range(n[_N] + 1))
1809 elif replacement:
1810 return n[_ITEMS]**k
1811 elif k == n[_N]:
1812 return factorial(k)/prod([factorial(i) for i in n[_M] if i > 1])
1813 elif k > n[_N]:
1814 return 0
1815 elif k == 1:
1816 return n[_ITEMS]
1817 else:
1818 # assert k >= 0
1819 tot = 0
1820 n = list(n)
1821 for i in range(len(n[_M])):
1822 if not n[i]:
1823 continue
1824 n[_N] -= 1
1825 if n[i] == 1:
1826 n[i] = 0
1827 n[_ITEMS] -= 1
1828 tot += _nP(_MultisetHistogram(n), k - 1)
1829 n[_ITEMS] += 1
1830 n[i] = 1
1831 else:
1832 n[i] -= 1
1833 tot += _nP(_MultisetHistogram(n), k - 1)
1834 n[i] += 1
1835 n[_N] += 1
1836 return tot
1839@cacheit
1840def _AOP_product(n):
1841 """for n = (m1, m2, .., mk) return the coefficients of the polynomial,
1842 prod(sum(x**i for i in range(nj + 1)) for nj in n); i.e. the coefficients
1843 of the product of AOPs (all-one polynomials) or order given in n. The
1844 resulting coefficient corresponding to x**r is the number of r-length
1845 combinations of sum(n) elements with multiplicities given in n.
1846 The coefficients are given as a default dictionary (so if a query is made
1847 for a key that is not present, 0 will be returned).
1849 Examples
1850 ========
1852 >>> from sympy.functions.combinatorial.numbers import _AOP_product
1853 >>> from sympy.abc import x
1854 >>> n = (2, 2, 3) # e.g. aabbccc
1855 >>> prod = ((x**2 + x + 1)*(x**2 + x + 1)*(x**3 + x**2 + x + 1)).expand()
1856 >>> c = _AOP_product(n); dict(c)
1857 {0: 1, 1: 3, 2: 6, 3: 8, 4: 8, 5: 6, 6: 3, 7: 1}
1858 >>> [c[i] for i in range(8)] == [prod.coeff(x, i) for i in range(8)]
1859 True
1861 The generating poly used here is the same as that listed in
1862 https://tinyurl.com/cep849r, but in a refactored form.
1864 """
1866 n = list(n)
1867 ord = sum(n)
1868 need = (ord + 2)//2
1869 rv = [1]*(n.pop() + 1)
1870 rv.extend((0,) * (need - len(rv)))
1871 rv = rv[:need]
1872 while n:
1873 ni = n.pop()
1874 N = ni + 1
1875 was = rv[:]
1876 for i in range(1, min(N, len(rv))):
1877 rv[i] += rv[i - 1]
1878 for i in range(N, need):
1879 rv[i] += rv[i - 1] - was[i - N]
1880 rev = list(reversed(rv))
1881 if ord % 2:
1882 rv = rv + rev
1883 else:
1884 rv[-1:] = rev
1885 d = defaultdict(int)
1886 for i, r in enumerate(rv):
1887 d[i] = r
1888 return d
1891def nC(n, k=None, replacement=False):
1892 """Return the number of combinations of ``n`` items taken ``k`` at a time.
1894 Possible values for ``n``:
1896 integer - set of length ``n``
1898 sequence - converted to a multiset internally
1900 multiset - {element: multiplicity}
1902 If ``k`` is None then the total of all combinations of length 0
1903 through the number of items represented in ``n`` will be returned.
1905 If ``replacement`` is True then a given item can appear more than once
1906 in the ``k`` items. (For example, for 'ab' sets of 2 would include 'aa',
1907 'ab', and 'bb'.) The multiplicity of elements in ``n`` is ignored when
1908 ``replacement`` is True but the total number of elements is considered
1909 since no element can appear more times than the number of elements in
1910 ``n``.
1912 Examples
1913 ========
1915 >>> from sympy.functions.combinatorial.numbers import nC
1916 >>> from sympy.utilities.iterables import multiset_combinations
1917 >>> nC(3, 2)
1918 3
1919 >>> nC('abc', 2)
1920 3
1921 >>> nC('aab', 2)
1922 2
1924 When ``replacement`` is True, each item can have multiplicity
1925 equal to the length represented by ``n``:
1927 >>> nC('aabc', replacement=True)
1928 35
1929 >>> [len(list(multiset_combinations('aaaabbbbcccc', i))) for i in range(5)]
1930 [1, 3, 6, 10, 15]
1931 >>> sum(_)
1932 35
1934 If there are ``k`` items with multiplicities ``m_1, m_2, ..., m_k``
1935 then the total of all combinations of length 0 through ``k`` is the
1936 product, ``(m_1 + 1)*(m_2 + 1)*...*(m_k + 1)``. When the multiplicity
1937 of each item is 1 (i.e., k unique items) then there are 2**k
1938 combinations. For example, if there are 4 unique items, the total number
1939 of combinations is 16:
1941 >>> sum(nC(4, i) for i in range(5))
1942 16
1944 See Also
1945 ========
1947 sympy.utilities.iterables.multiset_combinations
1949 References
1950 ==========
1952 .. [1] https://en.wikipedia.org/wiki/Combination
1953 .. [2] https://tinyurl.com/cep849r
1955 """
1957 if isinstance(n, SYMPY_INTS):
1958 if k is None:
1959 if not replacement:
1960 return 2**n
1961 return sum(nC(n, i, replacement) for i in range(n + 1))
1962 if k < 0:
1963 raise ValueError("k cannot be negative")
1964 if replacement:
1965 return binomial(n + k - 1, k)
1966 return binomial(n, k)
1967 if isinstance(n, _MultisetHistogram):
1968 N = n[_N]
1969 if k is None:
1970 if not replacement:
1971 return prod(m + 1 for m in n[_M])
1972 return sum(nC(n, i, replacement) for i in range(N + 1))
1973 elif replacement:
1974 return nC(n[_ITEMS], k, replacement)
1975 # assert k >= 0
1976 elif k in (1, N - 1):
1977 return n[_ITEMS]
1978 elif k in (0, N):
1979 return 1
1980 return _AOP_product(tuple(n[_M]))[k]
1981 else:
1982 return nC(_multiset_histogram(n), k, replacement)
1985def _eval_stirling1(n, k):
1986 if n == k == 0:
1987 return S.One
1988 if 0 in (n, k):
1989 return S.Zero
1991 # some special values
1992 if n == k:
1993 return S.One
1994 elif k == n - 1:
1995 return binomial(n, 2)
1996 elif k == n - 2:
1997 return (3*n - 1)*binomial(n, 3)/4
1998 elif k == n - 3:
1999 return binomial(n, 2)*binomial(n, 4)
2001 return _stirling1(n, k)
2004@cacheit
2005def _stirling1(n, k):
2006 row = [0, 1]+[0]*(k-1) # for n = 1
2007 for i in range(2, n+1):
2008 for j in range(min(k,i), 0, -1):
2009 row[j] = (i-1) * row[j] + row[j-1]
2010 return Integer(row[k])
2013def _eval_stirling2(n, k):
2014 if n == k == 0:
2015 return S.One
2016 if 0 in (n, k):
2017 return S.Zero
2019 # some special values
2020 if n == k:
2021 return S.One
2022 elif k == n - 1:
2023 return binomial(n, 2)
2024 elif k == 1:
2025 return S.One
2026 elif k == 2:
2027 return Integer(2**(n - 1) - 1)
2029 return _stirling2(n, k)
2032@cacheit
2033def _stirling2(n, k):
2034 row = [0, 1]+[0]*(k-1) # for n = 1
2035 for i in range(2, n+1):
2036 for j in range(min(k,i), 0, -1):
2037 row[j] = j * row[j] + row[j-1]
2038 return Integer(row[k])
2041def stirling(n, k, d=None, kind=2, signed=False):
2042 r"""Return Stirling number $S(n, k)$ of the first or second (default) kind.
2044 The sum of all Stirling numbers of the second kind for $k = 1$
2045 through $n$ is ``bell(n)``. The recurrence relationship for these numbers
2046 is:
2048 .. math :: {0 \brace 0} = 1; {n \brace 0} = {0 \brace k} = 0;
2050 .. math :: {{n+1} \brace k} = j {n \brace k} + {n \brace {k-1}}
2052 where $j$ is:
2053 $n$ for Stirling numbers of the first kind,
2054 $-n$ for signed Stirling numbers of the first kind,
2055 $k$ for Stirling numbers of the second kind.
2057 The first kind of Stirling number counts the number of permutations of
2058 ``n`` distinct items that have ``k`` cycles; the second kind counts the
2059 ways in which ``n`` distinct items can be partitioned into ``k`` parts.
2060 If ``d`` is given, the "reduced Stirling number of the second kind" is
2061 returned: $S^{d}(n, k) = S(n - d + 1, k - d + 1)$ with $n \ge k \ge d$.
2062 (This counts the ways to partition $n$ consecutive integers into $k$
2063 groups with no pairwise difference less than $d$. See example below.)
2065 To obtain the signed Stirling numbers of the first kind, use keyword
2066 ``signed=True``. Using this keyword automatically sets ``kind`` to 1.
2068 Examples
2069 ========
2071 >>> from sympy.functions.combinatorial.numbers import stirling, bell
2072 >>> from sympy.combinatorics import Permutation
2073 >>> from sympy.utilities.iterables import multiset_partitions, permutations
2075 First kind (unsigned by default):
2077 >>> [stirling(6, i, kind=1) for i in range(7)]
2078 [0, 120, 274, 225, 85, 15, 1]
2079 >>> perms = list(permutations(range(4)))
2080 >>> [sum(Permutation(p).cycles == i for p in perms) for i in range(5)]
2081 [0, 6, 11, 6, 1]
2082 >>> [stirling(4, i, kind=1) for i in range(5)]
2083 [0, 6, 11, 6, 1]
2085 First kind (signed):
2087 >>> [stirling(4, i, signed=True) for i in range(5)]
2088 [0, -6, 11, -6, 1]
2090 Second kind:
2092 >>> [stirling(10, i) for i in range(12)]
2093 [0, 1, 511, 9330, 34105, 42525, 22827, 5880, 750, 45, 1, 0]
2094 >>> sum(_) == bell(10)
2095 True
2096 >>> len(list(multiset_partitions(range(4), 2))) == stirling(4, 2)
2097 True
2099 Reduced second kind:
2101 >>> from sympy import subsets, oo
2102 >>> def delta(p):
2103 ... if len(p) == 1:
2104 ... return oo
2105 ... return min(abs(i[0] - i[1]) for i in subsets(p, 2))
2106 >>> parts = multiset_partitions(range(5), 3)
2107 >>> d = 2
2108 >>> sum(1 for p in parts if all(delta(i) >= d for i in p))
2109 7
2110 >>> stirling(5, 3, 2)
2111 7
2113 See Also
2114 ========
2115 sympy.utilities.iterables.multiset_partitions
2118 References
2119 ==========
2121 .. [1] https://en.wikipedia.org/wiki/Stirling_numbers_of_the_first_kind
2122 .. [2] https://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind
2124 """
2125 # TODO: make this a class like bell()
2127 n = as_int(n)
2128 k = as_int(k)
2129 if n < 0:
2130 raise ValueError('n must be nonnegative')
2131 if k > n:
2132 return S.Zero
2133 if d:
2134 # assert k >= d
2135 # kind is ignored -- only kind=2 is supported
2136 return _eval_stirling2(n - d + 1, k - d + 1)
2137 elif signed:
2138 # kind is ignored -- only kind=1 is supported
2139 return S.NegativeOne**(n - k)*_eval_stirling1(n, k)
2141 if kind == 1:
2142 return _eval_stirling1(n, k)
2143 elif kind == 2:
2144 return _eval_stirling2(n, k)
2145 else:
2146 raise ValueError('kind must be 1 or 2, not %s' % k)
2149@cacheit
2150def _nT(n, k):
2151 """Return the partitions of ``n`` items into ``k`` parts. This
2152 is used by ``nT`` for the case when ``n`` is an integer."""
2153 # really quick exits
2154 if k > n or k < 0:
2155 return 0
2156 if k in (1, n):
2157 return 1
2158 if k == 0:
2159 return 0
2160 # exits that could be done below but this is quicker
2161 if k == 2:
2162 return n//2
2163 d = n - k
2164 if d <= 3:
2165 return d
2166 # quick exit
2167 if 3*k >= n: # or, equivalently, 2*k >= d
2168 # all the information needed in this case
2169 # will be in the cache needed to calculate
2170 # partition(d), so...
2171 # update cache
2172 tot = partition._partition(d)
2173 # and correct for values not needed
2174 if d - k > 0:
2175 tot -= sum(_npartition[:d - k])
2176 return tot
2177 # regular exit
2178 # nT(n, k) = Sum(nT(n - k, m), (m, 1, k));
2179 # calculate needed nT(i, j) values
2180 p = [1]*d
2181 for i in range(2, k + 1):
2182 for m in range(i + 1, d):
2183 p[m] += p[m - i]
2184 d -= 1
2185 # if p[0] were appended to the end of p then the last
2186 # k values of p are the nT(n, j) values for 0 < j < k in reverse
2187 # order p[-1] = nT(n, 1), p[-2] = nT(n, 2), etc.... Instead of
2188 # putting the 1 from p[0] there, however, it is simply added to
2189 # the sum below which is valid for 1 < k <= n//2
2190 return (1 + sum(p[1 - k:]))
2193def nT(n, k=None):
2194 """Return the number of ``k``-sized partitions of ``n`` items.
2196 Possible values for ``n``:
2198 integer - ``n`` identical items
2200 sequence - converted to a multiset internally
2202 multiset - {element: multiplicity}
2204 Note: the convention for ``nT`` is different than that of ``nC`` and
2205 ``nP`` in that
2206 here an integer indicates ``n`` *identical* items instead of a set of
2207 length ``n``; this is in keeping with the ``partitions`` function which
2208 treats its integer-``n`` input like a list of ``n`` 1s. One can use
2209 ``range(n)`` for ``n`` to indicate ``n`` distinct items.
2211 If ``k`` is None then the total number of ways to partition the elements
2212 represented in ``n`` will be returned.
2214 Examples
2215 ========
2217 >>> from sympy.functions.combinatorial.numbers import nT
2219 Partitions of the given multiset:
2221 >>> [nT('aabbc', i) for i in range(1, 7)]
2222 [1, 8, 11, 5, 1, 0]
2223 >>> nT('aabbc') == sum(_)
2224 True
2226 >>> [nT("mississippi", i) for i in range(1, 12)]
2227 [1, 74, 609, 1521, 1768, 1224, 579, 197, 50, 9, 1]
2229 Partitions when all items are identical:
2231 >>> [nT(5, i) for i in range(1, 6)]
2232 [1, 2, 2, 1, 1]
2233 >>> nT('1'*5) == sum(_)
2234 True
2236 When all items are different:
2238 >>> [nT(range(5), i) for i in range(1, 6)]
2239 [1, 15, 25, 10, 1]
2240 >>> nT(range(5)) == sum(_)
2241 True
2243 Partitions of an integer expressed as a sum of positive integers:
2245 >>> from sympy import partition
2246 >>> partition(4)
2247 5
2248 >>> nT(4, 1) + nT(4, 2) + nT(4, 3) + nT(4, 4)
2249 5
2250 >>> nT('1'*4)
2251 5
2253 See Also
2254 ========
2255 sympy.utilities.iterables.partitions
2256 sympy.utilities.iterables.multiset_partitions
2257 sympy.functions.combinatorial.numbers.partition
2259 References
2260 ==========
2262 .. [1] https://web.archive.org/web/20210507012732/https://teaching.csse.uwa.edu.au/units/CITS7209/partition.pdf
2264 """
2266 if isinstance(n, SYMPY_INTS):
2267 # n identical items
2268 if k is None:
2269 return partition(n)
2270 if isinstance(k, SYMPY_INTS):
2271 n = as_int(n)
2272 k = as_int(k)
2273 return Integer(_nT(n, k))
2274 if not isinstance(n, _MultisetHistogram):
2275 try:
2276 # if n contains hashable items there is some
2277 # quick handling that can be done
2278 u = len(set(n))
2279 if u <= 1:
2280 return nT(len(n), k)
2281 elif u == len(n):
2282 n = range(u)
2283 raise TypeError
2284 except TypeError:
2285 n = _multiset_histogram(n)
2286 N = n[_N]
2287 if k is None and N == 1:
2288 return 1
2289 if k in (1, N):
2290 return 1
2291 if k == 2 or N == 2 and k is None:
2292 m, r = divmod(N, 2)
2293 rv = sum(nC(n, i) for i in range(1, m + 1))
2294 if not r:
2295 rv -= nC(n, m)//2
2296 if k is None:
2297 rv += 1 # for k == 1
2298 return rv
2299 if N == n[_ITEMS]:
2300 # all distinct
2301 if k is None:
2302 return bell(N)
2303 return stirling(N, k)
2304 m = MultisetPartitionTraverser()
2305 if k is None:
2306 return m.count_partitions(n[_M])
2307 # MultisetPartitionTraverser does not have a range-limited count
2308 # method, so need to enumerate and count
2309 tot = 0
2310 for discard in m.enum_range(n[_M], k-1, k):
2311 tot += 1
2312 return tot
2315#-----------------------------------------------------------------------------#
2316# #
2317# Motzkin numbers #
2318# #
2319#-----------------------------------------------------------------------------#
2322class motzkin(Function):
2323 """
2324 The nth Motzkin number is the number
2325 of ways of drawing non-intersecting chords
2326 between n points on a circle (not necessarily touching
2327 every point by a chord). The Motzkin numbers are named
2328 after Theodore Motzkin and have diverse applications
2329 in geometry, combinatorics and number theory.
2331 Motzkin numbers are the integer sequence defined by the
2332 initial terms `M_0 = 1`, `M_1 = 1` and the two-term recurrence relation
2333 `M_n = \frac{2*n + 1}{n + 2} * M_{n-1} + \frac{3n - 3}{n + 2} * M_{n-2}`.
2336 Examples
2337 ========
2339 >>> from sympy import motzkin
2341 >>> motzkin.is_motzkin(5)
2342 False
2343 >>> motzkin.find_motzkin_numbers_in_range(2,300)
2344 [2, 4, 9, 21, 51, 127]
2345 >>> motzkin.find_motzkin_numbers_in_range(2,900)
2346 [2, 4, 9, 21, 51, 127, 323, 835]
2347 >>> motzkin.find_first_n_motzkins(10)
2348 [1, 1, 2, 4, 9, 21, 51, 127, 323, 835]
2351 References
2352 ==========
2354 .. [1] https://en.wikipedia.org/wiki/Motzkin_number
2355 .. [2] https://mathworld.wolfram.com/MotzkinNumber.html
2357 """
2359 @staticmethod
2360 def is_motzkin(n):
2361 try:
2362 n = as_int(n)
2363 except ValueError:
2364 return False
2365 if n > 0:
2366 if n in (1, 2):
2367 return True
2369 tn1 = 1
2370 tn = 2
2371 i = 3
2372 while tn < n:
2373 a = ((2*i + 1)*tn + (3*i - 3)*tn1)/(i + 2)
2374 i += 1
2375 tn1 = tn
2376 tn = a
2378 if tn == n:
2379 return True
2380 else:
2381 return False
2383 else:
2384 return False
2386 @staticmethod
2387 def find_motzkin_numbers_in_range(x, y):
2388 if 0 <= x <= y:
2389 motzkins = []
2390 if x <= 1 <= y:
2391 motzkins.append(1)
2392 tn1 = 1
2393 tn = 2
2394 i = 3
2395 while tn <= y:
2396 if tn >= x:
2397 motzkins.append(tn)
2398 a = ((2*i + 1)*tn + (3*i - 3)*tn1)/(i + 2)
2399 i += 1
2400 tn1 = tn
2401 tn = int(a)
2403 return motzkins
2405 else:
2406 raise ValueError('The provided range is not valid. This condition should satisfy x <= y')
2408 @staticmethod
2409 def find_first_n_motzkins(n):
2410 try:
2411 n = as_int(n)
2412 except ValueError:
2413 raise ValueError('The provided number must be a positive integer')
2414 if n < 0:
2415 raise ValueError('The provided number must be a positive integer')
2416 motzkins = [1]
2417 if n >= 1:
2418 motzkins.append(1)
2419 tn1 = 1
2420 tn = 2
2421 i = 3
2422 while i <= n:
2423 motzkins.append(tn)
2424 a = ((2*i + 1)*tn + (3*i - 3)*tn1)/(i + 2)
2425 i += 1
2426 tn1 = tn
2427 tn = int(a)
2429 return motzkins
2431 @staticmethod
2432 @recurrence_memo([S.One, S.One])
2433 def _motzkin(n, prev):
2434 return ((2*n + 1)*prev[-1] + (3*n - 3)*prev[-2]) // (n + 2)
2436 @classmethod
2437 def eval(cls, n):
2438 try:
2439 n = as_int(n)
2440 except ValueError:
2441 raise ValueError('The provided number must be a positive integer')
2442 if n < 0:
2443 raise ValueError('The provided number must be a positive integer')
2444 return Integer(cls._motzkin(n - 1))
2447def nD(i=None, brute=None, *, n=None, m=None):
2448 """return the number of derangements for: ``n`` unique items, ``i``
2449 items (as a sequence or multiset), or multiplicities, ``m`` given
2450 as a sequence or multiset.
2452 Examples
2453 ========
2455 >>> from sympy.utilities.iterables import generate_derangements as enum
2456 >>> from sympy.functions.combinatorial.numbers import nD
2458 A derangement ``d`` of sequence ``s`` has all ``d[i] != s[i]``:
2460 >>> set([''.join(i) for i in enum('abc')])
2461 {'bca', 'cab'}
2462 >>> nD('abc')
2463 2
2465 Input as iterable or dictionary (multiset form) is accepted:
2467 >>> assert nD([1, 2, 2, 3, 3, 3]) == nD({1: 1, 2: 2, 3: 3})
2469 By default, a brute-force enumeration and count of multiset permutations
2470 is only done if there are fewer than 9 elements. There may be cases when
2471 there is high multiplicity with few unique elements that will benefit
2472 from a brute-force enumeration, too. For this reason, the `brute`
2473 keyword (default None) is provided. When False, the brute-force
2474 enumeration will never be used. When True, it will always be used.
2476 >>> nD('1111222233', brute=True)
2477 44
2479 For convenience, one may specify ``n`` distinct items using the
2480 ``n`` keyword:
2482 >>> assert nD(n=3) == nD('abc') == 2
2484 Since the number of derangments depends on the multiplicity of the
2485 elements and not the elements themselves, it may be more convenient
2486 to give a list or multiset of multiplicities using keyword ``m``:
2488 >>> assert nD('abc') == nD(m=(1,1,1)) == nD(m={1:3}) == 2
2490 """
2491 from sympy.integrals.integrals import integrate
2492 from sympy.functions.special.polynomials import laguerre
2493 from sympy.abc import x
2494 def ok(x):
2495 if not isinstance(x, SYMPY_INTS):
2496 raise TypeError('expecting integer values')
2497 if x < 0:
2498 raise ValueError('value must not be negative')
2499 return True
2501 if (i, n, m).count(None) != 2:
2502 raise ValueError('enter only 1 of i, n, or m')
2503 if i is not None:
2504 if isinstance(i, SYMPY_INTS):
2505 raise TypeError('items must be a list or dictionary')
2506 if not i:
2507 return S.Zero
2508 if type(i) is not dict:
2509 s = list(i)
2510 ms = multiset(s)
2511 elif type(i) is dict:
2512 all(ok(_) for _ in i.values())
2513 ms = {k: v for k, v in i.items() if v}
2514 s = None
2515 if not ms:
2516 return S.Zero
2517 N = sum(ms.values())
2518 counts = multiset(ms.values())
2519 nkey = len(ms)
2520 elif n is not None:
2521 ok(n)
2522 if not n:
2523 return S.Zero
2524 return subfactorial(n)
2525 elif m is not None:
2526 if isinstance(m, dict):
2527 all(ok(i) and ok(j) for i, j in m.items())
2528 counts = {k: v for k, v in m.items() if k*v}
2529 elif iterable(m) or isinstance(m, str):
2530 m = list(m)
2531 all(ok(i) for i in m)
2532 counts = multiset([i for i in m if i])
2533 else:
2534 raise TypeError('expecting iterable')
2535 if not counts:
2536 return S.Zero
2537 N = sum(k*v for k, v in counts.items())
2538 nkey = sum(counts.values())
2539 s = None
2540 big = int(max(counts))
2541 if big == 1: # no repetition
2542 return subfactorial(nkey)
2543 nval = len(counts)
2544 if big*2 > N:
2545 return S.Zero
2546 if big*2 == N:
2547 if nkey == 2 and nval == 1:
2548 return S.One # aaabbb
2549 if nkey - 1 == big: # one element repeated
2550 return factorial(big) # e.g. abc part of abcddd
2551 if N < 9 and brute is None or brute:
2552 # for all possibilities, this was found to be faster
2553 if s is None:
2554 s = []
2555 i = 0
2556 for m, v in counts.items():
2557 for j in range(v):
2558 s.extend([i]*m)
2559 i += 1
2560 return Integer(sum(1 for i in multiset_derangements(s)))
2561 from sympy.functions.elementary.exponential import exp
2562 return Integer(abs(integrate(exp(-x)*Mul(*[
2563 laguerre(i, x)**m for i, m in counts.items()]), (x, 0, oo))))