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

    >>> from mpmath import mp, fac
    >>> mp.dps = 15
    >>> mp.pretty = True
    >>> fac(100)
    9.33262154439442e+157
    >>> print(int(_))    # most digits are wrong
    93326215443944150965646704795953882578400970373184098831012889540582227238570431295066113089288327277825849664006524270554535976289719382852181865895959724032
    >>> mp.dps = 160
    >>> fac(100)
    93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000.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

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, fac, sqrt, inf, pi, exp, nsum
      >>> 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

  `factorial()` 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)

mpmath.fac2(x)

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

  \[\begin{split}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}\end{split}\]

  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, fac2, nprint, mpf, exp, pi, nsum, sqrt, e,
      ...                     gamma, inf, erf, fac)
      >>> 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

  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. [WolframFunctions] http://functions.wolfram.com/GammaBetaErf/Factorial2/27/01/0002/

  2. [Weisstein] http://mathworld.wolfram.com/DoubleFactorial.html

## Binomial coefficients

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, binomial, nprint, exp, taylor, j, chop, quad, pi
      >>> 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

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, gamma, inf, sqrt, pi, exp, quad
      >>> 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

  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

  References

  - [Spouge]

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, rgamma, inf, pi, log, exp, quad, e
      >>> mp.dps = 25
      >>> mp.pretty = True
      >>> rgamma(1)
      1.0
      >>> rgamma(4)
      0.1666666666666666666666667
      >>> rgamma(0)
      0.0
      >>> rgamma(-1)
      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

mpmath.gammaprod(a, b, _infsign=False)

  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, gammaprod, limit, gamma
      >>> mp.pretty = True
      >>> 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

mpmath.loggamma(x, **kwargs)

  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, log, loggamma, pi, j, ln2, sqrt, inf, quad,
      ...                     diff, psi)
      >>> mp.dps = 25
      >>> mp.pretty = True
      >>> loggamma('13.2')
      20.49400419456603678498394
      >>> log(gamma('13.2'))
      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)
      (1.265512123484645396488946 - 3.141592653589793238462643j)
      >>> loggamma(-1.5)
      (0.8600470153764810145109327 - 6.283185307179586476925287j)
      >>> loggamma(-2.5)
      (-0.05624371649767405067259453 - 9.42477796076937971538793j)

  Some special values:

      >>> loggamma(1)
      0.0
      >>> loggamma(2)
      0.0
      >>> loggamma(3)
      0.6931471805599453094172321
      >>> +ln2
      0.6931471805599453094172321
      >>> loggamma(3.5)
      1.200973602347074224816022
      >>> log(15*sqrt(pi)/8)
      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])
      0.9189385332046727417803297
      >>> log(2*pi)/2
      0.9189385332046727417803297
      >>> z = 3+4j
      >>> quad(loggamma, [z,z+1])
      (-0.9619286014994750641314421 + 5.219637303741238195688575j)
      >>> (log(z)-1)*z + log(2*pi)/2
      (-0.9619286014994750641314421 + 5.219637303741238195688575j)

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

      >>> diff(loggamma, -4+3j)
      (1.688493531222971393607153 + 2.554898911356806978892748j)
      >>> psi(0, -4+3j)
      (1.688493531222971393607153 + 2.554898911356806978892748j)
      >>> diff(loggamma, -4+3j, 2)
      (-0.1539414829219882371561038 - 0.1020485197430267719746479j)
      >>> psi(1, -4+3j)
      (-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)
      (-2.092851753092733349564189 + 2.302396543466867626153708j)
      >>> loggamma(z+1) - log(z)
      (-2.092851753092733349564189 + 2.302396543466867626153708j)

## Rising and falling factorials

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 rf, mp, nprint, taylor
      >>> 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)

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, ff, nprint, taylor
      >>> 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

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, beta, inf, j, pi, e, quad, sqrt, sin, cos
      >>> 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
      >>> quad(lambda t: t**2.5*(1-t)**2, [0, 1])
      0.0230880230880231
      >>> beta(3.5, 3)
      0.0230880230880231
      >>> quad(lambda t: sin(t)**4 * sqrt(cos(t)), [0, pi/2])
      0.319504062596158
      >>> beta(2.5, 0.75)/2
      0.319504062596158

