Relation between harmonic series $H(m)$ and polygamma function? I have the following formula:
$$h(x)=\sum_{R=1}^m \frac{1}{R+1}-1$$
I have re-expressed this (correctly, I hope!) in terms of the harmonic number $H(m)=\sum_{R=1}^m \frac{1}{R}$:
$$h(x)=H(m)+\frac{1}{m+1}-2$$
However, Mathematica insists on simplifying $h(x)$ to
$$\gamma+\psi_0(m+2)-H(m)-2$$
where $\gamma$ is the Euler-Mascheroni constant and $\psi$ is the polygamma function, which various websites tell me is given by $\psi_0(m+2)=\frac{d^1}{d(m+2)^1} \ln \Gamma(m+2)$.
Is there a proof for Mathematica's simplification? And have I summarised it correctly?
I've never worked with gamma functions, so I'd be grateful for help.
 A: The Digamma Function
Differentiating the logarithm of $\Gamma(x+1)=x\,\Gamma(x)$ gives
$$
\frac{\Gamma'(x+1)}{\Gamma(x+1)}=\frac1x+\frac{\Gamma'(x)}{\Gamma(x)}\tag1
$$
that is, the digamma function satisfies
$$
\psi(x+1)=\frac1x+\psi(x)\tag2
$$
Since $\log(\Gamma(x+1))-\log(\Gamma(x))=\log(x)$, the Mean Value Theorem says that there is a $\xi\in(x,x+1)$ so that
$$
\frac{\Gamma'(\xi)}{\Gamma(\xi)}=\log(x)\tag3
$$
Because $\Gamma$ is log-convex, we get that
$$
\log(x-1)\lt\psi(x)\lt\log(x)\tag4
$$
Therefore,
$$
\lim_{x\to\infty}(\psi(x)-\log(x))=0\tag5
$$

The Extended Harmonic Numbers
If we define
$$
H(x)=\sum_{k=1}^\infty\left(\frac1k-\frac1{k+x}\right)\tag6
$$
then we get by using Telescoping Series
$$
\begin{align}
H(x+1)-H(x)
&=\sum_{k=1}^\infty\left(\frac1{k+x}-\frac1{k+x+1}\right)\\
&=\frac1{x+1}\tag7
\end{align}
$$
and for $n\in\mathbb{N}$, $H(n)=H_n$, the $n^\text{th}$ Harmonic Number.
As is shown in this answer,
$$
\lim_{n\to\infty}(H_n-\log(n))=\gamma\doteq0.57721566490153286060651209\tag8
$$
If $|x-n|\lt1$, then $|H_n-H(x)|\le\frac1{n-1}$ and $|\log(n)-\log(x)|\le\frac1{n-1}$. Therefore, we get that
$$
\lim_{x\to\infty}(H(x)-\log(x))=\gamma\tag9
$$

Relation Between $\boldsymbol{H(x)}$ and $\boldsymbol{\psi(x+1)}$
Equations $(2)$ and $(7)$ say that for $f(x)=H(x)-\psi(x+1)$, we have
$$
f(x+1)=f(x)\tag{10}
$$
Equations $(5)$ and $(9)$ say that
$$
\lim_{x\to\infty}f(x)=\gamma\tag{11}
$$
Therefore, we get that $f(x)=\gamma$; that is,
$$
H(x)=\psi(x+1)+\gamma\tag{12}
$$

To The Question
In your question you seem to be equating $x$ and $m$. I will use only $m$ since it is more commonly used for an integer.
Your equation
$$
h(m)=H_m+\frac1{m+1}-2\tag{13}
$$
is correct. $(13)$ can also be written as
$$
h(m)=H_{m+1}-2\tag{14}
$$
As shown in $(12)$, $H_{m+1}=\psi(m+2)+\gamma$, so $(14)$ becomes
$$
h(m)=\psi(m+2)+\gamma-2\tag{15}
$$
It would seem that the formula that Mathematica gives would be for
$$
\gamma+\psi(m+2)-H_m-2=h(m)-H_m\tag{16}
$$
A:  You can find this in my notes too, in the section about special functions.
The Mittag-Leffler and Weierstrass theorems about the factorization of entire functions equip us with the nice identity
$$ \Gamma(z+1) = e^{-\gamma z}\prod_{n\geq 1}\left(1+\frac{z}{n}\right)^{-1} e^{z/n} \tag{1}$$
which holds uniformly over any compact subset of $\mathbb{C}\setminus\{0,-1,-2,-3,\ldots\}$. It can be proved also through Euler's product formula, which on its turn, if restricted to the positive real line, is a consequence of the definition of the $\Gamma$ function provided by the Bohr-Mollerup theorem, i.e. $\Gamma(x+1)$ is the only log-convex function which fulfills the functional equation $f(x+1)=(x+1)\cdot f(x),\;f(0)=1$, i.e. the "most natural" extension of the factorial function to the positive real numbers. The constant $\gamma$ appearing above stands for
$$ \lim_{n\to +\infty} H_n-\log(n) = \sum_{n\geq 1}\left(\frac{1}{n}-\log\left(1+\frac{1}{n}\right)\right)\stackrel{\mathcal{L}^{-1}}{=}\int_{0}^{+\infty}\frac{1}{e^x-1}-\frac{1}{x e^x}\,dx\approx\frac{1}{\sqrt{3}}. \tag{2}$$
If we apply $\frac{d}{dz}\log(\cdot)$ to both sides of $(1)$ and define the Digamma function as $\frac{d}{dz}\log\Gamma(z)=\frac{\Gamma'(z)}{\Gamma(z)}$ we have:
$$ \psi(z+1)+\gamma = \sum_{n\geq 1}\left(\frac{1}{n}-\frac{1}{n+z}\right) \tag{3} $$
and if $z\in\mathbb{N}$ the RHS of $(3)$ is a telescopic series, equal to $H_z$ by the combinatorial definition of harmonic numbers. Clearly $H_z = \gamma+\psi(z+1)$ allows to define harmonic numbers with non-integer parameters, and the functional identities for the $\Gamma$ function
$$ \Gamma(s)\Gamma(1-s)=\frac{\pi}{\sin(\pi s)},\qquad \Gamma(x+1)=x\cdot\Gamma(x) $$
imply the following functional identities for the Digamma function, by logarithmic differentiation:
$$ \psi(s)-\psi(1-s) = -\pi\cot(\pi s),\qquad \psi(s+1)=\frac{1}{s}+\psi(s). \tag{4}$$
We also have duplication and multiplication formulas, and a little gem from Gauss. They lead to a number of interesting facts, among them:
$$ H_{1/2}=2-2\log 2,\qquad \int_{0}^{+\infty}\frac{dx}{1+x^a}=\frac{\pi}{a\sin\frac{\pi}{a}}\quad\forall a>1. \tag{5}$$
By differentiating $(3)$ and invoking creative telescoping, a short proof of Stirling's approximation follows.
