Combinatorial

This module implements various combinatorial functions.

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

See also

bernoulli, catalan, euler, fibonacci, harmonic, lucas, genocchi, partition, tribonacci

References

class sympy.functions.combinatorial.numbers.bernoulli(n, x=None)[source]

Bernoulli numbers / Bernoulli polynomials / Bernoulli function

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

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

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

The Bernoulli polynomials satisfy the analogous formula:

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

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

The generalized Bernoulli function \(\operatorname{B}(s, a)\) is defined for any complex \(s\) and \(a\), except where \(a\) is a nonpositive integer and \(s\) is not a nonnegative integer. It is an entire function of \(s\) for fixed \(a\), related to the Hurwitz zeta function by

\[\begin{split}\operatorname{B}(s, a) = \begin{cases} -s \zeta(1-s, a) & s \ne 0 \\ 1 & s = 0 \end{cases}\end{split}\]

When \(s\) is a nonnegative integer this function reduces to the Bernoulli polynomials: \(\operatorname{B}(n, x) = B_n(x)\). When \(a\) is omitted it is assumed to be 1, yielding the (ordinary) Bernoulli function which interpolates the Bernoulli numbers and is related to the Riemann zeta function.

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 \(\frac{2}{3}\) of the terms. For Bernoulli polynomials, we use Appell sequences.

For \(n\) a nonnegative integer and \(s\), \(a\), \(x\) arbitrary complex numbers,

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

  • bernoulli(s) gives the Bernoulli function \(\operatorname{B}(s)\)

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

  • bernoulli(s, a) gives the generalized Bernoulli function \(\operatorname{B}(s, a)\)

Changed in version 1.12: bernoulli(1) gives \(+\frac{1}{2}\) instead of \(-\frac{1}{2}\). This choice of value confers several theoretical advantages [R214], including the extension to complex parameters described above which this function now implements. The previous behavior, defined only for nonnegative integers \(n\), can be obtained with (-1)**n*bernoulli(n).

Examples