mpmath.betainc(a, b, x1=0, x2=1, regularized=False)

  `betainc(a, b, x1=0, x2=1, regularized=False)` gives the generalized incomplete beta function,

  \[I_{x_1}^{x_2}(a,b) = \int_{x_1}^{x_2} t^{a-1} (1-t)^{b-1} dt.\]

  When \(x_1 = 0, x_2 = 1\), this reduces to the ordinary (complete) beta function \(B(a,b)\); see `beta()`.

  With the keyword argument `regularized=True`, `betainc()` computes the regularized incomplete beta function \(I_{x_1}^{x_2}(a,b) / B(a,b)\). This is the cumulative distribution of the beta distribution with parameters \(a\), \(b\).

  Examples

  Verifying that `betainc()` computes the integral in the definition:

      >>> from mpmath import mp, betainc, quad, identify, chop, pi, e
      >>> mp.dps = 25
      >>> mp.pretty = True
      >>> x,y,a,b = 3, 4, 0, 6
      >>> betainc(x, y, a, b)
      -4010.4
      >>> quad(lambda t: t**(x-1) * (1-t)**(y-1), [a, b])
      -4010.4

  The arguments may be arbitrary complex numbers:

      >>> betainc(0.75, 1-4j, 0, 2+3j)
      (0.2241657956955709603655887 + 0.3619619242700451992411724j)

  With regularization:

      >>> betainc(1, 2, 0, 0.25, regularized=True)
      0.4375
      >>> betainc(pi, e, 0, 1, regularized=True)   # Complete
      1.0

  The beta integral satisfies some simple argument transformation symmetries:

      >>> mp.dps = 15
      >>> betainc(2,3,4,5), -betainc(2,3,5,4), betainc(3,2,1-5,1-4)
      (56.0833333333333, 56.0833333333333, 56.0833333333333)

  The beta integral can often be evaluated analytically. For integer and rational arguments, the incomplete beta function typically reduces to a simple algebraic-logarithmic expression:

      >>> mp.dps = 25
      >>> identify(chop(betainc(0, 0, 3, 4)))
      '-(log((9/8)))'
      >>> identify(betainc(2, 3, 4, 5))
      '(673/12)'
      >>> identify(betainc(1.5, 1, 1, 2))
      '((-12+sqrt(1152))/18)'

## Super- and hyperfactorials

mpmath.superfac(z)

  Computes the superfactorial, defined as the product of consecutive factorials

  \[\mathrm{sf}(n) = \prod_{k=1}^n k!\]

  For general complex \(z\), \(\mathrm{sf}(z)\) is defined in terms of the Barnes G-function (see `barnesg()`).

  Examples

  The first few superfactorials are (OEIS A000178):

      >>> from mpmath import mp, superfac, pi, diff
      >>> mp.pretty = True
      >>> for n in range(10):
      ...     print("%s %s" % (n, superfac(n)))
      ...
      0 1.0
      1 1.0
      2 2.0
      3 12.0
      4 288.0
      5 34560.0
      6 24883200.0
      7 125411328000.0
      8 5.05658474496e+15
      9 1.83493347225108e+21

  Superfactorials grow very rapidly:

      >>> superfac(1000)
      3.24570818422368e+1177245
      >>> superfac(10**10)
      2.61398543581249e+467427913956904067453

  Evaluation is supported for arbitrary arguments:

      >>> mp.dps = 25
      >>> superfac(pi)
      17.20051550121297985285333
      >>> superfac(2+3j)
      (-0.005915485633199789627466468 + 0.008156449464604044948738263j)
      >>> diff(superfac, 1)
      0.2645072034016070205673056

  References

  1. [OEIS] http://oeis.org/A000178

