Franklin’s Blog

The Gamma and Digamma functions

2018 March 10

Before I get into the definitions and the properties of these two functions, I must first introduce a new mathematical constant: the Euler-Mascheroni constant
\[\gamma\approx 0.57721566...\]
This constant is defined as follows:
\[\gamma := \lim_{n\to\infty} (H_n-\ln(n))\]
where \(H_n\) are the harmonic numbers.

I will now establish an integral representation of the Euler-Mascheroni constant which will prove to be useful later on. Recall that the harmonic numbers \(H_n\) are given by
\[H_n=\int_0^1 \frac{1-x^n}{1-x}dx=\int_0^1\frac{1-(1-x)^n}{x}dx\]
and recall that
\[\ln(n)=\int_{1/n}^1 \frac{dx}{x}\]
Recall also the value of the limit
\[\lim_{n\to\infty}\bigg(1+\frac{x}{n}\bigg)^n=e^x\]
which will come into play shortly.

By definition, we have
\[\begin{align}\gamma &= \lim_{n\to\infty} \bigg(\int_0^1\frac{1-(1-x)^n}{x}dx-\int_{1/n}^1 \frac{dx}{x}\bigg)\\ &= \lim_{n\to\infty} \bigg(\int_0^{1/n}\frac{1-(1-x)^n}{x}dx+\int_{1/n}^1\frac{1-(1-x)^n}{x}dx-\int_{1/n}^1 \frac{dx}{x}\bigg)\\ &= \lim_{n\to\infty} \bigg(\int_0^{1/n}\frac{1-(1-x)^n}{x}dx-\int_{1/n}^1\frac{(1-x)^n}{x}dx\bigg)\\ &= \lim_{n\to\infty} \bigg(\int_0^{1}\frac{1-(1-\frac{x}{n})^n}{x}dx-\int_{1}^n\frac{(1-\frac{x}{n})^n}{x}dx\bigg)\\ &= \lim_{n\to\infty} \bigg(\int_0^{1}\frac{1-e^{-x}}{x}dx-\int_{1}^n\frac{e^{-x}}{x}dx\bigg)\\ &= \int_0^{1}\frac{1-e^{-x}}{x}dx-\int_{1}^\infty \frac{e^{-x}}{x}dx\\ \end{align}\]

and so

\[\bbox[lightgray,5px] { \gamma=\int_0^{1}\frac{1-e^{-x}}{x}dx-\int_{1}^\infty \frac{e^{-x}}{x}dx }\]

This is the integral representation that I was looking for! Now I will get into the definitions and properties of the gamma and digamma functions, and the Euler-Mascheroni constant will show up eventually - just wait.

In my previous post about the Gamma and Lambert-W functions, I introduced the Gamma function, but did not explore its properties in much depth. Here is its definition:

\[\Gamma(z):=\int_0^\infty x^{z-1}e^{-x}dx\]

For positive integer \(z\), it holds that
\[\Gamma(z)=(z-1)!\]
and, in general, the gamma function has the property
\[\Gamma(z+1)=z\Gamma(z)\]
which I demonstrated in my other post by integrating by parts. It also has the special value
\[\Gamma(1/2)=\sqrt{\pi}\]
which is the famed “Gaussian integral.”

Before moving on to the digamma function, I would like to prove the following property of the gamma function, relating it to the Beta function \(\mathrm{B}(a,b)\):
\[\mathrm{B}(a,b):=\int_{0}^1 x^{a-1}(1-x)^{b-1}dx=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}\]
Here is the proof:
\[\begin{align} \Gamma(a)\Gamma(b) &=\int_0^\infty \int_0^\infty x^{a-1}y^{b-1}e^{-x-y}dydx\\ &=4\int_0^\infty \int_0^\infty x^{2a-1}y^{2b-1}e^{-x^2-y^2}dydx\\ &=4\int_0^{\pi/2} \int_0^\infty (r\cos\theta)^{2a-1}(r\sin\theta)^{2b-1}e^{-r^2} rdrd\theta\\ &=4\int_0^{\pi/2} \int_0^\infty r^{2a+2b-1}\cos^{2a-1}(\theta)\sin^{2b-1}(\theta)e^{-r^2} drd\theta\\ &=2\int_0^{\pi/2} \int_0^\infty r^{a+b-1}\cos^{2a-1}(\theta)\sin^{2b-1}(\theta)e^{-r} drd\theta\\ &=2\bigg(\int_0^{\pi/2} \cos^{2a-1}(\theta)\sin^{2b-1}(\theta) d\theta\bigg) \bigg(\int_0^\infty r^{a+b-1}e^{-r}dr\bigg)\\ &=2\bigg(\int_0^{\pi/2} \cos^{2a-1}(\theta)\sin^{2b-1}(\theta) d\theta\bigg) \Gamma(a+b)\\ &=2\Gamma(a+b)\int_0^{1} \cos^{2a-1}(\arcsin x)\sin^{2b-1}(\arcsin x) \frac{dx}{\sqrt{1-x^2}}\\ &=2\Gamma(a+b)\int_0^{1} (1-x^2)^{a-1}x^{2b-1} dx\\ &=\Gamma(a+b)\int_0^{1} (1-x)^{a-1}x^{b-1} dx\\ &=\Gamma(a+b)\mathrm{B}(a,b)\\ \end{align}\]
and so
\[\Gamma(a)\Gamma(b)=\Gamma(a+b)\mathrm{B}(a,b)\]
and
\[\bbox[lightgray,5px] { \mathrm{B}(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)} }\]