>>> from sympy import bernoulli
>>> from sympy.abc import x
>>> [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
>>> bernoulli(3, x)
x**3 - 3*x**2/2 + x/2

See also

andre, bell, catalan, euler, fibonacci, harmonic, lucas, genocchi, partition, tribonacci, sympy.polys.appellseqs.bernoulli_poly

References

[R214] (1,2)

Peter Luschny, “The Bernoulli Manifesto”, https://luschny.de/math/zeta/The-Bernoulli-Manifesto.html

[R215]

Peter Luschny, “An introduction to the Bernoulli function”, https://arxiv.org/abs/2009.06743

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]
[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

In many cases, we can also compute binomial coefficients modulo a prime p quickly using Lucas’ Theorem [R217], though we need to include \(evaluate=False\) to postpone evaluation:

>>> from sympy import Mod
>>> Mod(binomial(156675, 4433, evaluate=False), 10**5 + 3)
28625

Using a generalisation of Lucas’s Theorem given by Granville [R218], we can extend this to arbitrary n:

>>> Mod(binomial(10**18, 10**12, evaluate=False), (10**5 + 3)**2)
3744312326

References

[R218] (1,2)

Binomial coefficients modulo prime powers, Andrew Granville, Available: https://web.archive.org/web/20170202003812/http://www.dms.umontreal.ca/~andrew/PDF/BinCoeff.pdf

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((-n, 1 - 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

See also

andre, bell, bernoulli, euler, fibonacci, harmonic, lucas, genocchi, partition, tribonacci, sympy.functions.combinatorial.factorials.binomial

References

class sympy.functions.combinatorial.numbers.euler(n, x=None)[source]

Euler numbers / Euler polynomials / Euler function

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 Appell sequences, but numerical evaluation of the Euler polynomial is computed more efficiently (and more accurately) using the mpmath library.

The Euler polynomials are special cases of the generalized Euler function, related to the Genocchi function as

\[\operatorname{E}(s, a) = -\frac{\operatorname{G}(s+1, a)}{s+1}\]

with the limit of \(\psi\left(\frac{a+1}{2}\right) - \psi\left(\frac{a}{2}\right)\) being taken when \(s = -1\). The (ordinary) Euler function interpolating the Euler numbers is then obtained as \(\operatorname{E}(s) = 2^s \operatorname{E}\left(s, \frac{1}{2}\right)\).

  • euler(n) gives the nth Euler number \(E_n\).

  • euler(s) gives the Euler function \(\operatorname{E}(s)\).

  • euler(n, x) gives the nth Euler polynomial \(E_n(x)\).

  • euler(s, a) gives the generalized Euler function \(\operatorname{E}(s, a)\).

Examples

>>> from sympy import euler, Symbol, S
>>> [euler(n) for n in range(10)]
[1, 0, -1, 0, 5, 0, -61, 0, 1385, 0]
>>> [2**n*euler(n,1) for n in range(10)]
[1, 1, 0, -2, 0, 16, 0, -272, 0, 7936]
>>> 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

See also

andre, bell, bernoulli, catalan, fibonacci, harmonic, lucas, genocchi, partition, tribonacci, sympy.polys.appellseqs.euler_poly

References

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)

See also

factorial2, RisingFactorial, FallingFactorial

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 [R228]

\[!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

See also

factorial, uppergamma, sympy.utilities.iterables.generate_derangements

References

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

See also

factorial, RisingFactorial, FallingFactorial

References

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 [R230].

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

See also

factorial, factorial2, RisingFactorial

References

[R231]

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

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

See also

bell, bernoulli, catalan, euler, harmonic, lucas, genocchi, partition, tribonacci

References

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

See also

bell, bernoulli, catalan, euler, fibonacci, harmonic, lucas, genocchi, partition

References

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)

This function can be extended to complex \(n\) and \(m\) where \(n\) is not a negative integer or \(m\) is a nonpositive integer as

\[\begin{split}\operatorname{H}_{n,m} = \begin{cases} \zeta(m) - \zeta(m, n+1) & m \ne 1 \\ \psi(n+1) + \gamma & m = 1 \end{cases}\end{split}\]

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", integer=True, positive=True)

>>> 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 + zeta(3)

>>> simplify(harmonic(n,m).rewrite(polygamma))
Piecewise((polygamma(0, n + 1) + EulerGamma, Eq(m, 1)),
(-(-1)**m*polygamma(m - 1, n + 1)/factorial(m - 1) + zeta(m), True))

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)
zeta(3)

For \(m > 1\), \(H_{n,m}\) tends to \(\zeta(m)\) in the limit of infinite \(n\):

>>> m = Symbol("m", positive=True)
>>> limit(harmonic(n, m+1), n, oo)
zeta(m + 1)

See also

bell, bernoulli, catalan, euler, fibonacci, lucas, genocchi, partition, tribonacci

References

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]

See also

bell, bernoulli, catalan, euler, fibonacci, harmonic, genocchi, partition, tribonacci

References

class sympy.functions.combinatorial.numbers.genocchi(n, x=None)[source]

Genocchi numbers / Genocchi polynomials / Genocchi function

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

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

They are related to the Bernoulli numbers by

\[G_n = 2 (1 - 2^n) B_n\]

and generalize like the Bernoulli numbers to the Genocchi polynomials and function as

\[\operatorname{G}(s, a) = 2 \left(\operatorname{B}(s, a) - 2^s \operatorname{B}\left(s, \frac{a+1}{2}\right)\right)\]

Changed in version 1.12: genocchi(1) gives \(-1\) instead of \(1\).

Examples

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

See also

bell, bernoulli, catalan, euler, fibonacci, harmonic, lucas, partition, tribonacci, sympy.polys.appellseqs.genocchi_poly

References

[R245]

Peter Luschny, “An introduction to the Bernoulli function”, https://arxiv.org/abs/2009.06743

class sympy.functions.combinatorial.numbers.andre(n)[source]

Andre numbers / Andre function