mpmath.hyperfac(z)

  Computes the hyperfactorial, defined for integers as the product

  \[H(n) = \prod_{k=1}^n k^k.\]

  The hyperfactorial satisfies the recurrence formula \(H(z) = z^z H(z-1)\). It can be defined more generally in terms of the Barnes G-function (see `barnesg()`) and the gamma function by the formula

  \[H(z) = \frac{\Gamma(z+1)^z}{G(z)}.\]

  The extension to complex numbers can also be done via the integral representation

  \[H(z) = (2\pi)^{-z/2} \exp \left[ {z+1 \choose 2} + \int_0^z \log(t!)\,dt \right].\]

  Examples

  The rapidly-growing sequence of hyperfactorials begins (OEIS A002109):

      >>> from mpmath import (mp, hyperfac, diff, pi, mpf, chop, exp, quad,
      ...                     binomial, loggamma, sqrt, j)
      >>> mp.pretty = True
      >>> for n in range(10):
      ...     print("%s %s" % (n, hyperfac(n)))
      ...
      0 1.0
      1 1.0
      2 4.0
      3 108.0
      4 27648.0
      5 86400000.0
      6 4031078400000.0
      7 3.3197663987712e+18
      8 5.56964379417266e+25
      9 2.15779412229419e+34

  Some even larger hyperfactorials are:

      >>> hyperfac(1000)
      5.46458120882585e+1392926
      >>> hyperfac(10**10)
      4.60408207642219e+489142638002418704309

  The hyperfactorial can be evaluated for arbitrary arguments:

      >>> hyperfac(0.5)
      0.880449235173423
      >>> diff(hyperfac, 1)
      0.581061466795327
      >>> hyperfac(pi)
      205.211134637462
      >>> hyperfac(-10+1j)
      (3.01144471378225e+46 - 2.45285242480185e+46j)

  The recurrence property of the hyperfactorial holds generally:

      >>> z = 3-4*j
      >>> hyperfac(z)
      (-4.49795891462086e-7 - 6.33262283196162e-7j)
      >>> z**z * hyperfac(z-1)
      (-4.49795891462086e-7 - 6.33262283196162e-7j)
      >>> z = mpf(-0.6)
      >>> chop(z**z * hyperfac(z-1))
      1.28170142849352
      >>> hyperfac(z)
      1.28170142849352

  The hyperfactorial may also be computed using the integral definition:

      >>> z = 2.5
      >>> hyperfac(z)
      15.9842119922237
      >>> (2*pi)**(-z/2)*exp(binomial(z+1,2) +
      ...     quad(lambda t: loggamma(t+1), [0, z]))
      15.9842119922237

  `hyperfac()` supports arbitrary-precision evaluation:

      >>> mp.dps = 50
      >>> hyperfac(10)
      215779412229418562091680268288000000000000000.0
      >>> hyperfac(1/sqrt(2))
      0.89404818005227001975423476035729076375705084390942

  References

  1. [OEIS] http://oeis.org/A002109

  2. [Weisstein] http://mathworld.wolfram.com/Hyperfactorial.html

