Combinatorial

This module implements various combinatorial functions.

bell

class sympy.functions.combinatorial.numbers.bell[source]

Bell numbers / Bell polynomials

The Bell numbers satisfy B0=1 and

B_n = \sum_{k=0}^{n-1} \binom{n-1}{k} B_k.

They are also given by:

B_n = \frac{1}{e} \sum_{k=0}^{\infty} \frac{k^n}{k!}.

The Bell polynomials are given by B_0(x) = 1 and

B_n(x) = x \sum_{k=1}^{n-1} \binom{n-1}{k-1} B_{k-1}(x).

The second kind of Bell polynomials (are sometimes called “partial” Bell polynomials or incomplete Bell polynomials) are defined as

B_{n,k}(x_1, x_2,\dotsc x_{n-k+1}) = \sum_{j_1+j_2+j_2+\dotsb=k \atop j_1+2j_2+3j_2+\dotsb=n} \frac{n!}{j_1!j_2!\dotsb j_{n-k+1}!} \left(\frac{x_1}{1!} \right)^{j_1} \left(\frac{x_2}{2!} \right)^{j_2} \dotsb \left(\frac{x_{n-k+1}}{(n-k+1)!} \right) ^{j_{n-k+1}}.
  • bell(n) gives the n^{th} Bell number, B_n.

  • bell(n, x) gives the n^{th} Bell polynomial, B_n(x).

  • bell(n, k, (x1, x2, ...)) gives Bell polynomials of the second kind, B_{n,k}(x_1, x_2, \dotsc, x_{n-k+1}).

Notes

Not to be confused with Bernoulli numbers and Bernoulli polynomials, which use the same notation.

Examples

>>> from sympy import bell, Symbol, symbols
>>> [bell(n) for n in range(11)]
[1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975]
>>> bell(30)
846749014511809332450147
>>> bell(4, Symbol('t'))
t**4 + 6*t**3 + 7*t**2 + t
>>> bell(6, 2, symbols('x:6')[1:])
6*x1*x5 + 15*x2*x4 + 10*x3**2

References

R138

https://en.wikipedia.org/wiki/Bell_number

R139

http://mathworld.wolfram.com/BellNumber.html

R140

http://mathworld.wolfram.com/BellPolynomial.html

bernoulli

class sympy.functions.combinatorial.numbers.bernoulli[source]

Bernoulli numbers / Bernoulli polynomials

The Bernoulli numbers are a sequence of rational numbers defined by B_0 = 1 and the recursive relation (n > 0):

0 = \sum_{k=0}^n \binom{n+1}{k} B_k

They are also commonly defined by their exponential generating function, which is \frac{x}{e^x - 1}. For odd indices > 1, the Bernoulli numbers are zero.

The Bernoulli polynomials satisfy the analogous formula:

B_n(x) = \sum_{k=0}^n \binom{n}{k} B_k x^{n-k}

Bernoulli numbers and Bernoulli polynomials are related as B_n(0) = B_n.

We compute Bernoulli numbers using Ramanujan’s formula:

B_n = \frac{A(n) - S(n)}{\binom{n+3}{n}}

where:

\begin{split}A(n) = \begin{cases} \frac{n+3}{3} & n \equiv 0\ \text{or}\ 2 \pmod{6} \\ -\frac{n+3}{6} & n \equiv 4 \pmod{6} \end{cases}\end{split}

and:

S(n) = \sum_{k=1}^{[n/6]} \binom{n+3}{n-6k} B_{n-6k}

This formula is similar to the sum given in the definition, but cuts 2/3 of the terms. For Bernoulli polynomials, we use the formula in the definition.

  • bernoulli(n) gives the nth Bernoulli number, B_n

  • bernoulli(n, x) gives the nth Bernoulli polynomial in x, B_n(x)

Examples

>>> from sympy import bernoulli
>>> [bernoulli(n) for n in range(11)]
[1, -1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0, 5/66]
>>> bernoulli(1000001)
0

References

R141

https://en.wikipedia.org/wiki/Bernoulli_number

R142

https://en.wikipedia.org/wiki/Bernoulli_polynomial

R143

http://mathworld.wolfram.com/BernoulliNumber.html

R144

http://mathworld.wolfram.com/BernoulliPolynomial.html

binomial

class sympy.functions.combinatorial.factorials.binomial[source]

