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
\[\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.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
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
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
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 >>> 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
betainc()¶
- 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.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¶
superfac()¶
- 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.dps = 15; 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
hyperfac()¶
- 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.dps = 15; 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
barnesg()¶
- 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.dps = 15; 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
- Whittaker & Watson, A Course of Modern Analysis, Cambridge University Press, 4th edition (1927), p.264
- http://en.wikipedia.org/wiki/Barnes_G-function
- http://mathworld.wolfram.com/BarnesG-Function.html
Polygamma functions and harmonic numbers¶
psi()/digamma()¶
- 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.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([1,-2,3])) [(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(z)¶
Shortcut for psi(0,z).
harmonic()¶
- mpmath.harmonic(z)¶
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.dps = 15; 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