mpmath.barnesg(z)

  Evaluates the Barnes G-function, which generalizes the superfactorial (`superfac()`) and by extension also the hyperfactorial (`hyperfac()`) to the complex numbers in an analogous way to how the gamma function generalizes the ordinary factorial.

  The Barnes G-function may be defined in terms of a Weierstrass product:

  \[G(z+1) = (2\pi)^{z/2} e^{-[z(z+1)+\gamma z^2]/2} \prod_{n=1}^\infty \left[\left(1+\frac{z}{n}\right)^ne^{-z+z^2/(2n)}\right]\]

  For positive integers \(n\), we have have relation to superfactorials \(G(n) = \mathrm{sf}(n-2) = 0! \cdot 1! \cdots (n-2)!\).

  Examples

  Some elementary values and limits of the Barnes G-function:

      >>> from mpmath import (mp, barnesg, sqrt, exp, log, glaisher, inf,
      ...                     catalan, pi, nthroot, gamma, limit, mpf, psi, diff)
      >>> mp.pretty = True
      >>> barnesg(1), barnesg(2), barnesg(3)
      (1.0, 1.0, 1.0)
      >>> barnesg(4)
      2.0
      >>> barnesg(5)
      12.0
      >>> barnesg(6)
      288.0
      >>> barnesg(7)
      34560.0
      >>> barnesg(8)
      24883200.0
      >>> barnesg(inf)
      inf
      >>> barnesg(0), barnesg(-1), barnesg(-2)
      (0.0, 0.0, 0.0)

  Closed-form values are known for some rational arguments:

      >>> barnesg('1/2')
      0.603244281209446
      >>> sqrt(exp(0.25+log(2)/12)/sqrt(pi)/glaisher**3)
      0.603244281209446
      >>> barnesg('1/4')
      0.29375596533861
      >>> nthroot(exp('3/8')/exp(catalan/pi)/
      ...      gamma(0.25)**3/sqrt(glaisher)**9, 4)
      0.29375596533861

  The Barnes G-function satisfies the functional equation \(G(z+1) = \Gamma(z) G(z)\):

      >>> z = pi
      >>> barnesg(z+1)
      2.39292119327948
      >>> gamma(z)*barnesg(z)
      2.39292119327948

  The asymptotic growth rate of the Barnes G-function is related to the Glaisher-Kinkelin constant:

      >>> limit(lambda n: barnesg(n+1)/(n**(n**2/2-mpf(1)/12)*
      ...     (2*pi)**(n/2)*exp(-3*n**2/4)), inf)
      0.847536694177301
      >>> exp('1/12')/glaisher
      0.847536694177301

  The Barnes G-function can be differentiated in closed form:

      >>> z = 3
      >>> diff(barnesg, z)
      0.264507203401607
      >>> barnesg(z)*((z-1)*psi(0,z)-z+(log(2*pi)+1)/2)
      0.264507203401607

  Evaluation is supported for arbitrary arguments and at arbitrary precision:

      >>> barnesg(6.5)
      2548.7457695685
      >>> barnesg(-pi)
      0.00535976768353037
      >>> barnesg(3+4j)
      (-0.000676375932234244 - 4.42236140124728e-5j)
      >>> mp.dps = 50
      >>> barnesg(1/sqrt(2))
      0.81305501090451340843586085064413533788206204124732
      >>> q = barnesg(10j)
      >>> q.real
      0.000000000021852360840356557241543036724799812371995850552234
      >>> q.imag
      -0.00000000000070035335320062304849020654215545839053210041457588
      >>> mp.dps = 15
      >>> barnesg(100)
      3.10361006263698e+6626
      >>> barnesg(-101)
      0.0
      >>> barnesg(-10.5)
      5.94463017605008e+25
      >>> barnesg(-10000.5)
      -6.14322868174828e+167480422
      >>> barnesg(1000j)
      (5.21133054865546e-1173597 + 4.27461836811016e-1173597j)
      >>> barnesg(-1000+1000j)
      (2.43114569750291e+1026623 + 2.24851410674842e+1026623j)

  References

  1. [WhittakerWatson], p.264

  2. [Wikipedia] http://en.wikipedia.org/wiki/Barnes_G-function

  3. [Weisstein] http://mathworld.wolfram.com/BarnesG-Function.html

## Polygamma functions and harmonic numbers