Implementation of the binomial coefficient. It can be defined in two ways depending on its desired interpretation:

\binom{n}{k} = \frac{n!}{k!(n-k)!}\ \text{or}\ \binom{n}{k} = \frac{ff(n, k)}{k!}

First, in a strict combinatorial sense it defines the number of ways we can choose k elements from a set of n elements. In this case both arguments are nonnegative integers and binomial is computed using an efficient algorithm based on prime factorization.

The other definition is generalization for arbitrary n, however k must also be nonnegative. This case is very useful when evaluating summations.

For the sake of convenience for negative integer k this function will return zero no matter what valued is the other argument.

To expand the binomial when n is a symbol, use either expand_func() or expand(func=True). The former will keep the polynomial in factored form while the latter will expand the polynomial itself. See examples for details.

Examples

>>> from sympy import Symbol, Rational, binomial, expand_func
>>> n = Symbol('n', integer=True, positive=True)
>>> binomial(15, 8)
6435
>>> binomial(n, -1)
0

Rows of Pascal’s triangle can be generated with the binomial function:

>>> for N in range(8):
...     print([binomial(N, i) for i in range(N + 1)])
...
[1]
[1, 1]
[1, 2, 1]
[1, 3, 3, 1]
[1, 4, 6, 4, 1]
[1, 5, 10, 10, 5, 1]
[1, 6, 15, 20, 15, 6, 1]
[1, 7, 21, 35, 35, 21, 7, 1]

As can a given diagonal, e.g. the 4th diagonal:

>>> N = -4
>>> [binomial(N, i) for i in range(1 - N)]
[1, -4, 10, -20, 35]
>>> binomial(Rational(5, 4), 3)
-5/128
>>> binomial(Rational(-5, 4), 3)
-195/128
>>> binomial(n, 3)
binomial(n, 3)
>>> binomial(n, 3).expand(func=True)
n**3/6 - n**2/2 + n/3
>>> expand_func(binomial(n, 3))
n*(n - 2)*(n - 1)/6

References

R145

https://www.johndcook.com/blog/binomial_coefficients/

catalan

class sympy.functions.combinatorial.numbers.catalan[source]

Catalan numbers

The n^{th} catalan number is given by:

C_n = \frac{1}{n+1} \binom{2n}{n}
  • catalan(n) gives the n^{th} Catalan number, C_n

Examples

>>> from sympy import (Symbol, binomial, gamma, hyper, polygamma,
...             catalan, diff, combsimp, Rational, I)
>>> [catalan(i) for i in range(1,10)]
[1, 2, 5, 14, 42, 132, 429, 1430, 4862]
>>> n = Symbol("n", integer=True)
>>> catalan(n)
catalan(n)

Catalan numbers can be transformed into several other, identical expressions involving other mathematical functions

>>> catalan(n).rewrite(binomial)
binomial(2*n, n)/(n + 1)
>>> catalan(n).rewrite(gamma)
4**n*gamma(n + 1/2)/(sqrt(pi)*gamma(n + 2))
>>> catalan(n).rewrite(hyper)
hyper((1 - n, -n), (2,), 1)

For some non-integer values of n we can get closed form expressions by rewriting in terms of gamma functions:

>>> catalan(Rational(1,2)).rewrite(gamma)
8/(3*pi)

We can differentiate the Catalan numbers C(n) interpreted as a continuous real function in n:

>>> diff(catalan(n), n)
(polygamma(0, n + 1/2) - polygamma(0, n + 2) + log(4))*catalan(n)

As a more advanced example consider the following ratio between consecutive numbers:

>>> combsimp((catalan(n + 1)/catalan(n)).rewrite(binomial))
2*(2*n + 1)/(n + 2)

The Catalan numbers can be generalized to complex numbers:

>>> catalan(I).rewrite(gamma)
4**I*gamma(1/2 + I)/(sqrt(pi)*gamma(2 + I))

and evaluated with arbitrary precision:

>>> catalan(I).evalf(20)
0.39764993382373624267 - 0.020884341620842555705*I

References

R146

https://en.wikipedia.org/wiki/Catalan_number

R147

http://mathworld.wolfram.com/CatalanNumber.html

R148

http://functions.wolfram.com/GammaBetaErf/CatalanNumber/

R149

http://geometer.org/mathcircles/catalan.pdf

euler

class sympy.functions.combinatorial.numbers.euler[source]

Euler numbers / Euler polynomials