This result can be used to prove Legendre’s famous duplication formula for the Gamma Function, which expresses \(\Gamma(2z)\) in terms of \(\Gamma(z)\) and \(\Gamma(z+1/2)\). Recall the definition of the Beta function:
\[\mathrm{B}(a,b):=\int_{0}^1 x^{a-1}(1-x)^{b-1}dx\]
By substituting \(x\to x^2\) in the definite integral, this is equivalent to
\[\mathrm{B}(a,b)=2\int_{0}^1 x^{2a-1}(1-x^2)^{b-1}dx=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}\]
Notice also that, by substituting \(x\to \frac{1+x}{2}\) in the original integral, we have
\[\mathrm{B}(a,b)=\frac{1}{2}\int_{-1}^1 \bigg(\frac{1+x}{2}\bigg)^{a-1}\bigg(\frac{1-x}{2}\bigg)^{b-1}dx\]
and, as a consequence,
\[\begin{align} \mathrm{B}(z,z) &=\frac{1}{2}\int_{-1}^1 \bigg(\frac{1+x}{2}\bigg)^{z-1}\bigg(\frac{1-x}{2}\bigg)^{z-1}dx\\ &=\frac{1}{2^{2z-1}}\int_{-1}^1 (1+x)^{z-1} (1-x)^{z-1}dx\\ &=2^{1-2z}\int_{-1}^1 (1-x^2)^{z-1}dx\\ &=2^{1-2z}\cdot 2\int_{0}^1 (1-x^2)^{z-1}dx\\ &=2^{1-2z}\mathrm{B}(1/2,z)\\ \end{align}\]
Thus, we have
\[\mathrm{B}(z,z)=2^{1-2z}\mathrm{B}(1/2,z)\]
Now, by using the representation of the Beta function in terms of the Gamma function, we have
\[\frac{\Gamma(z)\Gamma(z)}{\Gamma(2z)}=\frac{2^{1-2z}\Gamma(1/2)\Gamma(z)}{\Gamma(z+1/2)}\]
which can be rearranged to obtain

\[\bbox[lightgray,5px] { \Gamma(2z)=\frac{2^{2z-1}\Gamma(z)\Gamma(z+1/2)}{\sqrt \pi} }\]

which is the final form of the Gamma function’s duplication formula.

