# Digamma function

(Redirected from Gauss's digamma theorem)
The color representation of the Digamma function, ${\displaystyle \psi (z)}$, in a rectangular region of the complex plane

In mathematics, the digamma function is defined as the logarithmic derivative of the gamma function:[1][2]

${\displaystyle \psi (x)={\frac {d}{dx}}\ln {\big (}\Gamma (x){\big )}={\frac {\Gamma '(x)}{\Gamma (x)}}.}$

It is the first of the polygamma functions.

The digamma function is often denoted as ψ0(x), ψ(0)(x) or Ϝ (the uppercase form of the archaic Greek consonant digamma) meaning Double-Gamma.

## Relation to harmonic numbers

The gamma function obeys the equation

${\displaystyle \Gamma (z+1)=z\Gamma (z).}$

Taking the derivative with respect to z gives:

${\displaystyle \Gamma '(z+1)=z\Gamma '(z)+\Gamma (z)}$

Dividing by Γ(z + 1) or the equivalent zΓ(z) gives:

${\displaystyle {\frac {\Gamma '(z+1)}{\Gamma (z+1)}}={\frac {\Gamma '(z)}{\Gamma (z)}}+{\frac {1}{z}}}$

or:

${\displaystyle \psi (z+1)=\psi (z)+{\frac {1}{z}}}$

Since the harmonic numbers are defined as

${\displaystyle H_{n}=\sum _{k=1}^{n}{\frac {1}{k}}}$

the digamma function is related to it by:

${\displaystyle \psi (n)=H_{n-1}-\gamma }$

where Hn is the nth harmonic number, and γ is the Euler–Mascheroni constant. For half-integer values, it may be expressed as

${\displaystyle \psi \left(n+{\tfrac {1}{2}}\right)=-\gamma -2\ln 2+\sum _{k=1}^{n}{\frac {2}{2k-1}}}$

## Integral representations

If the real part of x is positive then the digamma function has the following integral representation

${\displaystyle \psi (x)=\int _{0}^{\infty }\left({\frac {e^{-t}}{t}}-{\frac {e^{-xt}}{1-e^{-t}}}\right)\,dt}$.

This may be written as

${\displaystyle \psi (s+1)=-\gamma +\int _{0}^{1}\left({\frac {1-x^{s}}{1-x}}\right)\,dx}$

which follows from Leonhard Euler's integral formula for the harmonic numbers.

## Infinite product representation

The function ${\displaystyle \psi (z)/\Gamma (z)}$ is an entire function[3], and it can be represented by the infinite product

${\displaystyle {\frac {\psi (z)}{\Gamma (z)}}=-e^{2\gamma z}\prod _{k=0}^{\infty }\left(1-{\frac {z}{x_{k}}}\right)e^{\frac {z}{x_{k}}}.}$

Here ${\displaystyle x_{k}}$ is the kth zero of ${\displaystyle \psi }$ (see below), and ${\displaystyle \gamma }$ is the Euler–Mascheroni constant.

## Series formula

The digamma function can be computed in the complex plane outside negative integers (Abramowitz and Stegun 6.3.16),[1] using

${\displaystyle \psi (z+1)=-\gamma +\sum _{n=1}^{\infty }{\frac {z}{n(n+z)}}\qquad z\neq -1,-2,-3,\ldots }$

or

${\displaystyle \psi (z)=-\gamma +\sum _{n=0}^{\infty }{\frac {z-1}{(n+1)(n+z)}}=-\gamma +\sum _{n=0}^{\infty }\left({\frac {1}{n+1}}-{\frac {1}{n+z}}\right)\qquad z\neq 0,-1,-2,-3,\ldots }$

This can be utilized to evaluate infinite sums of rational functions, i.e.,

${\displaystyle \sum _{n=0}^{\infty }u_{n}=\sum _{n=0}^{\infty }{\frac {p(n)}{q(n)}},}$

where p(n) and q(n) are polynomials of n.

Performing partial fraction on un in the complex field, in the case when all roots of q(n) are simple roots,

${\displaystyle u_{n}={\frac {p(n)}{q(n)}}=\sum _{k=1}^{m}{\frac {a_{k}}{n+b_{k}}}.}$

For the series to converge,

${\displaystyle \lim _{n\to \infty }nu_{n}=0,}$

otherwise the series will be greater than the harmonic series and thus diverge. Hence

${\displaystyle \sum _{k=1}^{m}a_{k}=0,}$

and

{\displaystyle {\begin{aligned}\sum _{n=0}^{\infty }u_{n}&=\sum _{n=0}^{\infty }\sum _{k=1}^{m}{\frac {a_{k}}{n+b_{k}}}\\&=\sum _{n=0}^{\infty }\sum _{k=1}^{m}a_{k}\left({\frac {1}{n+b_{k}}}-{\frac {1}{n+1}}\right)\\&=\sum _{k=1}^{m}\left(a_{k}\sum _{n=0}^{\infty }\left({\frac {1}{n+b_{k}}}-{\frac {1}{n+1}}\right)\right)\\&=-\sum _{k=1}^{m}a_{k}{\big (}\psi (b_{k})+\gamma {\big )}\\&=-\sum _{k=1}^{m}a_{k}\psi (b_{k}).\end{aligned}}}

With the series expansion of higher rank polygamma function a generalized formula can be given as

${\displaystyle \sum _{n=0}^{\infty }u_{n}=\sum _{n=0}^{\infty }\sum _{k=1}^{m}{\frac {a_{k}}{(n+b_{k})^{r_{k}}}}=\sum _{k=1}^{m}{\frac {(-1)^{r_{k}}}{(r_{k}-1)!}}a_{k}\psi ^{r_{k}-1}(b_{k}),}$

provided the series on the left converges.

## Taylor series

The digamma has a rational zeta series, given by the Taylor series at z = 1. This is

${\displaystyle \psi (z+1)=-\gamma -\sum _{k=1}^{\infty }\zeta (k+1)(-z)^{k}}$,

which converges for |z| < 1. Here, ζ(n) is the Riemann zeta function. This series is easily derived from the corresponding Taylor's series for the Hurwitz zeta function.

## Newton series

The Newton series for the digamma follows from Euler's integral formula:

${\displaystyle \psi (s+1)=-\gamma -\sum _{k=1}^{\infty }{\frac {(-1)^{k}}{k}}{\binom {s}{k}}}$

where (s
k
)
is the binomial coefficient.

## Reflection formula

The digamma function satisfies a reflection formula similar to that of the gamma function:

${\displaystyle \psi (1-x)-\psi (x)=\pi \cot \pi x}$

## Recurrence formula and characterization

The digamma function satisfies the recurrence relation

${\displaystyle \psi (x+1)=\psi (x)+{\frac {1}{x}}.}$

Thus, it can be said to "telescope" 1 / x, for one has

${\displaystyle \Delta [\psi ](x)={\frac {1}{x}}}$

where Δ is the forward difference operator. This satisfies the recurrence relation of a partial sum of the harmonic series, thus implying the formula

${\displaystyle \psi (n)=H_{n-1}-\gamma }$

where γ is the Euler–Mascheroni constant.

More generally, one has

${\displaystyle \psi (x+1)=-\gamma +\sum _{k=1}^{\infty }\left({\frac {1}{k}}-{\frac {1}{x+k}}\right).}$

Actually, ψ is the only solution of the functional equation

${\displaystyle F(x+1)=F(x)+{\frac {1}{x}}}$

that is monotonic on + and satisfies F(1) = −γ. This fact follows immediately from the uniqueness of the Γ function given its recurrence equation and convexity restriction. This implies the useful difference equation:

${\displaystyle \psi (x+N)-\psi (x)=\sum _{k=0}^{N-1}{\frac {1}{x+k}}}$

## Some finite sums involving the digamma function

There are numerous finite summation formulas for the digamma function. Basic summation formulas, such as

${\displaystyle \sum _{r=1}^{m}\psi \left({\frac {r}{m}}\right)=-m(\gamma +\ln m),}$
${\displaystyle \sum _{r=1}^{m}\psi \left({\frac {r}{m}}\right)\cdot \exp {\dfrac {2\pi rki}{m}}=m\ln \left(1-\exp {\frac {2\pi ki}{m}}\right),\qquad k\in \mathbb {Z} \ ,\ m\in \mathbb {N} \ ,\ k\neq m.}$
${\displaystyle \sum _{r=1}^{m-1}\psi \left({\frac {r}{m}}\right)\cdot \cos {\dfrac {2\pi rk}{m}}=m\ln \left(2\sin {\frac {k\pi }{m}}\right)+\gamma ,\qquad k=1,2,\ldots ,m-1}$
${\displaystyle \sum _{r=1}^{m-1}\psi \left({\frac {r}{m}}\right)\cdot \sin {\frac {2\pi rk}{m}}={\frac {\pi }{2}}(2k-m),\qquad k=1,2,\ldots ,m-1}$

are due to Gauss.[4][5] More complicated formulas, such as

${\displaystyle \sum _{r=0}^{m-1}\psi \left({\frac {2r+1}{2m}}\right)\cdot \cos {\frac {(2r+1)k\pi }{m}}=m\ln \left(\tan {\frac {\pi k}{2m}}\right),\qquad k=1,2,\ldots ,m-1}$
${\displaystyle \sum _{r=0}^{m-1}\psi \left({\frac {2r+1}{2m}}\right)\cdot \sin {\dfrac {(2r+1)k\pi }{m}}=-{\frac {\pi m}{2}},\qquad k=1,2,\ldots ,m-1}$
${\displaystyle \sum _{r=1}^{m-1}\psi \left({\frac {r}{m}}\right)\cdot \cot {\frac {\pi r}{m}}=-{\frac {\pi (m-1)(m-2)}{6}}}$
${\displaystyle \sum _{r=1}^{m-1}\psi \left({\frac {r}{m}}\right)\cdot {\frac {r}{m}}=-{\frac {\gamma }{2}}(m-1)-{\frac {m}{2}}\ln m-{\frac {\pi }{2}}\sum _{r=1}^{m-1}{\frac {r}{m}}\cdot \cot {\frac {\pi r}{m}}}$
${\displaystyle \sum _{r=1}^{m-1}\psi \left({\frac {r}{m}}\right)\cdot \cos {\dfrac {(2l+1)\pi r}{m}}=-{\frac {\pi }{m}}\sum _{r=1}^{m-1}{\frac {r\cdot \sin {\dfrac {2\pi r}{m}}}{\cos {\dfrac {2\pi r}{m}}-\cos {\dfrac {(2l+1)\pi }{m}}}},\qquad l\in \mathbb {Z} }$
${\displaystyle \sum _{r=1}^{m-1}\psi \left({\frac {r}{m}}\right)\cdot \sin {\dfrac {(2l+1)\pi r}{m}}=-(\gamma +\ln 2m)\cot {\frac {(2l+1)\pi }{2m}}+\sin {\dfrac {(2l+1)\pi }{m}}\sum _{r=1}^{m-1}{\frac {\ln \sin {\dfrac {\pi r}{m}}}{\cos {\dfrac {2\pi r}{m}}-\cos {\dfrac {(2l+1)\pi }{m}}}},\qquad l\in \mathbb {Z} }$
${\displaystyle \sum _{r=1}^{m-1}\psi ^{2}\left({\frac {r}{m}}\right)=(m-1)\gamma ^{2}+m(2\gamma +\ln 4m)\ln {m}-m(m-1)\ln ^{2}2+{\frac {\pi ^{2}(m^{2}-3m+2)}{12}}+m\sum _{l=1}^{m-1}\ln ^{2}\sin {\frac {\pi l}{m}}}$

are due to works of certain modern authors (see e.g. Appendix B in Blagouchine (2014)[6]).

## Gauss's digamma theorem

For positive integers r and m (r < m), the digamma function may be expressed in terms of Euler's constant and a finite number of elementary functions

${\displaystyle \psi \left({\frac {r}{m}}\right)=-\gamma -\ln(2m)-{\frac {\pi }{2}}\cot \left({\frac {r\pi }{m}}\right)+2\sum _{n=1}^{\left\lfloor {\frac {m-1}{2}}\right\rfloor }\cos \left({\frac {2\pi nr}{m}}\right)\ln \sin \left({\frac {\pi n}{m}}\right)}$

which holds, because of its recurrence equation, for all rational arguments.

## Computation and approximation

According to the Euler–Maclaurin formula applied to[7]

${\displaystyle \sum _{n=1}^{x}{\frac {1}{n}}}$

the digamma function for x, a real number, can be approximated by

${\displaystyle \psi (x)\approx \ln(x)-{\frac {1}{2x}}-{\frac {1}{12x^{2}}}+{\frac {1}{120x^{4}}}-{\frac {1}{252x^{6}}}+{\frac {1}{240x^{8}}}-{\frac {5}{660x^{10}}}+{\frac {691}{32760x^{12}}}-{\frac {1}{12x^{14}}}}$

which is the beginning of the asymptotical expansion of ψ(x). The full asymptotic series of this expansions is

${\displaystyle \psi (x)\sim \ln(x)-{\frac {1}{2x}}+\sum _{n=1}^{\infty }{\frac {\zeta (1-2n)}{x^{2n}}}=\ln(x)-{\frac {1}{2x}}-\sum _{n=1}^{\infty }{\frac {B_{2n}}{2nx^{2n}}}}$

where Bk is the kth Bernoulli number and ζ is the Riemann zeta function. Although the infinite sum does not converge for any x, this expansion becomes more accurate for larger values of x and any finite partial sum cut off from the full series. To compute ψ(x) for small x, the recurrence relation

${\displaystyle \psi (x+1)={\frac {1}{x}}+\psi (x)}$

can be used to shift the value of x to a higher value. Beal[8] suggests using the above recurrence to shift x to a value greater than 6 and then applying the above expansion with terms above x14 cut off, which yields "more than enough precision" (at least 12 digits except near the zeroes).

As x goes to infinity, ψ(x) gets arbitrarily close to both ln(x − 1/2) and ln x. Going down from x + 1 to x, ψ decreases by 1 / x, ln(x − 1/2) decreases by ln (x + 1/2) / (x − 1/2), which is more than 1 / x, and ln x decreases by ln (1 + 1 / x), which is less than 1 / x. From this we see that for any positive x greater than 1/2,

${\displaystyle \psi (x)\in \left(\ln \left(x-{\tfrac {1}{2}}\right),\ln x\right)}$

or, for any positive x,

${\displaystyle \exp \psi (x)\in \left(x-{\tfrac {1}{2}},x\right).}$

The exponential exp ψ(x) is approximately x − 1/2 for large x, but gets closer to x at small x, approaching 0 at x = 0.

For x < 1, we can calculate limits based on the fact that between 1 and 2, ψ(x) ∈ [−γ, 1 − γ], so

${\displaystyle \psi (x)\in \left(-{\frac {1}{x}}-\gamma ,1-{\frac {1}{x}}-\gamma \right),\quad x\in (0,1)}$

or

${\displaystyle \exp \psi (x)\in \left(\exp \left(-{\frac {1}{x}}-\gamma \right),e\exp \left(-{\frac {1}{x}}-\gamma \right)\right).}$

From the above asymptotic series for ψ, one can derive an asymptotic series for exp(−ψ(x)). The series matches the overall behaviour well, that is, it behaves asymptotically as it should for large arguments, and has a zero of unbounded multiplicity at the origin too.

${\displaystyle {\frac {1}{\exp \psi (x)}}\sim {\frac {1}{x}}+{\frac {1}{2\cdot x^{2}}}+{\frac {5}{4\cdot 3!\cdot x^{3}}}+{\frac {3}{2\cdot 4!\cdot x^{4}}}+{\frac {47}{48\cdot 5!\cdot x^{5}}}-{\frac {5}{16\cdot 6!\cdot x^{6}}}+\cdots }$

This is similar to a Taylor expansion of exp(−ψ(1 / y)) at y = 0, but it does not converge.[9] (The function is not analytic at infinity.) A similar series exists for exp(ψ(x)) which starts with ${\displaystyle \exp \psi (x)\sim x-{\frac {1}{2}}.}$

If one calculates the asymptotic series for ψ(x+1/2) it turns out that there are no odd powers of x (there is no x−1 term). This leads to the following asymptotic expansion, which saves computing terms of even order.

${\displaystyle \exp \psi \left(x+{\tfrac {1}{2}}\right)\sim x+{\frac {1}{4!\cdot x}}-{\frac {37}{8\cdot 6!\cdot x^{3}}}+{\frac {10313}{72\cdot 8!\cdot x^{5}}}-{\frac {5509121}{384\cdot 10!\cdot x^{7}}}+\cdots }$

## Special values

The digamma function has values in closed form for rational numbers, as a result of Gauss's digamma theorem. Some are listed below:

{\displaystyle {\begin{aligned}\psi (1)&=-\gamma \\\psi \left({\tfrac {1}{2}}\right)&=-2\ln {2}-\gamma \\\psi \left({\tfrac {1}{3}}\right)&=-{\frac {\pi }{2{\sqrt {3}}}}-{\frac {3\ln {3}}{2}}-\gamma \\\psi \left({\tfrac {1}{4}}\right)&=-{\frac {\pi }{2}}-3\ln {2}-\gamma \\\psi \left({\tfrac {1}{6}}\right)&=-{\frac {\pi {\sqrt {3}}}{2}}-2\ln {2}-{\frac {3\ln {3}}{2}}-\gamma \\\psi \left({\tfrac {1}{8}}\right)&=-{\frac {\pi }{2}}-4\ln {2}-{\frac {\pi +\ln \left(2+{\sqrt {2}}\right)-\ln \left(2-{\sqrt {2}}\right)}{\sqrt {2}}}-\gamma .\end{aligned}}}

Moreover, by the series representation, it can easily be deduced that at the imaginary unit,

{\displaystyle {\begin{aligned}\operatorname {Re} \psi (i)&=-\gamma -\sum _{n=0}^{\infty }{\frac {n-1}{n^{3}+n^{2}+n+1}},\\[8px]\operatorname {Im} \psi (i)&=\sum _{n=0}^{\infty }{\frac {1}{n^{2}+1}}\\&={\frac {1}{2}}+{\frac {\pi }{2}}\coth \pi .\end{aligned}}}

## Roots of the digamma function

The roots of the digamma function are the saddle points of the complex-valued gamma function. Thus they lie all on the real axis. The only one on the positive real axis is the unique minimum of the real-valued gamma function on + at x0 = 1.461632144968.... All others occur single between the poles on the negative axis:

{\displaystyle {\begin{aligned}x_{1}&=-0.504\,083\,008...,\\x_{2}&=-1.573\,498\,473...,\\x_{3}&=-2.610\,720\,868...,\\x_{4}&=-3.635\,293\,366...,\\&\qquad \vdots \end{aligned}}}

Already in 1881, Charles Hermite observed[10] that

${\displaystyle x_{n}=-n+{\frac {1}{\ln n}}+O\left({\frac {1}{(\ln n)^{2}}}\right)}$

holds asymptotically. A better approximation of the location of the roots is given by

${\displaystyle x_{n}\approx -n+{\frac {1}{\pi }}\arctan \left({\frac {\pi }{\ln n}}\right)\qquad n\geq 2}$

and using a further term it becomes still better

${\displaystyle x_{n}\approx -n+{\frac {1}{\pi }}\arctan \left({\frac {\pi }{\ln n+{\frac {1}{8n}}}}\right)\qquad n\geq 1}$

which both spring off the reflection formula via

${\displaystyle 0=\psi (1-x_{n})=\psi (x_{n})+{\frac {\pi }{\tan \pi x_{n}}}}$

and substituting ψ(xn) by its not convergent asymptotic expansion. The correct second term of this expansion is 1 / 2n, where the given one works good to approximate roots with small n.

Another improvement of Hermite's formula can be given[3]:

${\displaystyle x_{n}=-n+{\frac {1}{\log n}}-{\frac {1}{2n(\log n)^{2}}}+O\left({\frac {1}{n^{2}(\log n)^{2}}}\right).}$

Regarding the zeros, the following infinite sum identities were recently proved by István Mező and Michael Hoffman[3]

{\displaystyle {\begin{aligned}\sum _{n=0}^{\infty }{\frac {1}{x_{n}^{2}}}&=\gamma ^{2}+{\frac {\pi ^{2}}{2}},\\\sum _{n=0}^{\infty }{\frac {1}{x_{n}^{3}}}&=-4\zeta (3)-\gamma ^{3}-{\frac {\gamma \pi ^{2}}{2}},\\\sum _{n=0}^{\infty }{\frac {1}{x_{n}^{4}}}&=\gamma ^{4}+{\frac {\pi ^{4}}{9}}+{\frac {2}{3}}\gamma ^{2}\pi ^{2}+4\gamma \zeta (3).\end{aligned}}}

In general, the function

${\displaystyle Z(k)=\sum _{n=0}^{\infty }{\frac {1}{x_{n}^{k}}}}$

can be determined and it is studied in detail by the cited authors.

The following results[3]

{\displaystyle {\begin{aligned}\sum _{n=0}^{\infty }{\frac {1}{x_{n}^{2}+x_{n}}}&=-2,\\\sum _{n=0}^{\infty }{\frac {1}{x_{n}^{2}-x_{n}}}&=\gamma +{\frac {\pi ^{2}}{6\gamma }}\\\end{aligned}}}

also hold true.

Here γ is the Euler–Mascheroni constant.

## Regularization

The digamma function appears in the regularization of divergent integrals

${\displaystyle \int _{0}^{\infty }{\frac {dx}{x+a}},}$

this integral can be approximated by a divergent general Harmonic series, but the following value can be attached to the series

${\displaystyle \sum _{n=0}^{\infty }{\frac {1}{n+a}}=-\psi (a).}$