The Euler numbers are given by:

E_{2n} = I \sum_{k=1}^{2n+1} \sum_{j=0}^k \binom{k}{j} \frac{(-1)^j (k-2j)^{2n+1}}{2^k I^k k}
E_{2n+1} = 0

Euler numbers and Euler polynomials are related by

E_n = 2^n E_n\left(\frac{1}{2}\right).

We compute symbolic Euler polynomials using [R154]

E_n(x) = \sum_{k=0}^n \binom{n}{k} \frac{E_k}{2^k} \left(x - \frac{1}{2}\right)^{n-k}.

However, numerical evaluation of the Euler polynomial is computed more efficiently (and more accurately) using the mpmath library.

  • euler(n) gives the n^{th} Euler number, E_n.

  • euler(n, x) gives the n^{th} Euler polynomial, E_n(x).

Examples

>>> from sympy import Symbol, S
>>> from sympy.functions import euler
>>> [euler(n) for n in range(10)]
[1, 0, -1, 0, 5, 0, -61, 0, 1385, 0]
>>> n = Symbol("n")
>>> euler(n + 2*n)
euler(3*n)
>>> x = Symbol("x")
>>> euler(n, x)
euler(n, x)
>>> euler(0, x)
1
>>> euler(1, x)
x - 1/2
>>> euler(2, x)
x**2 - x
>>> euler(3, x)
x**3 - 3*x**2/2 + 1/4
>>> euler(4, x)
x**4 - 2*x**3 + x
>>> euler(12, S.Half)
2702765/4096
>>> euler(12)
2702765

References

R150

https://en.wikipedia.org/wiki/Euler_numbers

R151

http://mathworld.wolfram.com/EulerNumber.html

R152

https://en.wikipedia.org/wiki/Alternating_permutation

R153

http://mathworld.wolfram.com/AlternatingPermutation.html

R154(1,2)

http://dlmf.nist.gov/24.2#ii

factorial

class sympy.functions.combinatorial.factorials.factorial[source]

Implementation of factorial function over nonnegative integers. By convention (consistent with the gamma function and the binomial coefficients), factorial of a negative integer is complex infinity.

The factorial is very important in combinatorics where it gives the number of ways in which n objects can be permuted. It also arises in calculus, probability, number theory, etc.

There is strict relation of factorial with gamma function. In fact n! = gamma(n+1) for nonnegative integers. Rewrite of this kind is very useful in case of combinatorial simplification.

Computation of the factorial is done using two algorithms. For small arguments a precomputed look up table is used. However for bigger input algorithm Prime-Swing is used. It is the fastest algorithm known and computes n! via prime factorization of special class of numbers, called here the ‘Swing Numbers’.

Examples

>>> from sympy import Symbol, factorial, S
>>> n = Symbol('n', integer=True)
>>> factorial(0)
1
>>> factorial(7)
5040
>>> factorial(-2)
zoo
>>> factorial(n)
factorial(n)
>>> factorial(2*n)
factorial(2*n)
>>> factorial(S(1)/2)
factorial(1/2)

subfactorial

class sympy.functions.combinatorial.factorials.subfactorial[source]

The subfactorial counts the derangements of n items and is defined for non-negative integers as:

\begin{split}!n = \begin{cases} 1 & n = 0 \\ 0 & n = 1 \\ (n-1)(!(n-1) + !(n-2)) & n > 1 \end{cases}\end{split}

It can also be written as int(round(n!/exp(1))) but the recursive definition with caching is implemented for this function.

An interesting analytic expression is the following [R156]

!x = \Gamma(x + 1, -1)/e

which is valid for non-negative integers x. The above formula is not very useful incase of non-integers. \Gamma(x + 1, -1) is single-valued only for integral arguments x, elsewhere on the positive real axis it has an infinite number of branches none of which are real.

Examples

>>> from sympy import subfactorial
>>> from sympy.abc import n
>>> subfactorial(n + 1)
subfactorial(n + 1)
>>> subfactorial(5)
44

References

R155

https://en.wikipedia.org/wiki/Subfactorial

R156(1,2)

http://mathworld.wolfram.com/Subfactorial.html

factorial2 / double factorial

class sympy.functions.combinatorial.factorials.factorial2[source]

The double factorial n!!, not to be confused with (n!)!

The double factorial is defined for nonnegative integers and for odd negative integers as:

\begin{split}n!! = \begin{cases} 1 & n = 0 \\ n(n-2)(n-4) \cdots 1 & n\ \text{positive odd} \\ n(n-2)(n-4) \cdots 2 & n\ \text{positive even} \\ (n+2)!!/(n+2) & n\ \text{negative odd} \end{cases}\end{split}

Examples

>>> from sympy import factorial2, var
>>> var('n')
n
>>> factorial2(n + 1)
factorial2(n + 1)
>>> factorial2(5)
15
>>> factorial2(-1)
1
>>> factorial2(-5)
1/3

References

R157

https://en.wikipedia.org/wiki/Double_factorial

FallingFactorial

class sympy.functions.combinatorial.factorials.FallingFactorial[source]

Falling factorial (related to rising factorial) is a double valued function arising in concrete mathematics, hypergeometric functions and series expansions. It is defined by

ff(x,k) = x \cdot (x-1) \cdots (x-k+1)

where x can be arbitrary expression and k is an integer. For more information check “Concrete mathematics” by Graham, pp. 66 or visit http://mathworld.wolfram.com/FallingFactorial.html page.

When x is a Poly instance of degree >= 1 with single variable, ff(x,k) = x(y) \cdot x(y-1) \cdots x(y-k+1), where y is the variable of x. This is as described in Peter Paule, “Greatest Factorial Factorization and Symbolic Summation”, Journal of Symbolic Computation, vol. 20, pp. 235-268, 1995.

>>> from sympy import ff, factorial, rf, gamma, polygamma, binomial, symbols, Poly
>>> from sympy.abc import x, k
>>> n, m = symbols('n m', integer=True)
>>> ff(x, 0)
1
>>> ff(5, 5)
120
>>> ff(x, 5) == x*(x-1)*(x-2)*(x-3)*(x-4)
True
>>> ff(Poly(x**2, x), 2)
Poly(x**4 - 2*x**3 + x**2, x, domain='ZZ')
>>> ff(n, n)
factorial(n)

Rewrite

>>> ff(x, k).rewrite(gamma)
(-1)**k*gamma(k - x)/gamma(-x)
>>> ff(x, k).rewrite(rf)
RisingFactorial(-k + x + 1, k)
>>> ff(x, m).rewrite(binomial)
binomial(x, m)*factorial(m)
>>> ff(n, m).rewrite(factorial)
factorial(n)/factorial(-m + n)

References

R158

http://mathworld.wolfram.com/FallingFactorial.html

fibonacci

class sympy.functions.combinatorial.numbers.fibonacci[source]

Fibonacci numbers / Fibonacci polynomials

The Fibonacci numbers are the integer sequence defined by the initial terms F_0 = 0, F_1 = 1 and the two-term recurrence relation F_n = F_{n-1} + F_{n-2}. This definition extended to arbitrary real and complex arguments using the formula

F_z = \frac{\phi^z - \cos(\pi z) \phi^{-z}}{\sqrt 5}

The Fibonacci polynomials are defined by F_1(x) = 1, F_2(x) = x, and F_n(x) = x*F_{n-1}(x) + F_{n-2}(x) for n > 2. For all positive integers n, F_n(1) = F_n.

  • fibonacci(n) gives the n^{th} Fibonacci number, F_n

  • fibonacci(n, x) gives the n^{th} Fibonacci polynomial in x, F_n(x)

Examples

>>> from sympy import fibonacci, Symbol
>>> [fibonacci(x) for x in range(11)]
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
>>> fibonacci(5, Symbol('t'))
t**4 + 3*t**2 + 1

References

R159

https://en.wikipedia.org/wiki/Fibonacci_number

R160

http://mathworld.wolfram.com/FibonacciNumber.html

tribonacci

class sympy.functions.combinatorial.numbers.tribonacci[source]

Tribonacci numbers / Tribonacci polynomials

The Tribonacci numbers are the integer sequence defined by the initial terms T_0 = 0, T_1 = 1, T_2 = 1 and the three-term recurrence relation T_n = T_{n-1} + T_{n-2} + T_{n-3}.

The Tribonacci polynomials are defined by T_0(x) = 0, T_1(x) = 1, T_2(x) = x^2, and T_n(x) = x^2 T_{n-1}(x) + x T_{n-2}(x) + T_{n-3}(x) for n > 2. For all positive integers n, T_n(1) = T_n.

  • tribonacci(n) gives the n^{th} Tribonacci number, T_n

  • tribonacci(n, x) gives the n^{th} Tribonacci polynomial in x, T_n(x)

