Navigation

Factorials and gamma functions

Factorials and factorial-like sums and products are basic tools of combinatorics and number theory. Much like the exponential function is fundamental to differential equations and analysis in general, the factorial function (and its extension to complex numbers, the gamma function) is fundamental to difference equations and functional equations.

A large selection of factorial-like functions is implemented in mpmath. All functions support complex arguments, and arguments may be arbitrarily large. Results are numerical approximations, so to compute exact values a high enough precision must be set manually:

>>> mp.dps = 15; mp.pretty = True
>>> fac(100)
9.33262154439442e+157
>>> print int(_)    # most digits are wrong
93326215443944150965646704795953882578400970373184098831012889540582227238570431
295066113089288327277825849664006524270554535976289719382852181865895959724032
>>> mp.dps = 160
>>> fac(100)
93326215443944152681699238856266700490715968264381621468592963895217599993229915
608941463976156518286253697920827223758251185210916864000000000000000000000000.0

The gamma and polygamma functions are closely related to Zeta functions, L-series and polylogarithms. See also q-functions for q-analogs of factorial-like functions.

Factorials

factorial()/fac()

mpmath.factorial(x, **kwargs)

Computes the factorial, x!. For integers n \ge 0, we have n! = 1 \cdot 2 \cdots (n-1) \cdot n and more generally the factorial is defined for real or complex x by x! = \Gamma(x+1).

Examples

Basic values and limits:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> for k in range(6):
...     print("%s %s" % (k, fac(k)))
...
0 1.0
1 1.0
2 2.0
3 6.0
4 24.0
5 120.0
>>> fac(inf)
+inf
>>> fac(0.5), sqrt(pi)/2
(0.886226925452758, 0.886226925452758)

For large positive x, x! can be approximated by Stirling’s formula:

>>> x = 10**10
>>> fac(x)
2.32579620567308e+95657055186
>>> sqrt(2*pi*x)*(x/e)**x
2.32579597597705e+95657055186

fac() supports evaluation for astronomically large values:

>>> fac(10**30)
6.22311232304258e+29565705518096748172348871081098

Reciprocal factorials appear in the Taylor series of the exponential function (among many other contexts):

>>> nsum(lambda k: 1/fac(k), [0, inf]), exp(1)
(2.71828182845905, 2.71828182845905)
>>> nsum(lambda k: pi**k/fac(k), [0, inf]), exp(pi)
(23.1406926327793, 23.1406926327793)

fac2()

mpmath.fac2(x)

Computes the double factorial x!!, defined for integers x > 0 by

x!! = \begin{cases} 1 \cdot 3 \cdots (x-2) \cdot x & x \;\mathrm{odd} \\ 2 \cdot 4 \cdots (x-2) \cdot x & x \;\mathrm{even} \end{cases}

and more generally by [1]

x!! = 2^{x/2} \left(\frac{\pi}{2}\right)^{(\cos(\pi x)-1)/4} \Gamma\left(\frac{x}{2}+1\right).

Examples

The integer sequence of double factorials begins:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> nprint([fac2(n) for n in range(10)])
[1.0, 1.0, 2.0, 3.0, 8.0, 15.0, 48.0, 105.0, 384.0, 945.0]

For large x, double factorials follow a Stirling-like asymptotic approximation:

>>> x = mpf(10000)
>>> fac2(x)
5.97272691416282e+17830
>>> sqrt(pi)*x**((x+1)/2)*exp(-x/2)
5.97262736954392e+17830

The recurrence formula x!! = x (x-2)!! can be reversed to define the double factorial of negative odd integers (but not negative even integers):

>>> fac2(-1), fac2(-3), fac2(-5), fac2(-7)
(1.0, -1.0, 0.333333333333333, -0.0666666666666667)
>>> fac2(-2)
Traceback (most recent call last):
  ...
ValueError: gamma function pole
Traceback (most recent call last):
  ...
ValueError: gamma function pole

With the exception of the poles at negative even integers, fac2() supports evaluation for arbitrary complex arguments. The recurrence formula is valid generally:

>>> fac2(pi+2j)
(-1.3697207890154e-12 + 3.93665300979176e-12j)
>>> (pi+2j)*fac2(pi-2+2j)
(-1.3697207890154e-12 + 3.93665300979176e-12j)

Double factorials should not be confused with nested factorials, which are immensely larger:

>>> fac(fac(20))
5.13805976125208e+43675043585825292774
>>> fac2(20)
3715891200.0

Double factorials appear, among other things, in series expansions of Gaussian functions and the error function. Infinite series include:

>>> nsum(lambda k: 1/fac2(k), [0, inf])
3.05940740534258
>>> sqrt(e)*(1+sqrt(pi/2)*erf(sqrt(2)/2))
3.05940740534258
>>> nsum(lambda k: 2**k/fac2(2*k-1), [1, inf])
4.06015693855741
>>> e * erf(1) * sqrt(pi)
4.06015693855741

A beautiful Ramanujan sum:

>>> nsum(lambda k: (-1)**k*(fac2(2*k-1)/fac2(2*k))**3, [0,inf])
0.90917279454693
>>> (gamma('9/8')/gamma('5/4')/gamma('7/8'))**2
0.90917279454693

References

  1. http://functions.wolfram.com/GammaBetaErf/Factorial2/27/01/0002/
  2. http://mathworld.wolfram.com/DoubleFactorial.html

Binomial coefficients

binomial()

mpmath.binomial(n, k)

Computes the binomial coefficient

{n \choose k} = \frac{n!}{k!(n-k)!}.

The binomial coefficient gives the number of ways that k items can be chosen from a set of n items. More generally, the binomial coefficient is a well-defined function of arbitrary real or complex n and k, via the gamma function.

Examples

Generate Pascal’s triangle:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> for n in range(5):
...     nprint([binomial(n,k) for k in range(n+1)])
...
[1.0]
[1.0, 1.0]
[1.0, 2.0, 1.0]
[1.0, 3.0, 3.0, 1.0]
[1.0, 4.0, 6.0, 4.0, 1.0]

There is 1 way to select 0 items from the empty set, and 0 ways to select 1 item from the empty set:

>>> binomial(0, 0)
1.0
>>> binomial(0, 1)
0.0

binomial() supports large arguments:

>>> binomial(10**20, 10**20-5)
8.33333333333333e+97
>>> binomial(10**20, 10**10)
2.60784095465201e+104342944813

Nonintegral binomial coefficients find use in series expansions:

>>> nprint(taylor(lambda x: (1+x)**0.25, 0, 4))
[1.0, 0.25, -0.09375, 0.0546875, -0.0375977]
>>> nprint([binomial(0.25, k) for k in range(5)])
[1.0, 0.25, -0.09375, 0.0546875, -0.0375977]

An integral representation:

>>> n, k = 5, 3
>>> f = lambda t: exp(-j*k*t)*(1+exp(j*t))**n
>>> chop(quad(f, [-pi,pi])/(2*pi))
10.0
>>> binomial(n,k)
10.0

Gamma function

gamma()

mpmath.gamma(x, **kwargs)

Computes the gamma function, \Gamma(x). The gamma function is a shifted version of the ordinary factorial, satisfying \Gamma(n) = (n-1)! for integers n > 0. More generally, it is defined by

\Gamma(x) = \int_0^{\infty} t^{x-1} e^{-t}\, dt

for any real or complex x with \Re(x) > 0 and for \Re(x) < 0 by analytic continuation.

Examples

Basic values and limits:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> for k in range(1, 6):
...     print("%s %s" % (k, gamma(k)))
...
1 1.0
2 1.0
3 2.0
4 6.0
5 24.0
>>> gamma(inf)
+inf
>>> gamma(0)
Traceback (most recent call last):
  ...
ValueError: gamma function pole
Traceback (most recent call last):
  ...
ValueError: gamma function pole

The gamma function of a half-integer is a rational multiple of \sqrt{\pi}:

>>> gamma(0.5), sqrt(pi)
(1.77245385090552, 1.77245385090552)
>>> gamma(1.5), sqrt(pi)/2
(0.886226925452758, 0.886226925452758)

We can check the integral definition:

>>> gamma(3.5)
3.32335097044784
>>> quad(lambda t: t**2.5*exp(-t), [0,inf])
3.32335097044784

gamma() supports arbitrary-precision evaluation and complex arguments:

>>> mp.dps = 50
>>> gamma(sqrt(3))
0.91510229697308632046045539308226554038315280564184
>>> mp.dps = 25
>>> gamma(2j)
(0.009902440080927490985955066 - 0.07595200133501806872408048j)

Arguments can also be large. Note that the gamma function grows very quickly:

>>> mp.dps = 15
>>> gamma(10**20)
1.9328495143101e+1956570551809674817225

rgamma()

mpmath.rgamma(x, **kwargs)

Computes the reciprocal of the gamma function, 1/\Gamma(z). This function evaluates to zero at the poles of the gamma function, z = 0, -1, -2, \ldots.

Examples

Basic examples:

>>> from mpmath import *
>>> mp.dps = 25; mp.pretty = True
>>> rgamma(1)
1.0
>>> rgamma(4)
0.1666666666666666666666667
>>> rgamma(0); rgamma(-1)
0.0
0.0
>>> rgamma(1000)
2.485168143266784862783596e-2565
>>> rgamma(inf)
0.0

A definite integral that can be evaluated in terms of elementary integrals:

>>> quad(rgamma, [0,inf])
2.807770242028519365221501
>>> e + quad(lambda t: exp(-t)/(pi**2+log(t)**2), [0,inf])
2.807770242028519365221501

gammaprod()

mpmath.gammaprod(a, b)

Given iterables a and b, gammaprod(a, b) computes the product / quotient of gamma functions:

\frac{\Gamma(a_0) \Gamma(a_1) \cdots \Gamma(a_p)} {\Gamma(b_0) \Gamma(b_1) \cdots \Gamma(b_q)}

Unlike direct calls to gamma(), gammaprod() considers the entire product as a limit and evaluates this limit properly if any of the numerator or denominator arguments are nonpositive integers such that poles of the gamma function are encountered. That is, gammaprod() evaluates

\lim_{\epsilon \to 0} \frac{\Gamma(a_0+\epsilon) \Gamma(a_1+\epsilon) \cdots \Gamma(a_p+\epsilon)} {\Gamma(b_0+\epsilon) \Gamma(b_1+\epsilon) \cdots \Gamma(b_q+\epsilon)}

In particular:

  • If there are equally many poles in the numerator and the denominator, the limit is a rational number times the remaining, regular part of the product.
  • If there are more poles in the numerator, gammaprod() returns +inf.
  • If there are more poles in the denominator, gammaprod() returns 0.

Examples

The reciprocal gamma function 1/\Gamma(x) evaluated at x = 0:

>>> from mpmath import *
>>> mp.dps = 15
>>> gammaprod([], [0])
0.0

A limit:

>>> gammaprod([-4], [-3])
-0.25
>>> limit(lambda x: gamma(x-1)/gamma(x), -3, direction=1)
-0.25
>>> limit(lambda x: gamma(x-1)/gamma(x), -3, direction=-1)
-0.25

loggamma()

mpmath.loggamma(x)

Computes the principal branch of the log-gamma function, \ln \Gamma(z). Unlike \ln(\Gamma(z)), which has infinitely many complex branch cuts, the principal log-gamma function only has a single branch cut along the negative half-axis. The principal branch continuously matches the asymptotic Stirling expansion

\ln \Gamma(z) \sim \frac{\ln(2 \pi)}{2} + \left(z-\frac{1}{2}\right) \ln(z) - z + O(z^{-1}).

The real parts of both functions agree, but their imaginary parts generally differ by 2 n \pi for some n \in \mathbb{Z}. They coincide for z \in \mathbb{R}, z > 0.

Computationally, it is advantageous to use loggamma() instead of gamma() for extremely large arguments.

Examples

Comparing with \ln(\Gamma(z)):

>>> from mpmath import *
>>> mp.dps = 25; mp.pretty = True
>>> loggamma('13.2'); log(gamma('13.2'))
20.49400419456603678498394
20.49400419456603678498394
>>> loggamma(3+4j)
(-1.756626784603784110530604 + 4.742664438034657928194889j)
>>> log(gamma(3+4j))
(-1.756626784603784110530604 - 1.540520869144928548730397j)
>>> log(gamma(3+4j)) + 2*pi*j
(-1.756626784603784110530604 + 4.742664438034657928194889j)

Note the imaginary parts for negative arguments:

>>> loggamma(-0.5); loggamma(-1.5); loggamma(-2.5)
(1.265512123484645396488946 - 3.141592653589793238462643j)
(0.8600470153764810145109327 - 6.283185307179586476925287j)
(-0.05624371649767405067259453 - 9.42477796076937971538793j)

Some special values:

>>> loggamma(1); loggamma(2)
0.0
0.0
>>> loggamma(3); +ln2
0.6931471805599453094172321
0.6931471805599453094172321
>>> loggamma(3.5); log(15*sqrt(pi)/8)
1.200973602347074224816022
1.200973602347074224816022
>>> loggamma(inf)
+inf

Huge arguments are permitted:

>>> loggamma('1e30')
6.807755278982137052053974e+31
>>> loggamma('1e300')
6.897755278982137052053974e+302
>>> loggamma('1e3000')
6.906755278982137052053974e+3003
>>> loggamma('1e100000000000000000000')
2.302585092994045684007991e+100000000000000000020
>>> loggamma('1e30j')
(-1.570796326794896619231322e+30 + 6.807755278982137052053974e+31j)
>>> loggamma('1e300j')
(-1.570796326794896619231322e+300 + 6.897755278982137052053974e+302j)
>>> loggamma('1e3000j')
(-1.570796326794896619231322e+3000 + 6.906755278982137052053974e+3003j)

