# Combinatorial#

This module implements various combinatorial functions.

## bell#

class sympy.functions.combinatorial.numbers.bell(n, k_sym=None, symbols=None)[source]#

Bell numbers / Bell polynomials

The Bell numbers satisfy $$B_0 = 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

## bernoulli#

class sympy.functions.combinatorial.numbers.bernoulli(n, sym=None)[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

## binomial#

class sympy.functions.combinatorial.factorials.binomial(n, k)[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{(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 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, 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

## catalan#

class sympy.functions.combinatorial.numbers.catalan(n)[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,
...     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

## euler#

class sympy.functions.combinatorial.numbers.euler(m, sym=None)[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 [R209]

$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 euler, Symbol, S
>>> [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

## factorial#

class sympy.functions.combinatorial.factorials.factorial(n)[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(arg)[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 [R211]

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

which is valid for non-negative integers $$x$$. The above formula is not very useful in case 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

## factorial2 / double factorial#

class sympy.functions.combinatorial.factorials.factorial2(arg)[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
>>> n = var('n')
>>> n
n
>>> factorial2(n + 1)
factorial2(n + 1)
>>> factorial2(5)
15
>>> factorial2(-1)
1
>>> factorial2(-5)
1/3


References

## FallingFactorial#

class sympy.functions.combinatorial.factorials.FallingFactorial(x, k)[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

$\texttt{ff(x, k)} = (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 [R213].

When $$x$$ is a $$~.Poly$$ instance of degree $$\ge 1$$ with single variable, $$(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

>>> from sympy import ff, Poly, Symbol
>>> from sympy.abc import x
>>> n = Symbol('n', 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)


Rewriting is complicated unless the relationship between the arguments is known, but falling factorial can be rewritten in terms of gamma, factorial and binomial and rising factorial.

>>> from sympy import factorial, rf, gamma, binomial, Symbol
>>> n = Symbol('n', integer=True, positive=True)
>>> F = ff(n, n - 2)
>>> for i in (rf, ff, factorial, binomial, gamma):
...  F.rewrite(i)
...
RisingFactorial(3, n - 2)
FallingFactorial(n, n - 2)
factorial(n)/2
binomial(n, n - 2)*factorial(n - 2)
gamma(n + 1)/2


References

[R214]

Peter Paule, “Greatest Factorial Factorization and Symbolic Summation”, Journal of Symbolic Computation, vol. 20, pp. 235-268, 1995.

## fibonacci#

class sympy.functions.combinatorial.numbers.fibonacci(n, sym=None)[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

## tribonacci#

class sympy.functions.combinatorial.numbers.tribonacci(n, sym=None)[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

## harmonic#

class sympy.functions.combinatorial.numbers.harmonic(n, m=None)[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(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
>>> 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 cannot compute the general relation yet:

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


which equals zeta(m) for m > 1.

References

## lucas#

class sympy.functions.combinatorial.numbers.lucas(n)[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

## genocchi#

class sympy.functions.combinatorial.numbers.genocchi(n)[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 genocchi, Symbol
>>> [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

## partition#

class sympy.functions.combinatorial.numbers.partition(n)[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 partition, Symbol
>>> [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

## MultiFactorial#

class sympy.functions.combinatorial.factorials.MultiFactorial(*args)[source]#

## RisingFactorial#

class sympy.functions.combinatorial.factorials.RisingFactorial(x, k)[source]#

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

$\texttt{rf(y, k)} = (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 $$\ge 1$$ with a single variable, $$(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 [R231].

Examples

>>> from sympy import rf, Poly
>>> from sympy.abc import x
>>> 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')


Rewriting is complicated unless the relationship between the arguments is known, but rising factorial can be rewritten in terms of gamma, factorial, binomial, and falling factorial.

>>> from sympy import Symbol, factorial, ff, binomial, gamma
>>> n = Symbol('n', integer=True, positive=True)
>>> R = rf(n, n + 2)
>>> for i in (rf, ff, factorial, binomial, gamma):
...  R.rewrite(i)
...
RisingFactorial(n, n + 2)
FallingFactorial(2*n + 1, n + 2)
factorial(2*n + 1)/factorial(n - 1)
binomial(2*n + 1, n + 2)*factorial(n + 2)
gamma(2*n + 2)/gamma(n)


References

[R231] (1,2)

Peter Paule, “Greatest Factorial Factorization and Symbolic Summation”, Journal of Symbolic Computation, vol. 20, pp. 235-268, 1995.

## 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 \ge k \ge 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 - i) 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

# 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#

sympy.functions.combinatorial.numbers.nC(n, k=None, replacement=False)[source]#

Return the number of combinations of n items taken k at a time.

Possible values for n:

integer - set of length n

sequence - converted to a multiset internally

multiset - {element: multiplicity}

If k is None then the total of all combinations of length 0 through the number of items represented in n will be returned.

If replacement is True then a given item can appear more than once in the k items. (For example, for ‘ab’ sets of 2 would include ‘aa’, ‘ab’, and ‘bb’.) The multiplicity of elements in n is ignored when replacement is True but the total number of elements is considered since no element can appear more times than the number of elements in n.

Examples

>>> from sympy.functions.combinatorial.numbers import nC
>>> from sympy.utilities.iterables import multiset_combinations
>>> nC(3, 2)
3
>>> nC('abc', 2)
3
>>> nC('aab', 2)
2


When replacement is True, each item can have multiplicity equal to the length represented by n:

>>> nC('aabc', replacement=True)
35
>>> [len(list(multiset_combinations('aaaabbbbcccc', i))) for i in range(5)]
[1, 3, 6, 10, 15]
>>> sum(_)
35


If there are k items with multiplicities m_1, m_2, ..., m_k then the total of all combinations of length 0 through k is the product, (m_1 + 1)*(m_2 + 1)*...*(m_k + 1). When the multiplicity of each item is 1 (i.e., k unique items) then there are 2**k combinations. For example, if there are 4 unique items, the total number of combinations is 16:

>>> sum(nC(4, i) for i in range(5))
16


References

## nP#

sympy.functions.combinatorial.numbers.nP(n, k=None, replacement=False)[source]#

Return the number of permutations of n items taken k at a time.

Possible values for n:

integer - set of length n

sequence - converted to a multiset internally

multiset - {element: multiplicity}

If k is None then the total of all permutations of length 0 through the number of items represented by n will be returned.

If replacement is True then a given item can appear more than once in the k items. (For example, for ‘ab’ permutations of 2 would include ‘aa’, ‘ab’, ‘ba’ and ‘bb’.) The multiplicity of elements in n is ignored when replacement is True but the total number of elements is considered since no element can appear more times than the number of elements in n.

Examples

>>> from sympy.functions.combinatorial.numbers import nP
>>> from sympy.utilities.iterables import multiset_permutations, multiset
>>> nP(3, 2)
6
>>> nP('abc', 2) == nP(multiset('abc'), 2) == 6
True
>>> nP('aab', 2)
3
>>> nP([1, 2, 2], 2)
3
>>> [nP(3, i) for i in range(4)]
[1, 3, 6, 6]
>>> nP(3) == sum(_)
True


When replacement is True, each item can have multiplicity equal to the length represented by n:

>>> nP('aabc', replacement=True)
121
>>> [len(list(multiset_permutations('aaaabbbbcccc', i))) for i in range(5)]
[1, 3, 9, 27, 81]
>>> sum(_)
121


References

## nT#

sympy.functions.combinatorial.numbers.nT(n, k=None)[source]#

Return the number of k-sized partitions of n items.

Possible values for n:

integer - n identical items

sequence - converted to a multiset internally

multiset - {element: multiplicity}

Note: the convention for nT is different than that of nC and nP in that here an integer indicates n identical items instead of a set of length n; this is in keeping with the partitions function which treats its integer-n input like a list of n 1s. One can use range(n) for n to indicate n distinct items.

If k is None then the total number of ways to partition the elements represented in n will be returned.

Examples

>>> from sympy.functions.combinatorial.numbers import nT


Partitions of the given multiset:

>>> [nT('aabbc', i) for i in range(1, 7)]
[1, 8, 11, 5, 1, 0]
>>> nT('aabbc') == sum(_)
True

>>> [nT("mississippi", i) for i in range(1, 12)]
[1, 74, 609, 1521, 1768, 1224, 579, 197, 50, 9, 1]


Partitions when all items are identical:

>>> [nT(5, i) for i in range(1, 6)]
[1, 2, 2, 1, 1]
>>> nT('1'*5) == sum(_)
True


When all items are different:

>>> [nT(range(5), i) for i in range(1, 6)]
[1, 15, 25, 10, 1]
>>> nT(range(5)) == sum(_)
True


Partitions of an integer expressed as a sum of positive integers:

>>> from sympy import partition
>>> partition(4)
5
>>> nT(4, 1) + nT(4, 2) + nT(4, 3) + nT(4, 4)
5
>>> nT('1'*4)
5


References