Examples

>>> from sympy import tribonacci, Symbol
>>> [tribonacci(x) for x in range(11)]
[0, 1, 1, 2, 4, 7, 13, 24, 44, 81, 149]
>>> tribonacci(5, Symbol('t'))
t**8 + 3*t**5 + 3*t**2

References

R161

https://en.wikipedia.org/wiki/Generalizations_of_Fibonacci_numbers#Tribonacci_numbers

R162

http://mathworld.wolfram.com/TribonacciNumber.html

R163

https://oeis.org/A000073

harmonic

class sympy.functions.combinatorial.numbers.harmonic[source]

Harmonic numbers

The nth harmonic number is given by \operatorname{H}_{n} = 1 + \frac{1}{2} + \frac{1}{3} + \ldots + \frac{1}{n}.

More generally:

\operatorname{H}_{n,m} = \sum_{k=1}^{n} \frac{1}{k^m}

As n \rightarrow \infty, \operatorname{H}_{n,m} \rightarrow \zeta(m), the Riemann zeta function.

  • harmonic(n) gives the nth harmonic number, \operatorname{H}_n

  • harmonic(n, m) gives the nth generalized harmonic number of order m, \operatorname{H}_{n,m}, where harmonic(n) == harmonic(n, 1)

Examples

>>> from sympy import harmonic, oo
>>> [harmonic(n) for n in range(6)]
[0, 1, 3/2, 11/6, 25/12, 137/60]
>>> [harmonic(n, 2) for n in range(6)]
[0, 1, 5/4, 49/36, 205/144, 5269/3600]
>>> harmonic(oo, 2)
pi**2/6
>>> from sympy import Symbol, Sum
>>> n = Symbol("n")
>>> harmonic(n).rewrite(Sum)
Sum(1/_k, (_k, 1, n))

We can evaluate harmonic numbers for all integral and positive rational arguments:

>>> from sympy import S, expand_func, simplify
>>> harmonic(8)
761/280
>>> harmonic(11)
83711/27720
>>> H = harmonic(1/S(3))
>>> H
harmonic(1/3)
>>> He = expand_func(H)
>>> He
-log(6) - sqrt(3)*pi/6 + 2*Sum(log(sin(_k*pi/3))*cos(2*_k*pi/3), (_k, 1, 1))
                       + 3*Sum(1/(3*_k + 1), (_k, 0, 0))
>>> He.doit()
-log(6) - sqrt(3)*pi/6 - log(sqrt(3)/2) + 3
>>> H = harmonic(25/S(7))
>>> He = simplify(expand_func(H).doit())
>>> He
log(sin(pi/7)**(-2*cos(pi/7))*sin(2*pi/7)**(2*cos(16*pi/7))*cos(pi/14)**(-2*sin(pi/14))/14)
+ pi*tan(pi/14)/2 + 30247/9900
>>> He.n(40)
1.983697455232980674869851942390639915940
>>> harmonic(25/S(7)).n(40)
1.983697455232980674869851942390639915940

We can rewrite harmonic numbers in terms of polygamma functions:

>>> from sympy import digamma, polygamma
>>> m = Symbol("m")
>>> harmonic(n).rewrite(digamma)
polygamma(0, n + 1) + EulerGamma
>>> harmonic(n).rewrite(polygamma)
polygamma(0, n + 1) + EulerGamma
>>> harmonic(n,3).rewrite(polygamma)
polygamma(2, n + 1)/2 - polygamma(2, 1)/2
>>> harmonic(n,m).rewrite(polygamma)
(-1)**m*(polygamma(m - 1, 1) - polygamma(m - 1, n + 1))/factorial(m - 1)

Integer offsets in the argument can be pulled out:

>>> from sympy import expand_func
>>> expand_func(harmonic(n+4))
harmonic(n) + 1/(n + 4) + 1/(n + 3) + 1/(n + 2) + 1/(n + 1)
>>> expand_func(harmonic(n-4))
harmonic(n) - 1/(n - 1) - 1/(n - 2) - 1/(n - 3) - 1/n

Some limits can be computed as well:

>>> from sympy import limit, oo
>>> limit(harmonic(n), n, oo)
oo
>>> limit(harmonic(n, 2), n, oo)
pi**2/6
>>> limit(harmonic(n, 3), n, oo)
-polygamma(2, 1)/2

However we can not compute the general relation yet:

>>> limit(harmonic(n, m), n, oo)
harmonic(oo, m)

which equals zeta(m) for m > 1.

References

R164

https://en.wikipedia.org/wiki/Harmonic_number

R165

http://functions.wolfram.com/GammaBetaErf/HarmonicNumber/

R166

http://functions.wolfram.com/GammaBetaErf/HarmonicNumber2/

lucas

class sympy.functions.combinatorial.numbers.lucas[source]

Lucas numbers

Lucas numbers satisfy a recurrence relation similar to that of the Fibonacci sequence, in which each term is the sum of the preceding two. They are generated by choosing the initial values L_0 = 2 and L_1 = 1.

  • lucas(n) gives the n^{th} Lucas number

Examples

>>> from sympy import lucas
>>> [lucas(x) for x in range(11)]
[2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123]

References

R167

https://en.wikipedia.org/wiki/Lucas_number

R168

http://mathworld.wolfram.com/LucasNumber.html

genocchi

class sympy.functions.combinatorial.numbers.genocchi[source]

Genocchi numbers

The Genocchi numbers are a sequence of integers G_n that satisfy the relation:

\frac{2t}{e^t + 1} = \sum_{n=1}^\infty \frac{G_n t^n}{n!}

Examples

>>> from sympy import Symbol
>>> from sympy.functions import genocchi
>>> [genocchi(n) for n in range(1, 9)]
[1, -1, 0, 1, 0, -3, 0, 17]
>>> n = Symbol('n', integer=True, positive=True)
>>> genocchi(2*n + 1)
0

References

R169

https://en.wikipedia.org/wiki/Genocchi_number

R170

http://mathworld.wolfram.com/GenocchiNumber.html

partition

class sympy.functions.combinatorial.numbers.partition[source]

Partition numbers

The Partition numbers are a sequence of integers p_n that represent the number of distinct ways of representing n as a sum of natural numbers (with order irrelevant). The generating function for p_n is given by:

\sum_{n=0}^\infty p_n x^n = \prod_{k=1}^\infty (1 - x^k)^{-1}

Examples

>>> from sympy import Symbol
>>> from sympy.functions import partition
>>> [partition(n) for n in range(9)]
[1, 1, 2, 3, 5, 7, 11, 15, 22]
>>> n = Symbol('n', integer=True, negative=True)
>>> partition(n)
0

References

R171

https://en.wikipedia.org/wiki/Partition_(number_theory%29

R172

https://en.wikipedia.org/wiki/Pentagonal_number_theorem

MultiFactorial

class sympy.functions.combinatorial.factorials.MultiFactorial[source]

RisingFactorial

class sympy.functions.combinatorial.factorials.RisingFactorial[source]

Rising factorial (also called Pochhammer symbol) is a double valued function arising in concrete mathematics, hypergeometric functions and series expansions. It is defined by:

rf(x,k) = x \cdot (x+1) \cdots (x+k-1)

where x can be arbitrary expression and k is an integer. For more information check “Concrete mathematics” by Graham, pp. 66 or visit http://mathworld.wolfram.com/RisingFactorial.html page.

When x is a Poly instance of degree >= 1 with a single variable, rf(x,k) = x(y) \cdot x(y+1) \cdots x(y+k-1), where y is the variable of x. This is as described in Peter Paule, “Greatest Factorial Factorization and Symbolic Summation”, Journal of Symbolic Computation, vol. 20, pp. 235-268, 1995.

Examples

>>> from sympy import rf, symbols, factorial, ff, binomial, Poly
>>> from sympy.abc import x
>>> n, k = symbols('n k', integer=True)
>>> rf(x, 0)
1
>>> rf(1, 5)
120
>>> rf(x, 5) == x*(1 + x)*(2 + x)*(3 + x)*(4 + x)
True
>>> rf(Poly(x**3, x), 2)
Poly(x**6 + 3*x**5 + 3*x**4 + x**3, x, domain='ZZ')

Rewrite

>>> rf(x, k).rewrite(ff)
FallingFactorial(k + x - 1, k)
>>> rf(x, k).rewrite(binomial)
binomial(k + x - 1, k)*factorial(k)
>>> rf(n, k).rewrite(factorial)
factorial(k + n - 1)/factorial(n - 1)