The log-gamma function can be integrated analytically on any interval of unit length:

>>> z = 0
>>> quad(loggamma, [z,z+1]); log(2*pi)/2
0.9189385332046727417803297
0.9189385332046727417803297
>>> z = 3+4j
>>> quad(loggamma, [z,z+1]); (log(z)-1)*z + log(2*pi)/2
(-0.9619286014994750641314421 + 5.219637303741238195688575j)
(-0.9619286014994750641314421 + 5.219637303741238195688575j)

The derivatives of the log-gamma function are given by the polygamma function (psi()):

>>> diff(loggamma, -4+3j); psi(0, -4+3j)
(1.688493531222971393607153 + 2.554898911356806978892748j)
(1.688493531222971393607153 + 2.554898911356806978892748j)
>>> diff(loggamma, -4+3j, 2); psi(1, -4+3j)
(-0.1539414829219882371561038 - 0.1020485197430267719746479j)
(-0.1539414829219882371561038 - 0.1020485197430267719746479j)

The log-gamma function satisfies an additive form of the recurrence relation for the ordinary gamma function:

>>> z = 2+3j
>>> loggamma(z); loggamma(z+1) - log(z)
(-2.092851753092733349564189 + 2.302396543466867626153708j)
(-2.092851753092733349564189 + 2.302396543466867626153708j)

Rising and falling factorials

rf()

mpmath.rf(x, n)

Computes the rising factorial or Pochhammer symbol,

x^{(n)} = x (x+1) \cdots (x+n-1) = \frac{\Gamma(x+n)}{\Gamma(x)}

where the rightmost expression is valid for nonintegral n.

Examples

For integral n, the rising factorial is a polynomial:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> for n in range(5):
...     nprint(taylor(lambda x: rf(x,n), 0, n))
...
[1.0]
[0.0, 1.0]
[0.0, 1.0, 1.0]
[0.0, 2.0, 3.0, 1.0]
[0.0, 6.0, 11.0, 6.0, 1.0]

Evaluation is supported for arbitrary arguments:

>>> rf(2+3j, 5.5)
(-7202.03920483347 - 3777.58810701527j)

ff()

mpmath.ff(x, n)

Computes the falling factorial,

(x)_n = x (x-1) \cdots (x-n+1) = \frac{\Gamma(x+1)}{\Gamma(x-n+1)}

where the rightmost expression is valid for nonintegral n.

Examples

For integral n, the falling factorial is a polynomial:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> for n in range(5):
...     nprint(taylor(lambda x: ff(x,n), 0, n))
...
[1.0]
[0.0, 1.0]
[0.0, -1.0, 1.0]
[0.0, 2.0, -3.0, 1.0]
[0.0, -6.0, 11.0, -6.0, 1.0]

Evaluation is supported for arbitrary arguments:

>>> ff(2+3j, 5.5)
(-720.41085888203 + 316.101124983878j)

Beta function

beta()

mpmath.beta(x, y)

Computes the beta function, B(x,y) = \Gamma(x) \Gamma(y) / \Gamma(x+y). The beta function is also commonly defined by the integral representation

B(x,y) = \int_0^1 t^{x-1} (1-t)^{y-1} \, dt

Examples

For integer and half-integer arguments where all three gamma functions are finite, the beta function becomes either rational number or a rational multiple of \pi:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> beta(5, 2)
0.0333333333333333
>>> beta(1.5, 2)
0.266666666666667
>>> 16*beta(2.5, 1.5)
3.14159265358979

Where appropriate, beta() evaluates limits. A pole of the beta function is taken to result in +inf:

>>> beta(-0.5, 0.5)
0.0
>>> beta(-3, 3)
-0.333333333333333
>>> beta(-2, 3)
+inf
>>> beta(inf, 1)
0.0
>>> beta(inf, 0)
nan

beta() supports complex numbers and arbitrary precision evaluation:

>>> beta(1, 2+j)
(0.4 - 0.2j)
>>> mp.dps = 25
>>> beta(j,0.5)
(1.079424249270925780135675 - 1.410032405664160838288752j)
>>> mp.dps = 50
>>> beta(pi, e)
0.037890298781212201348153837138927165984170287886464

Various integrals can be computed by means of the beta function:

>>> mp.dps = 15
>&