The Andre number \(\mathcal{A}_n\) is Luschny’s name for half the number of alternating permutations on \(n\) elements, where a permutation is alternating if adjacent elements alternately compare “greater” and “smaller” going from left to right. For example, \(2 < 3 > 1 < 4\) is an alternating permutation.

This sequence is A000111 in the OEIS, which assigns the names up/down numbers and Euler zigzag numbers. It satisfies a recurrence relation similar to that for the Catalan numbers, with \(\mathcal{A}_0 = 1\) and

\[2 \mathcal{A}_{n+1} = \sum_{k=0}^n \binom{n}{k} \mathcal{A}_k \mathcal{A}_{n-k}\]

The Bernoulli and Euler numbers are signed transformations of the odd- and even-indexed elements of this sequence respectively:

\[\operatorname{B}_{2k} = \frac{2k \mathcal{A}_{2k-1}}{(-4)^k - (-16)^k}\]
\[\operatorname{E}_{2k} = (-1)^k \mathcal{A}_{2k}\]

Like the Bernoulli and Euler numbers, the Andre numbers are interpolated by the entire Andre function:

\[\begin{split}\mathcal{A}(s) = (-i)^{s+1} \operatorname{Li}_{-s}(i) + i^{s+1} \operatorname{Li}_{-s}(-i) = \\ \frac{2 \Gamma(s+1)}{(2\pi)^{s+1}} (\zeta(s+1, 1/4) - \zeta(s+1, 3/4) \cos{\pi s})\end{split}\]

Examples

>>> from sympy import andre, euler, bernoulli
>>> [andre(n) for n in range(11)]
[1, 1, 1, 2, 5, 16, 61, 272, 1385, 7936, 50521]
>>> [(-1)**k * andre(2*k) for k in range(7)]
[1, -1, 5, -61, 1385, -50521, 2702765]
>>> [euler(2*k) for k in range(7)]
[1, -1, 5, -61, 1385, -50521, 2702765]
>>> [andre(2*k-1) * (2*k) / ((-4)**k - (-16)**k) for k in range(1, 8)]
[1/6, -1/30, 1/42, -1/30, 5/66, -691/2730, 7/6]
>>> [bernoulli(2*k) for k in range(1, 8)]
[1/6, -1/30, 1/42, -1/30, 5/66, -691/2730, 7/6]

See also

bernoulli, catalan, euler, sympy.polys.appellseqs.andre_poly

References

[R248]

Peter Luschny, “An introduction to the Bernoulli function”, https://arxiv.org/abs/2009.06743

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

See also

bell, bernoulli, catalan, euler, fibonacci, harmonic, lucas, genocchi, tribonacci

References

class sympy.functions.combinatorial.numbers.divisor_sigma(n, k=1)[source]

Calculate the divisor function \(\sigma_k(n)\) for positive integer n

divisor_sigma(n, k) is equal to sum([x**k for x in divisors(n)])

If n’s prime factorization is:

\[n = \prod_{i=1}^\omega p_i^{m_i},\]

then

\[\sigma_k(n) = \prod_{i=1}^\omega (1+p_i^k+p_i^{2k}+\cdots + p_i^{m_ik}).\]

Examples

>>> from sympy.functions.combinatorial.numbers import divisor_sigma
>>> divisor_sigma(18, 0)
6
>>> divisor_sigma(39, 1)
56
>>> divisor_sigma(12, 2)
210
>>> divisor_sigma(37)
38

See also

sympy.ntheory.factor_.divisor_count, totient, sympy.ntheory.factor_.divisors, sympy.ntheory.factor_.factorint

References

class sympy.functions.combinatorial.numbers.udivisor_sigma(n, k=1)[source]

Calculate the unitary divisor function \(\sigma_k^*(n)\) for positive integer n

udivisor_sigma(n, k) is equal to sum([x**k for x in udivisors(n)])

If n’s prime factorization is:

\[n = \prod_{i=1}^\omega p_i^{m_i},\]

then

\[\sigma_k^*(n) = \prod_{i=1}^\omega (1+ p_i^{m_ik}).\]
Parameters:

k : power of divisors in the sum

