thumb|300px|The digamma function <math>\psi(z)</math>,<br>visualized using [[domain coloring]]
thumb|300px|Plots of the digamma and the next three polygamma functions along the real line (they are real-valued on the real line)
In mathematics, the digamma function is defined as the logarithmic derivative of the gamma function: and it asymptotically behaves as
:<math>\psi(z) \sim \ln{z} - \frac{1}{2z},</math>
for complex numbers with large modulus (<math>|z|\rightarrow\infty</math>) in the sector <math>\left|\arg z\right|<\pi-\varepsilon</math> for any <math>\varepsilon > 0</math>.
The digamma function is often denoted as <math>\psi_0(x), \psi^{(0)}(x) </math> or (the uppercase form of the archaic Greek letter digamma meaning double-gamma).
Relation to harmonic numbers
The gamma function obeys the equation
:<math>\Gamma(z+1)=z\Gamma(z). \, </math>
Taking the logarithm on both sides and using the functional equation property of the log-gamma function gives:
:<math>\log \Gamma(z+1)=\log(z)+\log \Gamma(z), </math>
Differentiating both sides with respect to gives:
:<math>\psi(z+1)=\psi(z)+\frac{1}{z}</math>
Since the harmonic numbers are defined for positive integers as
:<math>H_n=\sum_{k=1}^n \frac 1 k, </math>
the digamma function is related to them by
:<math>\psi(n)=H_{n-1}-\gamma,</math>
where and is the Euler–Mascheroni constant. For half-integer arguments the digamma function takes the values
:<math> \psi \left(n+\tfrac12\right)=-\gamma-2\ln 2 +\sum_{k=1}^n \frac 2 {2k-1} = -\gamma-2\ln 2 + 2H_{2n}-H_n.</math>
Integral representations
If the real part of is positive then the digamma function has the following integral representation due to Gauss:
:<math>\psi(z) = \int_0^\infty \left(\frac{e^{-t{t} - \frac{e^{-zt{1-e^{-t\right)\,dt.</math>
Combining this expression with an integral identity for the Euler–Mascheroni constant <math>\gamma</math> gives:
:<math>\psi(z + 1) = -\gamma + \int_0^1 \left(\frac{1-t^z}{1-t}\right)\,dt.</math>
The integral is Euler's harmonic number <math>H_z</math>, so the previous formula may also be written
:<math>\psi(z + 1) = \psi(1) + H_z.</math>
A consequence is the following generalization of the recurrence relation:
:<math>\psi(w + 1) - \psi(z + 1) = H_w - H_z.</math>
An integral representation due to Dirichlet is:
:<math>\psi(z) = \log z - \frac{1}{2z} - \int_0^\infty \left(\frac{1}{2} - \frac{1}{t} + \frac{1}{e^t - 1}\right)e^{-tz}\,dt.</math>
This formula is also a consequence of Binet's first integral for the gamma function. The integral may be recognized as a Laplace transform.
Binet's second integral for the gamma function gives a different formula for <math>\psi</math> which also gives the first few terms of the asymptotic expansion:
:<math>\psi(z) = \log z - \frac{1}{2z} - 2\int_0^\infty \frac{t\,dt}{(t^2 + z^2)(e^{2\pi t} - 1)}.</math>
From the definition of <math>\psi</math> and the integral representation of the gamma function, one obtains
:<math>\psi(z) = \frac{1}{\Gamma(z)} \int_0^\infty t^{z-1} \ln (t) e^{-t}\,dt,</math>
with <math>\Re z > 0</math>.
Infinite product representation
The function <math>\psi(z)/\Gamma(z)</math> is an entire function, reads
:<math>\begin{align}
\psi(s)
&= -\gamma + (s-1) - \frac{(s-1)(s-2)}{2\cdot2!} + \frac{(s-1)(s-2)(s-3)}{3\cdot3!}\cdots,\quad\Re(s)> 0, \\
&= -\gamma - \sum_{k=1}^\infty \frac{(-1)^k}{k} \binom{s-1}{k}\cdots,\quad\Re(s)> 0.
\end{align}</math>
where is the binomial coefficient. It may also be generalized to
:<math>\psi(s+1) = -\gamma - \frac{1}{m} \sum_{k=1}^{m-1}\frac{m-k}{s+k} - \frac{1}{m}\sum_{k=1}^\infty\frac{(-1)^k}{k}\left\{\binom{s+m}{k+1}-\binom{s}{k+1}\right\},\qquad \Re(s)>-1,</math>
where More complicated formulas, such as
: <math>\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</math>
: <math>\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</math>
:<math>\sum_{r=1}^{m-1} \psi\left(\frac{r}{m}\right)\cdot\cot\frac{\pi r}{m}= -\frac{\pi(m-1)(m-2)}{6}</math>
:<math>\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} </math>
:<math>\sum_{r=1}^{m-1}\psi \left(\frac{r}{m}\right) \cdot\cos\dfrac{(2\ell+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{(2\ell+1)\pi }{m} }, \qquad \ell\in\mathbb{Z} </math>
:<math>\sum_{r=1}^{m-1}\psi \left(\frac{r}{m}\right) \cdot\sin\dfrac{(2\ell+1)\pi r}{m}=-(\gamma+\ln2m)\cot\frac{(2\ell+1)\pi}{2m} + \sin\dfrac{(2\ell+1)\pi }{m}\sum_{r=1}^{m-1} \frac{\ln\sin\dfrac{\pi r}{m {\cos\dfrac{2\pi r}{m} -\cos\dfrac{(2\ell+1)\pi }{m} } , \qquad \ell\in\mathbb{Z}</math>
:<math>\sum_{r=1}^{m-1} \psi^2\left(\frac{r}{m}\right)= (m-1)\gamma^2 + m(2\gamma+\ln4m)\ln{m} -m(m-1)\ln^2 2 +\frac{\pi^2 (m^2-3m+2)}{12} +m\sum_{\ell=1}^{ m-1 } \ln^2 \sin\frac{\pi\ell}{m}</math>
are due to works of certain modern authors (see e.g. Appendix B in Blagouchine (2014)).
We also have
:<math> 1+\frac{1}{2}+\frac{1}{3}+...+\frac{1}{k-1}-\gamma-\ln k=\frac{1}{k}\sum_{n=0}^{k-1}\psi\left(1+\frac{n}{k}\right), k=2,3, ...</math>
Gauss's digamma theorem
For positive integers and (), the digamma function may be expressed in terms of Euler's constant and a finite number of elementary functions
:<math>\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) </math>
which holds, because of its recurrence equation, for all rational arguments.
Multiplication theorem
The multiplication theorem of the <math>\Gamma</math>-function is equivalent to
:<math>\psi(nz)=\frac{1}{n}\sum_{k=0}^{n-1} \psi\left(z+\frac{k}{n}\right) +\ln n .</math>
Asymptotic expansion
The digamma function has the asymptotic expansion
:<math>\psi(z) \sim \ln z + \sum_{n=1}^\infty \frac{\zeta(1-n)}{z^n} = \ln z - \sum_{n=1}^\infty \frac{B_n}{nz^n},</math>
where is the th Bernoulli number and is the Riemann zeta function. The first few terms of this expansion are:
:<math>\psi(z) \sim \ln z - \frac{1}{2z} - \frac{1}{12z^2} + \frac{1}{120z^4} - \frac{1}{252z^6} + \frac{1}{240z^8} - \frac{1}{132z^{10 + \frac{691}{32760z^{12 - \frac{1}{12z^{14 + \cdots.</math>
Although the infinite sum does not converge for any , any finite partial sum becomes increasingly accurate as increases.
The expansion can be found by applying the Euler–Maclaurin formula to the sum
:<math>\sum_{n=1}^\infty \left(\frac{1}{n} - \frac{1}{z + n}\right)</math>
The expansion can also be derived from the integral representation coming from Binet's second integral formula for the gamma function. Expanding <math>t / (t^2 + z^2)</math> as a geometric series and substituting an integral representation of the Bernoulli numbers leads to the same asymptotic series as above. Furthermore, expanding only finitely many terms of the series gives a formula with an explicit error term:
:<math>\psi(z) = \ln z - \frac{1}{2z} - \sum_{n=1}^N \frac{B_{2n{2nz^{2n + (-1)^{N+1}\frac{2}{z^{2N \int_0^\infty \frac{t^{2N+1}\,dt}{(t^2 + z^2)(e^{2\pi t} - 1)}.</math>
Inequalities
When , the function
:<math>\ln x - \frac{1}{2x} - \psi(x)</math>
is completely monotonic and in particular positive. This is a consequence of Bernstein's theorem on monotone functions applied to the integral representation coming from Binet's first integral for the gamma function. Additionally, by the convexity inequality <math>1 + t \le e^t</math>, the integrand in this representation is bounded above by <math>e^{-tz}/2</math>.
:<math>\frac{1}{x} - \ln x + \psi(x)</math>
is also completely monotonic. It follows that, for all ,
:<math>\ln x - \frac{1}{x} \le \psi(x) \le \ln x - \frac{1}{2x}.</math>
This recovers a theorem of Horst Alzer. Alzer also proved that, for ,
:<math>\frac{1 - s}{x + s} < \psi(x + 1) - \psi(x + s),</math>
Related bounds were obtained by Elezovic, Giordano, and Pecaric, who proved that, for ,
:<math>\ln(x + \tfrac{1}{2}) - \frac{1}{x} < \psi(x) < \ln(x + e^{-\gamma}) - \frac{1}{x},</math>
where <math>\gamma=-\psi(1)</math> is the Euler–Mascheroni constant. The constants (<math>0.5</math> and <math>e^{-\gamma}\approx0.56</math>) appearing in these bounds are the best possible.
The mean value theorem implies the following analog of Gautschi's inequality: If , where is the unique positive real root of the digamma function, and if , then
:<math>\exp\left((1 - s)\frac{\psi'(x + 1)}{\psi(x + 1)}\right) \le \frac{\psi(x + 1)}{\psi(x + s)} \le \exp\left((1 - s)\frac{\psi'(x + s)}{\psi(x + s)}\right).</math>
Moreover, equality holds if and only if .
Inspired by the harmonic mean value inequality for the classical gamma function, Horzt Alzer and Graham Jameson proved, among other things, a harmonic mean-value inequality for the digamma function:
<math> -\gamma \leq \frac{2 \psi(x) \psi(\frac{1}{x})}{\psi(x)+\psi(\frac{1}{x})} </math> for <math>x>0</math>
Equality holds if and only if <math>x=1</math>.
Computation and approximation
The asymptotic expansion gives an easy way to compute when the real part of is large. To compute for small , the recurrence relation
:<math> \psi(x+1) = \frac{1}{x} + \psi(x)</math>
can be used to shift the value of to a higher value. Beal suggests using the above recurrence to shift to a value greater than 6 and then applying the above expansion with terms above cut off, which yields "more than enough precision" (at least 12 digits except near the zeroes).
As goes to infinity, gets arbitrarily close to both and . Going down from to , decreases by , decreases by , which is more than , and decreases by , which is less than . From this we see that for any positive greater than ,
:<math>\psi(x)\in \left(\ln\left(x-\tfrac12\right), \ln x\right)</math>
or, for any positive ,
:<math>\exp \psi(x)\in\left(x-\tfrac12,x\right).</math>
The exponential is approximately for large , but gets closer to at small , approaching 0 at .
For , we can calculate limits based on the fact that between 1 and 2, , so
:<math>\psi(x)\in\left(-\frac{1}{x}-\gamma, 1-\frac{1}{x}-\gamma\right),\quad x\in(0, 1)</math>
or
:<math>\exp \psi(x)\in\left(\exp\left(-\frac{1}{x}-\gamma\right),e\exp\left(-\frac{1}{x}-\gamma\right)\right).</math>
From the above asymptotic series for , one can derive an asymptotic series for . 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.
:<math> \frac{1}{\exp \psi(x)} \sim \frac{1}{x}+\frac{1}{2\cdot x^2}+\frac{5}{4\cdot3!\cdot x^3}+\frac{3}{2\cdot4!\cdot x^4}+\frac{47}{48\cdot5!\cdot x^5} - \frac{5}{16\cdot6!\cdot x^6} + \cdots</math>
This is similar to a Taylor expansion of at , but it does not converge. (The function is not analytic at infinity.) A similar series exists for which starts with <math>\exp \psi(x) \sim x- \frac 12.</math>
If one calculates the asymptotic series for it turns out that there are no odd powers of (there is no <sup>−1</sup> term). This leads to the following asymptotic expansion, which saves computing terms of even order.
:<math> \exp \psi\left(x+\tfrac{1}{2}\right) \sim x + \frac{1}{4!\cdot x} - \frac{37}{8\cdot6!\cdot x^3} + \frac{10313}{72\cdot8!\cdot x^5} - \frac{5509121}{384\cdot10!\cdot x^7} + \cdots</math>
Similar in spirit to the Lanczos approximation of the <math>\Gamma</math>-function is Spouge's approximation.
Another alternative is to use the recurrence relation or the multiplication formula to shift the argument of <math>\psi(x)</math> into the range <math>1\le x\le 3</math> and to evaluate the Chebyshev series there.
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:
:<math>\begin{align}
\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 (\sqrt{2} + 1 \right ) - \ln \left (\sqrt{2} - 1 \right ) }{\sqrt{2 - \gamma.
\end{align}</math>
Moreover, by taking the logarithmic derivative of <math>|\Gamma (bi)|^2</math> or <math>|\Gamma (\tfrac{1}{2}+bi)|^2</math> where <math>b</math> is real-valued, it can easily be deduced that
:<math>\operatorname{Im} \psi(bi) = \frac{1}{2b}+\frac{\pi}{2}\coth (\pi b),</math>
:<math>\operatorname{Im} \psi(\tfrac{1}{2}+bi) = \frac{\pi}{2}\tanh (\pi b).</math>
Apart from Gauss's digamma theorem, no such closed formula is known for the real part in general. We have, for example, at the imaginary unit the numerical approximation
:<math>\operatorname{Re} \psi(i) = -\gamma-\sum_{n=0}^\infty\frac{n-1}{n^3+n^2+n+1} \approx 0.09465.</math>
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 . All others occur single between the poles on the negative axis:
:
:
:
:
:<math>\vdots</math>
Already in 1881, Charles Hermite observed that
:<math>x_n = -n + \frac{1}{\ln n} + O\left(\frac{1}{(\ln n)^2}\right)</math>
holds asymptotically. A better approximation of the location of the roots is given by
:<math>x_n \approx -n + \frac{1}{\pi}\arctan\left(\frac{\pi}{\ln n}\right)\qquad n \ge 2</math>
and using a further term it becomes still better
:<math>x_n \approx -n + \frac{1}{\pi}\arctan\left(\frac{\pi}{\ln n + \frac{1}{8n\right)\qquad n \ge 1</math>
which both spring off the reflection formula via
:<math>0 = \psi(1-x_n) = \psi(x_n) + \frac{\pi}{\tan \pi x_n}</math>
and substituting by its not convergent asymptotic expansion. The correct second term of this expansion is , where the given one works well to approximate roots with small .
Another improvement of Hermite's formula can be given:
:<math>\begin{align}
\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} + \frac23 \gamma^2 \pi^2 + 4\gamma\zeta(3).
\end{align}</math>
In general, the function
:<math>
Z(k)=\sum_{n=0}^\infty\frac{1}{x_n^k}
</math>
can be determined and it is studied in detail by the cited authors.
The following results
</references>
External links
- —psi(1/2)
: psi(1/3), psi(2/3), psi(1/4), psi(3/4), to psi(1/5) to psi(4/5).
- Implementation in the GNU Scientific library
- C function with the Chebyshev series
- Java implementation with Taylor series