mpmath.psi(m, z)

  Gives the polygamma function of order \(m\) of \(z\), \(\psi^{(m)}(z)\). Special cases are known as the digamma function (\(\psi^{(0)}(z)\)), the trigamma function (\(\psi^{(1)}(z)\)), etc. The polygamma functions are defined as the logarithmic derivatives of the gamma function:

  \[\psi^{(m)}(z) = \left(\frac{d}{dz}\right)^{m+1} \log \Gamma(z)\]

  In particular, \(\psi^{(0)}(z) = \Gamma'(z)/\Gamma(z)\). In the present implementation of `psi()`, the order \(m\) must be a nonnegative integer, while the argument \(z\) may be an arbitrary complex number (with exception for the polygamma function’s poles at \(z = 0, -1, -2, \ldots\)).

  Examples

  For various rational arguments, the polygamma function reduces to a combination of standard mathematical constants:

      >>> from mpmath import (mp, psi, euler, catalan, pi, apery, quad, diff,
      ...                     sqrt, nsum, inf, j, nprint, polyroots)
      >>> mp.dps = 25
      >>> mp.pretty = True
      >>> psi(0, 1), -euler
      (-0.5772156649015328606065121, -0.5772156649015328606065121)
      >>> psi(1, '1/4'), pi**2+8*catalan
      (17.19732915450711073927132, 17.19732915450711073927132)
      >>> psi(2, '1/2'), -14*apery
      (-16.82879664423431999559633, -16.82879664423431999559633)

  The polygamma functions are derivatives of each other:

      >>> diff(lambda x: psi(3, x), pi), psi(4, pi)
      (-0.1105749312578862734526952, -0.1105749312578862734526952)
      >>> quad(lambda x: psi(4, x), [2, 3]), psi(3,3)-psi(3,2)
      (-0.375, -0.375)

  The digamma function diverges logarithmically as \(z \to \infty\), while higher orders tend to zero:

      >>> psi(0,inf), psi(1,inf), psi(2,inf)
      (inf, 0.0, 0.0)

  Evaluation for a complex argument:

      >>> psi(2, -1-2j)
      (0.03902435405364952654838445 + 0.1574325240413029954685366j)

  Evaluation is supported for large orders \(m\) and/or large arguments \(z\):

      >>> psi(3, 10**100)
      2.0e-300
      >>> psi(250, 10**30+10**20*j)
      (-1.293142504363642687204865e-7010 + 3.232856260909107391513108e-7018j)

  Application to infinite series

  Any infinite series where the summand is a rational function of the index \(k\) can be evaluated in closed form in terms of polygamma functions of the roots and poles of the summand:

      >>> a = sqrt(2)
      >>> b = sqrt(3)
      >>> nsum(lambda k: 1/((k+a)**2*(k+b)), [0, inf])
      0.4049668927517857061917531
      >>> (psi(0,a)-psi(0,b)-a*psi(1,a)+b*psi(1,a))/(a-b)**2
      0.4049668927517857061917531

  This follows from the series representation (\(m > 0\))

  \[\psi^{(m)}(z) = (-1)^{m+1} m! \sum_{k=0}^{\infty} \frac{1}{(z+k)^{m+1}}.\]

  Since the roots of a polynomial may be complex, it is sometimes necessary to use the complex polygamma function to evaluate an entirely real-valued sum:

      >>> nsum(lambda k: 1/(k**2-2*k+3), [0, inf])
      1.694361433907061256154665
      >>> nprint(polyroots([3,-2,1], asc=True))
      [(1.0 - 1.41421j), (1.0 + 1.41421j)]
      >>> r1 = 1-sqrt(2)*j
      >>> r2 = r1.conjugate()
      >>> (psi(0,-r2)-psi(0,-r1))/(r1-r2)
      (1.694361433907061256154665 + 0.0j)

mpmath.digamma(x, **kwargs)

  Shortcut for `psi(0,z)`.

mpmath.harmonic(x, **kwargs)

  If \(n\) is an integer, `harmonic(n)` gives a floating-point approximation of the \(n\)-th harmonic number \(H(n)\), defined as

  \[H(n) = 1 + \frac{1}{2} + \frac{1}{3} + \ldots + \frac{1}{n}\]

  The first few harmonic numbers are:

      >>> from mpmath import mp, harmonic, inf, pi, findroot, ceil
      >>> mp.pretty = True
      >>> for n in range(8):
      ...     print("%s %s" % (n, harmonic(n)))
      ...
      0 0.0
      1 1.0
      2 1.5
      3 1.83333333333333
      4 2.08333333333333
      5 2.28333333333333
      6 2.45
      7 2.59285714285714

  The infinite harmonic series \(1 + 1/2 + 1/3 + \ldots\) diverges:

      >>> harmonic(inf)
      inf

  `harmonic()` is evaluated using the digamma function rather than by summing the harmonic series term by term. It can therefore be computed quickly for arbitrarily large \(n\), and even for nonintegral arguments:

      >>> harmonic(10**100)
      230.835724964306
      >>> harmonic(0.5)
      0.613705638880109
      >>> harmonic(3+4j)
      (2.24757548223494 + 0.850502209186044j)

  `harmonic()` supports arbitrary precision evaluation:

      >>> mp.dps = 50
      >>> harmonic(11)
      3.0198773448773448773448773448773448773448773448773
      >>> harmonic(pi)
      1.8727388590273302654363491032336134987519132374152

  The harmonic series diverges, but at a glacial pace. It is possible to calculate the exact number of terms required before the sum exceeds a given amount, say 100:

      >>> mp.dps = 50
      >>> v = 10**findroot(lambda x: harmonic(10**x) - 100, 10)
      >>> v
      15092688622113788323693563264538101449859496.864101
      >>> v = int(ceil(v))
      >>> print(v)
      15092688622113788323693563264538101449859497
      >>> harmonic(v-1)
      99.999999999999999999999999999999999999999999942747
      >>> harmonic(v)
      100.000000000000000000000000000000000000000000009