Now I will introduce the Digamma function, defined as
\[\psi(z):=\frac{d}{dz}\ln \Gamma(z)=\frac{\Gamma'(z)}{\Gamma(z)}\]
I will start by calculating a few values of the digamma function - this is where the Euler-Mascheroni constant starts to appear. Since the Gamma function is defined as
\[\Gamma(z):=\int_0^\infty x^{z-1}e^{-x}dx\]
we have
\[\Gamma'(z)=\int_0^\infty x^{z-1}e^{-x}\ln(x)dx\]
and so
\[\psi(z)=\Gamma(z)\int_0^\infty x^{z-1}e^{-x}\ln(x)dx\]
Let us first endeavour to calculate the value of \(\psi(1)\):
\[\psi(1)=\int_0^\infty e^{-x}\ln(x)dx\]
Using integration by parts, we have
\[\begin{align} \psi(1) &=\int_0^\infty e^{-x}\ln(x)dx\\ &=\int_0^1 e^{-x}\ln(x)dx+\int_1^\infty e^{-x}\ln(x)dx\\ &=\bigg[(1-e^{-x})\ln(x)\bigg]_0^1-\int_0^1 \frac{1-e^{-x}}{x}dx+\int_1^\infty e^{-x}\ln(x)dx\\ &=-\int_0^1 \frac{1-e^{-x}}{x}dx+\int_1^\infty e^{-x}\ln(x)dx\\ &=-\int_0^1 \frac{1-e^{-x}}{x}dx+\bigg[-e^{-x}\ln(x)\bigg]_1^\infty +\int_1^\infty \frac{e^{-x}}{x}dx\\ &=-\int_0^1 \frac{1-e^{-x}}{x}dx+\int_1^\infty \frac{e^{-x}}{x}dx\\ &=-\int_0^1 \frac{1-e^{-x}}{x}dx+\int_1^\infty \frac{e^{-x}}{x}dx\\ &=-\gamma\\ \end{align}\]

The last step in this equality is justified by the integral representation of \(\gamma\) that we obtained earlier. Thus, we have
\[\psi(1)=-\gamma\]

We could go back and calculate \(\psi(2)\), \(\psi(3)\), and so on, but there’s a shortcut to do this. Remember this recurrence for the gamma function that I mentioned earlier?
\[\Gamma(z+1)=z\Gamma(z)\]
By taking the natural logarithm of both sides, we have
\[\ln\Gamma(z+1)=\ln\Gamma(z)+\ln(z)\]
and by differentiating both sides with respect to \(z\),
\[\psi(z+1)=\psi(z)+\frac{1}{z}\]
This tells us that, for any positive integer \(n\),
\[\psi(n)=-\gamma+H_{n-1}\]
where \(H_n\) are the harmonic numbers.

Observe the following. Since
\[\psi(z+1)-\psi(z)=\frac{1}{z}\]
as previously proven, we have that
\[\sum_{n=1}^N \frac{1}{n}=\psi(N+1)-\psi(1)=\psi(N+1)+\gamma\]
Furthermore, we also have that
\[\sum_{n=1}^N \frac{1}{n+z}=\psi(N+z+1)-\psi(z+1)\]
By subtracting these two equations, we have
\[\sum_{n=1}^N \frac{1}{n}-\frac{1}{n+z}=\psi(N+1)+\gamma-\psi(N+z+1)+\psi(z+1)\]
\[\sum_{n=1}^N \frac{1}{n}-\frac{1}{n+z}=(\psi(N+1)-\psi(N+z+1))+\psi(z+1)+\gamma\]
Now we let \(N\to\infty\). As \(N\) grows large, the term on the RHS inside of the parentheses approaches zero, and so we are left with
\[\sum_{n=1}^\infty \frac{1}{n}-\frac{1}{n+z}=\psi(z+1)+\gamma\]
Thus, we are left with the series representation for the digamma function:

\[\bbox[lightgray,5px] { \psi(z)=-\gamma+\sum_{n=0}^\infty\frac{1}{n+1}-\frac{1}{n+z} }\]

Believe it or not, a closed form for the digamma function can be found for rational arguments. Though its derivation is a bit of a slog, I shall power through it for you. Using the series representation that we just derived, we have
\[\begin{align} \psi(p/q) &= -\gamma+\sum_{n=0}^\infty \frac{1}{n+1}-\frac{1}{n+\frac{p}{q}}\\ &= -\gamma+\sum_{n=0}^\infty \frac{1}{n+1}-\frac{q}{qn+p}\\ &= -\gamma+\sum_{n=0}^\infty\frac{1}{n+1}-\frac{q}{qn+p}\\ &= -\gamma+\lim_{x\to 1^-}\sum_{n=0}^\infty\frac{x^{n}}{n+1}-\frac{qx^{n}}{qn+p}\\ \end{align}\]
Now recall my old blog post in which I devised a strategy for evaluating infinite series in the form
\[\sum_{n=0}^\infty \frac{x^{qn+p}}{qn+p}\]
for natural numbers \(p\) and \(q\). I will not explain the methods I use again in this post, so if you do not understand, I encourage you to go back to the old post. Using the complex roots of unity \(\omega_n\), defined as
\[\omega_n:=e^{2\pi i/n}=\cos\frac{2\pi}{n}+i\sin\frac{2\pi}{n}\]
and using the fact that
\[\sum_{n=0}^\infty \frac{x^{n}}{n+1}=-\frac{\ln(1-x)}{x}\]
one can derive the following formula, for natural numbers \(q\) and \(p\) with \(p\lt q\):
\[\sum_{n=0}^\infty \frac{x^{qn+p-1}}{qn+p}=-\frac{1}{q}\sum_{d=0}^{q-1} \frac{\omega_q^{(1-p)d}\ln(1-\omega_q^d x)}{\omega_q^d x}\]
Thus, it follows from this that
\[\sum_{n=0}^\infty \frac{x^{n}}{qn+p}=-\frac{x^{(1-p)/q}}{q}\sum_{d=0}^{q-1} \frac{\omega_q^{(1-p)d}\ln(1-\omega_q^d x^{1/q})}{\omega_q^d x^{1/q}}\]
and so, continuing the chain of equalities that I began earlier,
\[\begin{align} \psi(p/q) &= -\gamma+\lim_{x\to 1^-}\sum_{n=0}^\infty\frac{x^{n}}{n+1}-\frac{qx^{n}}{qn+p}\\ &= -\gamma+\lim_{x\to 1^-}\bigg[-\frac{\ln(1-x)}{x}+x^{(1-p)/q}\sum_{d=0}^{q-1} \frac{\omega_q^{(1-p)d}\ln(1-\omega_q^d x^{1/q})}{\omega_q^d x^{1/q}}\bigg]\\ &= -\gamma+\lim_{x\to 1^-}\bigg[-\frac{\ln(1-x)}{x}+x^{(1-p)/q}\ln(1-x^{1/q})+x^{(1-p)/q}\sum_{d=1}^{q-1} \frac{\omega_q^{(1-p)d}\ln(1-\omega_q^d x^{1/q})}{\omega_q^d x^{1/q}}\bigg]\\ &= -\gamma+\lim_{x\to 1^-}\bigg[-\frac{\ln(1-x)}{1}+1^{(1-p)/q}\ln(1-x^{1/q})+x^{(1-p)/q}\sum_{d=1}^{q-1} \frac{\omega_q^{(1-p)d}\ln(1-\omega_q^d x^{1/q})}{\omega_q^d x^{1/q}}\bigg]\\ &= -\gamma+\lim_{x\to 1^-}\bigg[-\ln(1-x)+\ln(1-x^{1/q})+x^{(1-p)/q}\sum_{d=1}^{q-1} \frac{\omega_q^{(1-p)d}\ln(1-\omega_q^d x^{1/q})}{\omega_q^d x^{1/q}}\bigg]\\ &= -\gamma+\lim_{x\to 1^-}\bigg[-\ln((1-x^{1/q})(1+x^{1/q}+x^{2/q}+...+x^{(q-1)/q}))+\ln(1-x^{1/q})+x^{(1-p)/q}\sum_{d=1}^{q-1} \frac{\omega_q^{(1-p)d}\ln(1-\omega_q^d x^{1/q})}{\omega_q^d x^{1/q}}\bigg]\\ &= -\gamma+\lim_{x\to 1^-}\bigg[-\ln(1+x^{1/q}+x^{2/q}+...+x^{(q-1)/q})+x^{(1-p)/q}\sum_{d=1}^{q-1} \frac{\omega_q^{(1-p)d}\ln(1-\omega_q^d x^{1/q})}{\omega_q^d x^{1/q}}\bigg]\\ &= -\gamma-\ln(1+1^{1/q}+1^{2/q}+...+1^{(q-1)/q})+\lim_{x\to 1^-}\bigg[x^{(1-p)/q}\sum_{d=1}^{q-1} \frac{\omega_q^{(1-p)d}\ln(1-\omega_q^d x^{1/q})}{\omega_q^d x^{1/q}}\bigg]\\ &= -\gamma-\ln(q)+\lim_{x\to 1^-}\bigg[x^{(1-p)/q}\sum_{d=1}^{q-1} \frac{\omega_q^{(1-p)d}\ln(1-\omega_q^d x^{1/q})}{\omega_q^d x^{1/q}}\bigg]\\ &= -\gamma-\ln(q)+\sum_{d=1}^{q-1} \frac{\omega_q^{(1-p)d}\ln(1-\omega_q^d)}{\omega_q^d}\\ &= -\gamma-\ln(q)+\sum_{d=1}^{q-1} \omega_q^{-pd}\ln(1-\omega_q^d)\\ \end{align}\]

Yuck, look at that nasty algebra. Just look away for a moment and take a little break. Here are a few more formulas that I am going to use in just a moment:
\[\ln(a+bi)=\frac{\ln(a^2+b^2)}{2}+i\arctan\frac{b}{a}\]
\[\cot \frac{x}{2}=\frac{\sin(x)}{1-\cos(x)}\]
\[2-2\cos(2x)=4\sin^2(x)\]
\[\arctan \cot x=\frac{\pi}{2}-x,\space\space\space x\in [0,\pi)\]
\[\sum_{d=1}^{q-1} \cos\frac{2\pi p d}{q}=-1,\space\space\space p,q\in\mathbb N\]
\[\sum_{d=1}^{q-1} \sin\frac{2\pi p d}{q}=0,\space\space\space p,q\in\mathbb N\]
\[\sum_{d=1}^{q-1} \frac{\pi d}{q}\sin\frac{2\pi p d}{q}=-\frac{\pi}{2}\cot\frac{\pi p}{q},\space\space\space p,q\in\mathbb N\]

These formulas are all pretty easy to derive, so I won’t derive them for you. Keep an eye out for them as I continue the derivation - I’m not going to point them out.

Ready? Here we go.

\[\begin{align} \psi(p/q) &= -\gamma-\ln(q)+\sum_{d=1}^{q-1} \omega_q^{-pd}\ln(1-\omega_q^d)\\ &= -\gamma-\ln(q)+\sum_{d=1}^{q-1} \bigg(\cos \frac{2\pi pd}{q}-i\sin \frac{2\pi pd}{q}\bigg)\ln\bigg(1-\cos \frac{2\pi d}{q}-i\sin \frac{2\pi d}{q}\bigg)\\ &= -\gamma-\ln(q)+\sum_{d=1}^{q-1} \bigg(\cos \frac{2\pi pd}{q}-i\sin \frac{2\pi pd}{q}\bigg) \bigg(\frac{\ln((1-\cos(2\pi d/q))^2+\sin^2 (2\pi d/q))}{2}-i\arctan \frac{\sin(2\pi d/q)}{1-\cos(2\pi d/q)} \bigg)\\ &= -\gamma-\ln(q)+\sum_{d=1}^{q-1} \bigg(\cos \frac{2\pi pd}{q}-i\sin \frac{2\pi pd}{q}\bigg) \bigg(\frac{\ln(2-2\cos(2\pi d/q))}{2}-i\arctan \frac{\sin(2\pi d/q)}{1-\cos(2\pi d/q)} \bigg)\\ &= -\gamma-\ln(q)+\sum_{d=1}^{q-1} \bigg(\cos \frac{2\pi pd}{q}-i\sin \frac{2\pi pd}{q}\bigg) \bigg(\frac{\ln(4\sin^2(\pi d/q))}{2}-i\arctan \frac{\sin(2\pi d/q)}{1-\cos(2\pi d/q)} \bigg)\\ &= -\gamma-\ln(q)+\sum_{d=1}^{q-1} \bigg(\cos \frac{2\pi pd}{q}-i\sin \frac{2\pi pd}{q}\bigg) \bigg(\ln(2\sin(\pi d/q))-i\arctan \frac{\sin(2\pi d/q)}{1-\cos(2\pi d/q)} \bigg)\\ \end{align}\]
Because \(\psi(p/q)\) is a real number for rational \(p/q\), it follows that the right side of the equality is a real number, and so we may ignore the imaginary parts of each term of the summation, since they are guaranteed to cancel out to zero. Thus, we have
\[\begin{align} \psi(p/q) &= -\gamma-\ln(q)+\sum_{d=1}^{q-1} \Re\bigg[\bigg(\cos \frac{2\pi pd}{q}-i\sin \frac{2\pi pd}{q}\bigg) \bigg(\ln(2\sin(\pi d/q))-i\arctan \frac{\sin(2\pi d/q)}{1-\cos(2\pi d/q)} \bigg)\bigg]\\ &= -\gamma-\ln(q)+\sum_{d=1}^{q-1} \bigg(\cos \frac{2\pi pd}{q} \ln(2\sin(\pi d/q))-\sin \frac{2\pi pd}{q} \arctan \frac{\sin(2\pi d/q)}{1-\cos(2\pi d/q)}\bigg)\\ &= -\gamma-\ln(q)+\sum_{d=1}^{q-1} \bigg(\cos \frac{2\pi pd}{q} \ln(2\sin(\pi d/q))-\sin \frac{2\pi pd}{q} \arctan \cot\frac{\pi d}{q}\bigg)\\ &= -\gamma-\ln(q)+\sum_{d=1}^{q-1} \bigg(\cos \frac{2\pi pd}{q} \ln(2\sin(\pi d/q))-\sin \frac{2\pi pd}{q} \arctan \cot\frac{\pi d}{q}\bigg)\\ &= -\gamma-\ln(q)+\sum_{d=1}^{q-1}\cos \frac{2\pi pd}{q} \ln(2)+\sum_{d=1}^{q-1} \bigg(\cos \frac{2\pi pd}{q} \ln(\sin(\pi d/q))-\sin \frac{2\pi pd}{q} \arctan \cot\frac{\pi d}{q}\bigg)\\ &= -\gamma-\ln(q)-\ln(2)+\sum_{d=1}^{q-1} \bigg(\cos \frac{2\pi pd}{q} \ln(\sin(\pi d/q))-\sin \frac{2\pi pd}{q} \arctan \cot\frac{\pi d}{q}\bigg)\\ &= -\gamma-\ln(2q)+\sum_{d=1}^{q-1} \bigg(\cos \frac{2\pi pd}{q} \ln(\sin(\pi d/q))-\sin \frac{2\pi pd}{q} \bigg(\frac{\pi}{2}-\frac{\pi d}{q}\bigg)\bigg)\\ &= -\gamma-\ln(2q)+\sum_{d=1}^{q-1}\sin \frac{2\pi pd}{q} \bigg(\frac{\pi d}{q}-\frac{\pi}{2}\bigg)+\sum_{d=1}^{q-1} \bigg(\cos \frac{2\pi pd}{q} \ln(\sin(\pi d/q))\bigg)\\ &= -\gamma-\ln(2q)-\frac{\pi}{2}\cot\frac{\pi p}{q}+\sum_{d=1}^{q-1} \cos \frac{2\pi pd}{q} \ln \sin \frac{\pi d}{q}\\ \end{align}\]

This is it! This is the result we’ve been working towards!

\[\bbox[lightgray,5px] { \psi(p/q)= -\gamma-\ln(2q)-\frac{\pi}{2}\cot\frac{\pi p}{q}+\sum_{d=1}^{q-1} \cos \frac{2\pi pd}{q} \ln \sin \frac{\pi d}{q} }\]

Of course, this only works for \(p/q\in (0,1)\). However, if we want it to work for all rational numbers, we can simply use its recurrence
\[\psi(z+1)-\psi(z)=\frac{1}{z}\]
to say that

\[\psi(n+p/q)= -\gamma-\ln(2q)-\frac{\pi}{2}\cot\frac{\pi p}{q}+\sum_{d=1}^{q-1} \cos \frac{2\pi pd}{q} \ln \sin \frac{\pi d}{q}+\sum_{k=0}^{n-1} \frac{1}{k+\frac{p}{q}}\]

or

\[\psi(-n+p/q)= -\gamma-\ln(2q)-\frac{\pi}{2}\cot\frac{\pi p}{q}+\sum_{d=1}^{q-1} \cos \frac{2\pi pd}{q} \ln \sin \frac{\pi d}{q}+\sum_{k=0}^{n-1} \frac{1}{k-\frac{p}{q}}\]

Using these two formulas, it is easy to verify that, for any rational number \(r\),

\[\psi(1-r)-\psi(r)=\pi\cos(\pi r)\]

Since the digamma function is continuous except at nonpositive integers, this property must necessarily extend to all real numbers, meaning that

\[\bbox[lightgray,5px] { \psi(1-z)-\psi(z)=\pi\cos(\pi z) }\]

By integrating both sides of this with respect to \(z\), we get
\[\ln\Gamma(1-z)+\ln\Gamma(z)=-\ln(\sin z)+C\]
By taking \(z=1/2\), we can solve and obtain \(C=\ln\pi\), giving us
\[\ln\Gamma(1-z)+\ln\Gamma(z)=-\ln(\sin z)+\ln\pi\]
and, by exponentiating both sides, we get the famous Gamma function reflection formula:
\[\bbox[lightgray,5px] { \Gamma(z)\Gamma(1-z)=\frac{\pi}{\sin(\pi z)} }\]

I conclude this post with an exercise for the reader. Using this reflection formula, can you calculate the value of the following integral?

\[\int_0^1 \ln \Gamma(z)dz\]

NOTE: This is called Raabe’s Integral, and the more general form
\[\int_a^{a+1} \ln \Gamma(z)dz\]
can also be evaluated.