for k = 0, 1: udivisor_sigma(n, 0) is equal to udivisor_count(n) udivisor_sigma(n, 1) is equal to sum(udivisors(n))

Default for k is 1.

Examples

>>> from sympy.functions.combinatorial.numbers import udivisor_sigma
>>> udivisor_sigma(18, 0)
4
>>> udivisor_sigma(74, 1)
114
>>> udivisor_sigma(36, 3)
47450
>>> udivisor_sigma(111)
152

See also

sympy.ntheory.factor_.divisor_count, totient, sympy.ntheory.factor_.divisors, sympy.ntheory.factor_.udivisors, sympy.ntheory.factor_.udivisor_count, divisor_sigma, sympy.ntheory.factor_.factorint

References

class sympy.functions.combinatorial.numbers.legendre_symbol(a, p)[source]

Returns the Legendre symbol \((a / p)\).

For an integer a and an odd prime p, the Legendre symbol is defined as

\[\begin{split}\genfrac(){}{}{a}{p} = \begin{cases} 0 & \text{if } p \text{ divides } a\\ 1 & \text{if } a \text{ is a quadratic residue modulo } p\\ -1 & \text{if } a \text{ is a quadratic nonresidue modulo } p \end{cases}\end{split}\]

Examples

>>> from sympy.functions.combinatorial.numbers import legendre_symbol
>>> [legendre_symbol(i, 7) for i in range(7)]
[0, 1, 1, -1, 1, -1, -1]
>>> sorted(set([i**2 % 7 for i in range(7)]))
[0, 1, 2, 4]

See also

sympy.ntheory.residue_ntheory.is_quad_residue, jacobi_symbol

class sympy.functions.combinatorial.numbers.jacobi_symbol(m, n)[source]

Returns the Jacobi symbol \((m / n)\).

For any integer m and any positive odd integer n the Jacobi symbol is defined as the product of the Legendre symbols corresponding to the prime factors of n:

\[\genfrac(){}{}{m}{n} = \genfrac(){}{}{m}{p^{1}}^{\alpha_1} \genfrac(){}{}{m}{p^{2}}^{\alpha_2} ... \genfrac(){}{}{m}{p^{k}}^{\alpha_k} \text{ where } n = p_1^{\alpha_1} p_2^{\alpha_2} ... p_k^{\alpha_k}\]

Like the Legendre symbol, if the Jacobi symbol \(\genfrac(){}{}{m}{n} = -1\) then m is a quadratic nonresidue modulo n.

But, unlike the Legendre symbol, if the Jacobi symbol \(\genfrac(){}{}{m}{n} = 1\) then m may or may not be a quadratic residue modulo n.

Examples

>>> from sympy.functions.combinatorial.numbers import jacobi_symbol, legendre_symbol
>>> from sympy import S
>>> jacobi_symbol(45, 77)
-1
>>> jacobi_symbol(60, 121)
1

The relationship between the jacobi_symbol and legendre_symbol can be demonstrated as follows:

>>> L = legendre_symbol
>>> S(45).factors()
{3: 2, 5: 1}
>>> jacobi_symbol(7, 45) == L(7, 3)**2 * L(7, 5)**1
True

See also

sympy.ntheory.residue_ntheory.is_quad_residue, legendre_symbol

class sympy.functions.combinatorial.numbers.kronecker_symbol(a, n)[source]

Returns the Kronecker symbol \((a / n)\).

Examples

>>> from sympy.functions.combinatorial.numbers import kronecker_symbol
>>> kronecker_symbol(45, 77)
-1
>>> kronecker_symbol(13, -120)
1

See also

jacobi_symbol, legendre_symbol

References

class sympy.functions.combinatorial.numbers.mobius(n)[source]

Mobius function maps natural number to {-1, 0, 1}

It is defined as follows:

  1. \(1\) if \(n = 1\).

  2. \(0\) if \(n\) has a squared prime factor.

  3. \((-1)^k\) if \(n\) is a square-free positive integer with \(k\) number of prime factors.

It is an important multiplicative function in number theory and combinatorics. It has applications in mathematical series, algebraic number theory and also physics (Fermion operator has very concrete realization with Mobius Function model).