References

R173

https://en.wikipedia.org/wiki/Pochhammer_symbol

stirling

sympy.functions.combinatorial.numbers.stirling(n, k, d=None, kind=2, signed=False)[source]

Return Stirling number S(n, k) of the first or second (default) kind.

The sum of all Stirling numbers of the second kind for k = 1 through n is bell(n). The recurrence relationship for these numbers is:

{0 \brace 0} = 1; {n \brace 0} = {0 \brace k} = 0;
{{n+1} \brace k} = j {n \brace k} + {n \brace {k-1}}
where j is:

n for Stirling numbers of the first kind -n for signed Stirling numbers of the first kind k for Stirling numbers of the second kind

The first kind of Stirling number counts the number of permutations of n distinct items that have k cycles; the second kind counts the ways in which n distinct items can be partitioned into k parts. If d is given, the “reduced Stirling number of the second kind” is returned: S^{d}(n, k) = S(n - d + 1, k - d + 1) with n >= k >= d. (This counts the ways to partition n consecutive integers into k groups with no pairwise difference less than d. See example below.)

To obtain the signed Stirling numbers of the first kind, use keyword signed=True. Using this keyword automatically sets kind to 1.

Examples

>>> from sympy.functions.combinatorial.numbers import stirling, bell
>>> from sympy.combinatorics import Permutation
>>> from sympy.utilities.iterables import multiset_partitions, permutations

First kind (unsigned by default):

>>> [stirling(6, i, kind=1) for i in range(7)]
[0, 120, 274, 225, 85, 15, 1]
>>> perms = list(permutations(range(4)))
>>> [sum(Permutation(p).cycles == i for p in perms) for i in range(5)]
[0, 6, 11, 6, 1]
>>> [stirling(4, i, kind=1) for i in range(5)]
[0, 6, 11, 6, 1]

First kind (signed):

>>> [stirling(4, i, signed=True) for i in range(5)]
[0, -6, 11, -6, 1]

Second kind:

>>> [stirling(10, i) for i in range(12)]
[0, 1, 511, 9330, 34105, 42525, 22827, 5880, 750, 45, 1, 0]
>>> sum(_) == bell(10)
True
>>> len(list(multiset_partitions(range(4), 2))) == stirling(4, 2)
True

Reduced second kind:

>>> from sympy import subsets, oo
>>> def delta(p):
...    if len(p) == 1:
...        return oo
...    return min(abs(i[0] - i[1]) for i in subsets(p, 2))
>>> parts = multiset_partitions(range(5), 3)
>>> d = 2
>>> sum(1 for p in parts if all(delta(i) >= d for i in p))
7
>>> stirling(5, 3, 2)
7

References

R174

https://en.wikipedia.org/wiki/Stirling_numbers_of_the_first_kind

R175

https://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind

Enumeration

Three functions are available. Each of them attempts to efficiently compute a given combinatorial quantity for a given set or multiset which can be entered as an integer, sequence or multiset (dictionary with elements as keys and multiplicities as values). The k parameter indicates the number of elements to pick (or the number of partitions to make). When k is None, the sum of the enumeration for all k (from 0 through the number of items represented by n) is returned. A replacement parameter is recognized for combinations and permutations; this indicates that any item may appear with multiplicity as high as the number of items in the original set.

>>> from sympy.functions.combinatorial.numbers import nC, nP, nT
>>> items = 'baby'

nC

Calculate the number of combinations of length k.

>>> [nC(items, k) for k in range(len(items) + 1)], nC(items)
([1, 3, 4, 3, 1], 12)
>>> nC('aaa', 2)
1
>>> nC('abc', 2)
3
>>> nC(3, 2)
3

nP

Calculate the number of permutations of length k.

>>> [nP(items, k) for k in range(len(items) + 1)], nP(items)
([1, 3, 7, 12, 12], 35)
>>> nC('aaa', 2)
1
>>> nC('abc', 2)
3
>>> nC(3, 2)
3

nT

Calculate the number of partitions that have k parts.

>>> [nT(items, k) for k in range(len(items) + 1)], nT(items)
([0, 1, 5, 4, 1], 11)
>>> nT('aaa', 2)
1
>>> nT('abc', 2)
3
>>> nT(3, 2)
1

Note that the integer for n indicates identical items for nT but indicates n different items for nC and nP.