Examples

>>> from sympy.functions.combinatorial.numbers import mobius
>>> mobius(13*7)
1
>>> mobius(1)
1
>>> mobius(13*7*5)
-1
>>> mobius(13**2)
0

Even in the case of a symbol, if it clearly contains a squared prime factor, it will be zero.

>>> from sympy import Symbol
>>> n = Symbol("n", integer=True, positive=True)
>>> mobius(4*n)
0
>>> mobius(n**2)
0

References

[R255]

Thomas Koshy “Elementary Number Theory with Applications”

class sympy.functions.combinatorial.numbers.primenu(n)[source]

Calculate the number of distinct prime factors for a positive integer n.

If n’s prime factorization is:

\[n = \prod_{i=1}^k p_i^{m_i},\]

then primenu(n) or \(\nu(n)\) is:

\[\nu(n) = k.\]

Examples

>>> from sympy.functions.combinatorial.numbers import primenu
>>> primenu(1)
0
>>> primenu(30)
3

See also

sympy.ntheory.factor_.factorint

References

class sympy.functions.combinatorial.numbers.primeomega(n)[source]

Calculate the number of prime factors counting multiplicities for a positive integer n.

If n’s prime factorization is:

\[n = \prod_{i=1}^k p_i^{m_i},\]

then primeomega(n) or \(\Omega(n)\) is:

\[\Omega(n) = \sum_{i=1}^k m_i.\]

Examples

>>> from sympy.functions.combinatorial.numbers import primeomega
>>> primeomega(1)
0
>>> primeomega(20)
3

See also

sympy.ntheory.factor_.factorint

References

class sympy.functions.combinatorial.numbers.totient(n)[source]

Calculate the Euler totient function phi(n)

totient(n) or \(\phi(n)\) is the number of positive integers \(\leq\) n that are relatively prime to n.

Examples

>>> from sympy.functions.combinatorial.numbers import totient
>>> totient(1)
1
>>> totient(25)
20
>>> totient(45) == totient(5)*totient(9)
True

See also

sympy.ntheory.factor_.divisor_count

References

class sympy.functions.combinatorial.numbers.reduced_totient(n)[source]

Calculate the Carmichael reduced totient function lambda(n)

reduced_totient(n) or \(\lambda(n)\) is the smallest m > 0 such that \(k^m \equiv 1 \mod n\) for all k relatively prime to n.

Examples

>>> from sympy.functions.combinatorial.numbers import reduced_totient
>>> reduced_totient(1)
1
>>> reduced_totient(8)
2
>>> reduced_totient(30)
4

See also

totient

References

class sympy.functions.combinatorial.numbers.primepi(n)[source]

Represents the prime counting function pi(n) = the number of prime numbers less than or equal to n.

Examples

>>> from sympy.functions.combinatorial.numbers import primepi
>>> from sympy import prime, prevprime, isprime
>>> primepi(25)
9

So there are 9 primes less than or equal to 25. Is 25 prime?

>>> isprime(25)
False

It is not. So the first prime less than 25 must be the 9th prime:

>>> prevprime(25) == prime(9)
True

See also

sympy.ntheory.primetest.isprime

Test if n is prime

sympy.ntheory.generate.primerange

Generate all primes in a given range

sympy.ntheory.generate.prime

Return the nth prime

References

class sympy.functions.combinatorial.factorials.MultiFactorial(*args)[source]
class sympy.functions.combinatorial.factorials.RisingFactorial(x, k)[source]

Rising factorial (also called Pochhammer symbol [R268]) 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 https://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 [R269].

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)

See also

factorial, factorial2, FallingFactorial

References

[R269] (1,2)

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

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

See also

sympy.utilities.iterables.multiset_partitions

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

See also

sympy.utilities.iterables.multiset_combinations

References

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

See also

sympy.utilities.iterables.multiset_permutations

References

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

See also

sympy.utilities.iterables.partitions, sympy.utilities.iterables.multiset_partitions, sympy.functions.combinatorial.numbers.partition

References